From 358d36843be4cb87e7fe0655d2f3ec62a8005cd6 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 3 Aug 2021 14:15:05 +0200
Subject: [PATCH] Kernel are now z-dependant
---
src/closure_mod.F90 | 16 ++++++------
src/collision_mod.F90 | 56 +++++++++++++++++++++---------------------
src/compute_Sapj.F90 | 4 +--
src/poisson.F90 | 4 +--
src/processing_mod.F90 | 22 ++++++++---------
5 files changed, 51 insertions(+), 51 deletions(-)
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 4185a1a1..1b4dabc5 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -32,8 +32,8 @@ SUBROUTINE apply_closure_model
moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
- kernel_e(ijsg_e,ikx,iky) = 0._dp
- kernel_e(ijeg_e,ikx,iky) = 0._dp
+ kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
+ kernel_e(ijeg_e,ikx,iky,iz) = 0._dp
DO ip = ipsg_i,ipeg_i
moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
@@ -45,8 +45,8 @@ SUBROUTINE apply_closure_model
moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
- kernel_i(ijsg_i,ikx,iky) = 0._dp
- kernel_i(ijeg_i,ikx,iky) = 0._dp
+ kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
+ kernel_i(ijeg_i,ikx,iky,iz) = 0._dp
ENDDO
ENDDO
ENDDO
@@ -66,8 +66,8 @@ SUBROUTINE apply_closure_model
moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
- kernel_e(ijsg_e,ikx,iky) = 0._dp
- kernel_e(ijeg_e,ikx,iky) = 0._dp
+ kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
+ kernel_e(ijeg_e,ikx,iky,iz) = 0._dp
DO ip = ipsg_i,ipeg_i
moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
@@ -79,8 +79,8 @@ SUBROUTINE apply_closure_model
moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
- kernel_i(ijsg_i,ikx,iky) = 0._dp
- kernel_i(ijeg_i,ikx,iky) = 0._dp
+ kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
+ kernel_i(ijeg_i,ikx,iky,iz) = 0._dp
ENDDO
ENDDO
ENDDO
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index b6dca646..5fe8fa39 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -43,7 +43,7 @@ CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
! Get adiabatic moment
- TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_)*phi(ikx_,iky_,iz_)
+ TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_,iz_)*phi(ikx_,iky_,iz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -51,9 +51,9 @@ CONTAINS
DO in_ = 1,jmaxe+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_e(in_,ikx_,iky_)
- Knp1 = Kernel_e(in_+1,ikx_,iky_)
- Knm1 = Kernel_e(in_-1,ikx_,iky_)
+ Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_)
+ Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_)
+ Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_e(1,in_ ,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -67,10 +67,10 @@ CONTAINS
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
! Add energy restoring term
- TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,ikx_,iky_)
- TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_)
- TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,ikx_,iky_)
- TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_) - Kernel_e(ij_-1,ikx_,iky_))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,ikx_,iky_,iz_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_,iz_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,ikx_,iky_,iz_)
+ TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_,iz_) - Kernel_e(ij_-1,ikx_,iky_,iz_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -79,9 +79,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxe+1
! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,ikx_,iky_) * moments_e(2,in_,ikx_,iky_,iz_,updatetlevel)
+ upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_) * moments_e(2,in_,ikx_,iky_,iz_,updatetlevel)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_)
+ TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -93,9 +93,9 @@ CONTAINS
DO in_ = 1,jmaxe+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_e(in_ ,ikx_,iky_)
- Knp1 = Kernel_e(in_+1,ikx_,iky_)
- Knm1 = Kernel_e(in_-1,ikx_,iky_)
+ Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_)
+ Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_)
+ Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_e(1,in_,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -106,7 +106,7 @@ CONTAINS
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
- TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
@@ -149,7 +149,7 @@ CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
! Get adiabatic moment
- TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_)*phi(ikx_,iky_,iz_)
+ TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_,iz_)*phi(ikx_,iky_,iz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -157,9 +157,9 @@ CONTAINS
DO in_ = 1,jmaxi+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_i(in_,ikx_,iky_)
- Knp1 = Kernel_i(in_+1,ikx_,iky_)
- Knm1 = Kernel_i(in_-1,ikx_,iky_)
+ Knp0 = Kernel_i(in_ ,ikx_,iky_,iz_)
+ Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_)
+ Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_i(1,in_ ,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -173,10 +173,10 @@ CONTAINS
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
! Add energy restoring term
- TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,ikx_,iky_)
- TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_)
- TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,ikx_,iky_)
- TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_) - Kernel_i(ij_-1,ikx_,iky_))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,ikx_,iky_,iz_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,ikx_,iky_,iz_)
+ TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_,iz_) - Kernel_i(ij_-1,ikx_,iky_,iz_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -185,9 +185,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxi+1
! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,ikx_,iky_) * moments_i(2,in_,ikx_,iky_,iz_,updatetlevel)
+ upar_ = upar_ + Kernel_i(in_,ikx_,iky_,iz_) * moments_i(2,in_,ikx_,iky_,iz_,updatetlevel)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_)
+ TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -199,9 +199,9 @@ CONTAINS
DO in_ = 1,jmaxi+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_i(in_ ,ikx_,iky_)
- Knp1 = Kernel_i(in_+1,ikx_,iky_)
- Knm1 = Kernel_i(in_-1,ikx_,iky_)
+ Knp0 = Kernel_i(in_ ,ikx_,iky_,iz_)
+ Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_)
+ Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_i(1,in_,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -212,7 +212,7 @@ CONTAINS
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
- TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index e66080a7..107e3349 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -54,7 +54,7 @@ zloop: DO iz = izs,ize
kyloope: DO iky = ikys,ikye ! Loop over ky
kx = kxarray(ikx)
ky = kyarray(iky)
- kerneln = kernel_e(in, ikx, iky)
+ kerneln = kernel_e(in, ikx, iky, iz)
! First convolution terms
Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
@@ -145,7 +145,7 @@ zloop: DO iz = izs,ize
kyloopi: DO iky = ikys,ikye ! Loop over ky
kx = kxarray(ikx)
ky = kyarray(iky)
- kerneln = kernel_i(in, ikx, iky)
+ kerneln = kernel_i(in, ikx, iky, iz)
! First convolution terms
Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
diff --git a/src/poisson.F90 b/src/poisson.F90
index 6ec3a2fd..00ad9de7 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -37,7 +37,7 @@ SUBROUTINE poisson
sum_kernel2_e = 0._dp
! loop over n only if the max polynomial degree
DO ine=1,jmaxe+1 ! ine = n+1
- Kne = kernel_e(ine, ikx, iky)
+ Kne = kernel_e(ine,ikx,iky,iz)
sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikx, iky, iz, updatetlevel)
sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
END DO
@@ -47,7 +47,7 @@ SUBROUTINE poisson
sum_kernel2_i = 0._dp
! loop over n only if the max polynomial degree
DO ini=1,jmaxi+1
- Kni = kernel_i(ini, ikx, iky)
+ Kni = kernel_i(ini,ikx,iky,iz)
sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikx, iky, iz, updatetlevel)
sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ...
END DO
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 129a4c17..b3e5dd94 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -34,7 +34,7 @@ SUBROUTINE compute_radial_ion_transport
imagu * ky_ * moments_i(1,1,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
DO ij = ijs_i, ije_i
pflux_local = pflux_local - &
- imagu * ky_ * kernel_i(ij,ikx,iky) * moments_i(1,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
+ imagu * ky_ * kernel_i(ij,ikx,iky,iz) * moments_i(1,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
ENDDO
ENDDO
ENDDO
@@ -81,14 +81,14 @@ SUBROUTINE compute_density
! electron density
dens_e(ikx,iky,iz) = 0._dp
DO ij = ijs_e, ije_e
- dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky) * &
- (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky)*phi(ikx,iky,iz))
+ dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky,iz) * &
+ (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz))
ENDDO
! ion density
dens_i(ikx,iky,iz) = 0._dp
DO ij = ijs_i, ije_i
- dens_i(ikx,iky,iz) = dens_i(ikx,iky,iz) + kernel_i(ij,ikx,iky) * &
- (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky)*phi(ikx,iky,iz))
+ dens_i(ikx,iky,iz) = dens_i(ikx,iky,iz) + kernel_i(ij,ikx,iky,iz) * &
+ (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iz))
ENDDO
ENDDO
ENDDO
@@ -118,9 +118,9 @@ SUBROUTINE compute_temperature
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
temp_e(ikx,iky,iz) = temp_e(ikx,iky,iz) + &
- 2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky) - (j_dp+1)*kernel_e(ij+1,ikx,iky) - j_dp*kernel_e(ij-1,ikx,iky))&
- * (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky)*phi(ikx,iky,iz)) + &
- SQRT2/3._dp * kernel_e(ij,ikx,iky) * moments_e(3,ij,ikx,iky,iz,updatetlevel)
+ 2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky,iz) - (j_dp+1)*kernel_e(ij+1,ikx,iky,iz) - j_dp*kernel_e(ij-1,ikx,iky,iz))&
+ * (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz)) &
+ + SQRT2/3._dp * kernel_e(ij,ikx,iky,iz) * moments_e(3,ij,ikx,iky,iz,updatetlevel)
ENDDO
! ion temperature
@@ -128,9 +128,9 @@ SUBROUTINE compute_temperature
DO ij = ijs_i, ije_i
j_dp = REAL(ij-1,dp)
temp_i(ikx,iky,iz) = temp_i(ikx,iky,iz) + &
- 2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky) - (j_dp+1)*kernel_i(ij+1,ikx,iky) - j_dp*kernel_i(ij-1,ikx,iky))&
- * (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky)*phi(ikx,iky,iz)) + &
- SQRT2/3._dp * kernel_i(ij,ikx,iky) * moments_i(3,ij,ikx,iky,iz,updatetlevel)
+ 2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz) - j_dp*kernel_i(ij-1,ikx,iky,iz))&
+ * (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iz)) &
+ + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz) * moments_i(3,ij,ikx,iky,iz,updatetlevel)
ENDDO
ENDDO
ENDDO
--
GitLab