diff --git a/Makefile b/Makefile
index c405b8e434025e986bb7984516a8f806f072c3d7..fa7aaed2a01db980a7a5a296d9d90d6bec65e834 100644
--- a/Makefile
+++ b/Makefile
@@ -72,6 +72,9 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@
+ $(OBJDIR)/calculus_mod.o : src/calculus_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/calculus_mod.F90 -o $@
+
$(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
@@ -102,6 +105,9 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
+ $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
+
$(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/ppinit.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
@@ -114,9 +120,6 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@
- $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
-
$(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
@@ -126,7 +129,7 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
- $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@
$(OBJDIR)/numerics_mod.o : src/numerics_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/utility_mod.o
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 668a95bcca7c4fa0b3467c4496db5de57af75eb2..b63b4ec08825faea0b5d9787b628171dcc7477b1 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -39,27 +39,23 @@ MODULE array
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
! Geoemtrical operators
! Curvature
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
! Jacobian
- REAL(dp), DIMENSION(:), ALLOCATABLE :: Jacobian ! dimensions: z
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
! Metric
- REAL(dp), DIMENSION(:), ALLOCATABLE :: gxx, gxy, gyy, gxz, gyz
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gyy, gxz, gyz ! dimensions: z, odd/even p
! derivatives of magnetic field strength
- REAL(dp), DIMENSION(:), allocatable :: gradzB ! dimensions: z
- REAL(dp), DIMENSION(:), allocatable :: gradxB
+ REAL(dp), DIMENSION(:,:), allocatable :: gradzB ! dimensions: z, odd/even p
+ REAL(dp), DIMENSION(:,:), allocatable :: gradxB
! Relative magnetic field strength
- REAL(dp), DIMENSION(:), allocatable :: hatB
+ REAL(dp), DIMENSION(:,:), allocatable :: hatB
! Relative strength of major radius
- REAL(dp), DIMENSION(:), allocatable :: hatR
- ! Geometrical factors
- REAL(dp), DIMENSION(:), allocatable :: Gamma1
- REAL(dp), DIMENSION(:), allocatable :: Gamma2
- REAL(dp), DIMENSION(:), allocatable :: Gamma3
+ REAL(dp), DIMENSION(:,:), allocatable :: hatR
! Some geometrical coefficients
- REAL(dp), DIMENSION(:) , allocatable :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
- ! Kernel function evaluation (ij,ikx,iky,iz)
- REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: kernel_e
- REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: kernel_i
+ REAL(dp), DIMENSION(:,:) , allocatable :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
+ ! Kernel function evaluation (ij,ikx,iky,iz,odd/even p)
+ REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_e
+ REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i
! Poisson operator (ikx,iky,iz)
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
diff --git a/src/auxval.F90 b/src/auxval.F90
index 8e642f6f8b7e744540429c938332ace97c5aedd3..c177f1f85456bb448b3b6f818562a289854c6c3e 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -29,6 +29,7 @@ subroutine auxval
CALL set_kygrid ! azymuthal modes
CALL set_zgrid ! field aligned angle
+ IF ((my_id .EQ. 0) .AND. SG) WRITE(*,*) '--2 staggered z grids--'
CALL memory ! Allocate memory for global arrays
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..d4ae848add5ff015324d0fa0af940a1ee3cc5686
--- /dev/null
+++ b/src/calculus_mod.F90
@@ -0,0 +1,158 @@
+MODULE calculus
+ ! Routine to evaluate gradients, interpolation schemes and integrals
+ USE basic
+ USE prec_const
+ USE grid
+ IMPLICIT NONE
+ REAL(dp), dimension(-2:2) :: dz_usu = &
+ (/ onetwelfth, -twothird, &
+ 0._dp, &
+ twothird, -onetwelfth /) ! fd4 centered stencil
+ REAL(dp), dimension(-2:1) :: dz_o2e = &
+ (/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
+ REAL(dp), dimension(-1:2) :: dz_e2o = &
+ (/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
+
+ PUBLIC :: simpson_rule_z, interp_z, grad_z
+
+CONTAINS
+
+SUBROUTINE grad_z(target,f,ddzf)
+ implicit none
+ ! Compute the periodic boundary condition 4 points centered finite differences
+ ! formula among staggered grid or not.
+ ! not staggered : the derivative results must be on the same grid as the field
+ ! staggered : the derivative is computed from a grid to the other
+ INTEGER, INTENT(IN) :: target
+ COMPLEX(dp),dimension( izs:ize ), intent(in) :: f
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzf
+ IF(SG) THEN
+ IF(TARGET .EQ. 0) THEN
+ CALL grad_z_o2e(f,ddzf)
+ ELSE
+ CALL grad_z_e2o(f,ddzf)
+ ENDIF
+ ELSE ! No staggered grid -> usual centered finite differences
+ DO iz = 3,Nz-2
+ ddzf(iz) = dz_usu(-2)*f(iz-2) + dz_usu(-1)*f(iz-1) &
+ +dz_usu( 1)*f(iz+1) + dz_usu( 2)*f(iz+2)
+ ENDDO
+ ! Periodic boundary conditions
+ ddzf(Nz-1) = dz_usu(-2)*f(Nz-3) + dz_usu(-1)*f(Nz-2) &
+ +dz_usu(+1)*f(Nz ) + dz_usu(+2)*f(1 )
+ ddzf(Nz ) = dz_usu(-2)*f(Nz-2) + dz_usu(-1)*f(Nz-1) &
+ +dz_usu(+1)*f(1 ) + dz_usu(+2)*f(2 )
+ ddzf(1 ) = dz_usu(-2)*f(Nz-1) + dz_usu(-1)*f(Nz ) &
+ +dz_usu(+1)*f(2 ) + dz_usu(+2)*f(3 )
+ ddzf(2 ) = dz_usu(-2)*f(Nz ) + dz_usu(-1)*f(1 ) &
+ +dz_usu(+1)*f(3 ) + dz_usu(+2)*f(4)
+ ENDIF
+END SUBROUTINE grad_z
+
+SUBROUTINE grad_z_o2e(fo,dzfe) ! Paruta 2018 eq (27)
+ ! gives the value of a field from the odd grid to the even one
+ implicit none
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: dzfe !
+ DO iz = 3,Nz-1
+ dzfe(iz) = dz_o2e(-2)*fo(iz-2) + dz_o2e(-1)*fo(iz-1) &
+ +dz_o2e( 0)*fo(iz ) + dz_o2e( 1)*fo(iz+1)
+ ENDDO
+ ! Periodic boundary conditions
+ dzfe(Nz) = dz_o2e(-2)*fo(Nz-2) + dz_o2e(-1)*fo(Nz-1) &
+ +dz_o2e( 0)*fo(Nz ) + dz_o2e( 1)*fo(1)
+ dzfe(1 ) = dz_o2e(-2)*fo(Nz-1) + dz_o2e(-1)*fo(Nz ) &
+ +dz_o2e( 0)*fo(1 ) + dz_o2e( 1)*fo(2)
+ dzfe(2 ) = dz_o2e(-2)*fo(Nz ) + dz_o2e(-1)*fo(1 ) &
+ +dz_o2e( 0)*fo(2 ) + dz_o2e( 1)*fo(3)
+END SUBROUTINE grad_z_o2e
+
+SUBROUTINE grad_z_e2o(fe,dzfo)
+ ! gives the value of a field from the even grid to the odd one
+ implicit none
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: dzfo
+ DO iz = 2,Nz-2
+ dzfo(iz) = dz_e2o(-1)*fe(iz-1) + dz_e2o(0)*fe(iz ) &
+ +dz_e2o( 1)*fe(iz+1) + dz_e2o(2)*fe(iz+2)
+ ENDDO
+ ! Periodic boundary conditions
+ dzfo(Nz-1) = dz_e2o(-1)*fe(Nz-2) + dz_e2o(0)*fe(Nz-1) &
+ +dz_e2o( 1)*fe(Nz ) + dz_e2o(2)*fe(1 )
+ dzfo(Nz ) = dz_e2o(-1)*fe(Nz-1) + dz_e2o(0)*fe(Nz ) &
+ +dz_e2o( 1)*fe(1 ) + dz_e2o(2)*fe(2 )
+ dzfo(1 ) = dz_e2o(-1)*fe(Nz ) + dz_e2o(0)*fe(1 ) &
+ +dz_e2o( 1)*fe(2 ) + dz_e2o(2)*fe(3 )
+END SUBROUTINE grad_z_e2o
+
+
+SUBROUTINE interp_z(target,f_in,f_out)
+ ! Function meant to interpolate one field defined on a even/odd z into
+ ! the other odd/even z grid.
+ ! If Staggered Grid flag (SG) is false, returns identity
+ implicit none
+ INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: f_in
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: f_out !
+ IF(SG) THEN
+ IF(TARGET .EQ. 0) THEN
+ CALL interp_o2e_z(f_in,f_out)
+ ELSE
+ CALL interp_e2o_z(f_in,f_out)
+ ENDIF
+ ELSE ! No staggered grid -> identity
+ f_out(:) = f_in(:)
+ ENDIF
+END SUBROUTINE interp_z
+
+SUBROUTINE interp_o2e_z(fo,fe)
+ ! gives the value of a field from the odd grid to the even one
+ implicit none
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: fe !
+ DO iz = 2,Nz
+ fe(iz) = 0.5_dp*(fo(iz)+fo(iz-1))
+ ENDDO
+ ! Periodic boundary conditions
+ fe(1) = 0.5_dp*(fo(1) + fo(Nz))
+END SUBROUTINE interp_o2e_z
+
+SUBROUTINE interp_e2o_z(fe,fo)
+ ! gives the value of a field from the even grid to the odd one
+ implicit none
+ COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
+ COMPLEX(dp),dimension(1:Nz), intent(out) :: fo
+ DO iz = 1,Nz-1
+ fo(iz) = 0.5_dp*(fe(iz+1)+fe(iz))
+ ENDDO
+ ! Periodic boundary conditions
+ fo(Nz) = 0.5_dp*(fe(1) + fe(Nz))
+END SUBROUTINE interp_e2o_z
+
+
+SUBROUTINE simpson_rule_z(f,intf)
+ ! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
+ !from molix BJ Frei
+ implicit none
+ complex(dp),dimension(izs:ize), intent(in) :: f
+ COMPLEX(dp), intent(out) :: intf
+ COMPLEX(dp) :: buff_
+ IF(Nz .GT. 1) THEN
+ IF(mod(Nz,2) .ne. 0 ) THEN
+ ERROR STOP 'Simpson rule: Nz must be an even number !!!!'
+ ENDIF
+ buff_ = 0._dp
+ DO iz = izs, Nz/2
+ IF(iz .eq. Nz/2) THEN ! ... iz = ize
+ buff_ = buff_ + (f(izs) + 4._dp*f(ize) + f(ize-1 ))
+ ELSE
+ buff_ = buff_ + (f(2*iz+1) + 4._dp*f(2*iz) + f(2*iz-1 ))
+ ENDIF
+ ENDDO
+ intf = buff_*deltaz/3._dp
+ ELSE
+ intf = f(izs)
+ ENDIF
+END SUBROUTINE simpson_rule_z
+
+END MODULE calculus
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 4b124fc32c7c2abfe9e397e69d27566bd26d9012..b1a5a7efad986d3f0f499e160304c3288c2b7fdf 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -80,8 +80,8 @@ SUBROUTINE ghosts_truncation
moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDIF
ENDDO
- kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
- kernel_e(ijeg_e,ikx,iky,iz) = 0._dp
+ kernel_e(ijsg_e,ikx,iky,iz,:) = 0._dp
+ kernel_e(ijeg_e,ikx,iky,iz,:) = 0._dp
ENDIF
DO ip = ipsg_i,ipeg_i
moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
@@ -95,8 +95,8 @@ SUBROUTINE ghosts_truncation
moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDIF
ENDDO
- kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
- kernel_i(ijeg_i,ikx,iky,iz) = 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 3d2bd674d1358deaf2d320569451d7bfa0430364..e81b6ed3d9cf085b3fccebc6e36b709ca43ae760 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -79,13 +79,14 @@ CONTAINS
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
COMPLEX(dp) :: nadiab_moment_0j
REAL(dp) :: Knp0, Knp1, Knm1, kp
- INTEGER :: in_
+ INTEGER :: in_, eo_
REAL(dp) :: n_dp, j_dp, p_dp, be_, be_2
!** Auxiliary variables **
p_dp = REAL(parray_e(ip_),dp)
+ eo_ = MODULO(parray_e(ip_),2)
j_dp = REAL(jarray_e(ij_),dp)
- kp = kparray(ikx_,iky_,iz_)
+ kp = kparray(ikx_,iky_,iz_,eo_)
be_2 = kp**2 * sigmae2_taue_o2 ! this is (be/2)^2
be_ = 2_dp*kp * sqrt_sigmae2_taue_o2 ! this is be
@@ -97,7 +98,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_,iz_)*phi(ikx_,iky_,iz_)
+ TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_,iz_,eo_)*phi(ikx_,iky_,iz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -105,9 +106,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_,iz_)
- Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_)
- Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_)
+ Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_,eo_)
+ Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_)
+ Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_e(ip0_e,in_ ,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -121,10 +122,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_,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_))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_,iz_,eo_) - Kernel_e(ij_-1,ikx_,iky_,iz_,eo_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -133,9 +134,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxe+1
! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_) * moments_e(ip1_e,in_,ikx_,iky_,iz_,updatetlevel)
+ upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_,eo_) * moments_e(ip1_e,in_,ikx_,iky_,iz_,updatetlevel)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_)
+ TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -147,9 +148,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_,iz_)
- Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_)
- Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_)
+ Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_,eo_)
+ Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_)
+ Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_e(ip0_e,in_,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -160,7 +161,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_,iz_)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
@@ -186,13 +187,14 @@ CONTAINS
COMPLEX(dp) :: bi_, bi_2
COMPLEX(dp) :: nadiab_moment_0j
REAL(dp) :: Knp0, Knp1, Knm1, kp
- INTEGER :: in_
+ INTEGER :: in_, eo_
REAL(dp) :: n_dp, j_dp, p_dp
!** Auxiliary variables **
p_dp = REAL(parray_i(ip_),dp)
+ eo_ = MODULO(parray_i(ip_),2)
j_dp = REAL(jarray_i(ij_),dp)
- kp = kparray(ikx_,iky_,iz_)
+ kp = kparray(ikx_,iky_,iz_,eo_)
bi_2 = kp**2 *sigmai2_taui_o2 ! this is (bi/2)^2
bi_ = 2_dp*kp*sqrt_sigmai2_taui_o2 ! this is bi
@@ -204,7 +206,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_,iz_)*phi(ikx_,iky_,iz_)
+ TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_,iz_,eo_)*phi(ikx_,iky_,iz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -212,9 +214,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_,iz_)
- Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_)
- Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_)
+ Knp0 = Kernel_i(in_ ,ikx_,iky_,iz_,eo_)
+ Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_)
+ Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_i(ip0_i,in_ ,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -228,10 +230,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_,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_))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_,iz_,eo_) - Kernel_i(ij_-1,ikx_,iky_,iz_,eo_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -240,9 +242,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxi+1
! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,ikx_,iky_,iz_) * moments_i(ip1_i,in_,ikx_,iky_,iz_,updatetlevel)
+ upar_ = upar_ + Kernel_i(in_,ikx_,iky_,iz_,eo_) * moments_i(ip1_i,in_,ikx_,iky_,iz_,updatetlevel)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_)
+ TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -254,9 +256,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_,iz_)
- Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_)
- Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_)
+ Knp0 = Kernel_i(in_ ,ikx_,iky_,iz_,eo_)
+ Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_)
+ Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
nadiab_moment_0j = moments_i(ip0_i,in_,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
! Density
@@ -267,7 +269,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_,iz_)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
@@ -313,11 +315,11 @@ CONTAINS
ENDDO
IF (num_procs_p .GT. 1) THEN
! Sum up all the sub collision terms on root 0
- CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 1003, comm_p, ierr)
+ CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
! distribute the sum over the process among p
CALL MPI_SCATTERV(buffer_e, counts_np_e, displs_np_e, MPI_DOUBLE_COMPLEX,&
TColl_distr_e, local_np_e, MPI_DOUBLE_COMPLEX,&
- 1003, comm_p, ierr)
+ 0, comm_p, ierr)
ELSE
TColl_distr_e = local_sum_e
ENDIF
@@ -335,12 +337,12 @@ CONTAINS
ENDDO
IF (num_procs_p .GT. 1) THEN
! Reduce the local_sums to root = 0
- CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 1004, comm_p, ierr)
+ CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
! buffer contains the entire collision term along p, we scatter it between
! the other processes (use of scatterv since Pmax/Np is not an integer)
CALL MPI_SCATTERV(buffer_i, counts_np_i, displs_np_i, MPI_DOUBLE_COMPLEX,&
TColl_distr_i, local_np_i, MPI_DOUBLE_COMPLEX, &
- 1004, comm_p, ierr)
+ 0, comm_p, ierr)
ELSE
TColl_distr_i = local_sum_i
ENDIF
@@ -694,7 +696,7 @@ CONTAINS
IF (my_id .EQ. 0 ) WRITE(*,*) '...Interpolation from matrices kperp to simulation kx,ky...'
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
- kperp_sim = SQRT(kxarray(ikx)**2+kyarray(iky)**2) ! current simulation kperp
+ kperp_sim = kparray(ikx,iky,1,0) ! current simulation kperp
! Find the interval in kp grid mat where kperp_sim is contained
! Loop over the whole kp mat grid to find the smallest kperp that is
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index bfc834960b7c5604d0729a9a6df3e0c3e20ca0e7..77c12b8a19b1dceb8ea5422d55b2e27950d345b8 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -29,20 +29,10 @@ SUBROUTINE compute_Sapj
IF(KIN_E) THEN
zloope: DO iz = izs,ize
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
- p_int = parray_e(ip)
+ eo = MODULO(parray_e(ip),2)
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
- j_int=jarray_e(ij)
- ! GF closure check (spare computations too)
- GF_CLOSURE_e: IF ((CLOS.EQ.1) .AND. (p_int+2*j_int .GT. dmaxe)) THEN
- ! Do nothing
- DO ikx = ikxs, ikxe
- DO iky = ikys, ikye
- Sepj(ip,ij,ikx,iky,iz) = 0._dp
- ENDDO
- ENDDO
- ELSE
+ j_int=jarray_e(ij)
real_data_c = 0._dp ! initialize sum over real nonlinear term
-
! Set non linear sum truncation
IF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxe
@@ -58,7 +48,7 @@ SUBROUTINE compute_Sapj
kyloope: DO iky = ikys,ikye ! Loop over ky
kx = kxarray(ikx)
ky = kyarray(iky)
- kerneln = kernel_e(in, ikx, iky, iz)
+ kerneln = kernel_e(in, ikx, iky, iz, eo)
! First convolution terms
Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
@@ -116,7 +106,6 @@ SUBROUTINE compute_Sapj
Sepj(ip,ij,ikx,iky,iz) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
- ENDIF GF_CLOSURE_e
ENDDO jloope
ENDDO ploope
ENDDO zloope
@@ -126,22 +115,12 @@ ENDIF
!!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
zloopi: DO iz = izs,ize
ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
-
+ eo = MODULO(parray_i(ip),2)
! we check if poly degree is even (eq to index is odd) to spare computation
!EVEN_P_i: IF (.TRUE. .OR. (MODULO(ip,2) .EQ. 1) .OR. (.NOT. COMPUTE_ONLY_EVEN_P)) THEN
jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
- j_int=jarray_i(ij)
- ! GF closure check (spare computations too)
- GF_CLOSURE_i: IF ((CLOS.EQ.1) .AND. (p_int+2*j_int .GT. dmaxi)) THEN
- ! Do nothing
- DO ikx = ikxs, ikxe
- DO iky = ikys, ikye
- Sipj(ip,ij,ikx,iky,iz) = 0._dp
- ENDDO
- ENDDO
- ELSE
+ j_int=jarray_i(ij)
real_data_c = 0._dp ! initialize sum over real nonlinear term
-
! Set non linear sum truncation
IF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxi
@@ -157,7 +136,7 @@ zloopi: DO iz = izs,ize
kyloopi: DO iky = ikys,ikye ! Loop over ky
kx = kxarray(ikx)
ky = kyarray(iky)
- kerneln = kernel_i(in, ikx, iky, iz)
+ kerneln = kernel_i(in, ikx, iky, iz, eo)
! First convolution terms
Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
@@ -215,7 +194,6 @@ zloopi: DO iz = izs,ize
Sipj(ip,ij,ikx,iky,iz) = cmpx_data_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
- ENDIF GF_CLOSURE_i
ENDDO jloopi
ENDDO ploopi
ENDDO zloopi
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 2ef6d63d89deb20184000084cd792aa9a8d09b75..b874797220b86f2e83a59596c503e683473bf03c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -8,7 +8,7 @@ module geometry
use array
use fields
use basic
- use utility, ONLY: simpson_rule_z
+ use calculus, ONLY: simpson_rule_z
implicit none
public
@@ -33,20 +33,22 @@ contains
! Evaluate perpendicular wavenumber
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
+ DO eo = 0,1
DO iz = izs,ize
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx = ikxs, ikxe
kx = kxarray(ikx)
- kparray(ikx, iky, iz) = &
- SQRT( gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2)/hatB(iz)
+ kparray(ikx, iky, iz, eo) = &
+ SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
! there is a factor 1/B from the normalization; important to match GENE
ENDDO
ENDDO
ENDDO
+ ENDDO
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
- integrant = Jacobian ! Convert into complex array
+ integrant = Jacobian(:,0) ! Convert into complex array
CALL simpson_rule_z(integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
@@ -59,39 +61,41 @@ contains
implicit none
REAL(dp) :: z, kx, ky
+ parity: DO eo = 0,1
zloop: DO iz = izs,ize
- z = zarray(iz)
+ z = zarray(iz,eo)
! metric
- gxx(iz) = 1._dp
- gxy(iz) = shear*z
- gyy(iz) = 1._dp + (shear*z)**2
+ gxx(iz,eo) = 1._dp
+ gxy(iz,eo) = shear*z
+ gyy(iz,eo) = 1._dp + (shear*z)**2
! Relative strengh of radius
- hatR(iz) = 1._dp + eps*COS(z)
+ hatR(iz,eo) = 1._dp + eps*COS(z)
! Jacobian
- Jacobian(iz) = q0*hatR(iz)
+ Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
- hatB(iz) = 1._dp / hatR(iz)
+ hatB(iz,eo) = 1._dp / hatR(iz,eo)
! Derivative of the magnetic field strenght
- gradxB(iz) = -COS(z)
- gradzB(iz) = eps * SIN(z) / hatR(iz)
+ gradxB(iz,eo) = -COS(z)
+ gradzB(iz,eo) = eps * SIN(z) / hatR(iz,eo)
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(ikx, iky, iz) = (-SIN(z)*kx - (COS(z) + shear* z* SIN(z))*ky) * hatB(iz) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(ikx, iky, iz,eo) = (-SIN(z)*kx - (COS(z) + shear* z* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
- gradz_coeff(iz) = 1._dp / Jacobian(iz) / hatB(iz)
+ gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
+ ENDDO parity
END SUBROUTINE eval_salphaB_geometry
!
@@ -103,34 +107,36 @@ contains
implicit none
REAL(dp) :: z, kx, ky
+ parity: DO eo = 0,1
zloop: DO iz = izs,ize
- z = zarray(iz)
+ z = zarray(iz,eo)
! metric
- gxx(iz) = 1._dp
- gxy(iz) = 0._dp
- gyy(iz) = 1._dp
+ gxx(iz,eo) = 1._dp
+ gxy(iz,eo) = 0._dp
+ gyy(iz,eo) = 1._dp
! Relative strengh of radius
- hatR(iz) = 1._dp
+ hatR(iz,eo) = 1._dp
! Jacobian
- Jacobian(iz) = 1._dp
+ Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
- hatB(iz) = 1._dp
+ hatB(iz,eo) = 1._dp
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(ikx, iky, iz) = -kx ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(ikx, iky, iz,eo) = -kx ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
- gradz_coeff(iz) = 1._dp
- ENDDO zloop
+ gradz_coeff(iz,eo) = 1._dp
+ ENDDO zloop
+ ENDDO parity
END SUBROUTINE eval_1D_geometry
!
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index a2736db8c102e24b5204804c3cd7216f61cc8157..0cf0f39471eda127acdb41f3492dcf03e7566060 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -37,13 +37,16 @@ MODULE grid
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_y
! Grids containing position in physical space
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: xarray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: yarray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray, zarray_full
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: izarray
- REAL(dp), PUBLIC, PROTECTED :: deltax, deltay, deltaz, inv_deltaz
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: xarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: yarray
+ ! Local and global z grids, 2D since it has to store odd and even grids
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: zarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray_full
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: izarray
+ REAL(dp), PUBLIC, PROTECTED :: deltax, deltay, deltaz, inv_deltaz
INTEGER, PUBLIC, PROTECTED :: ixs, ixe, iys, iye, izs, ize
INTEGER, PUBLIC, PROTECTED :: izgs, izge ! ghosts
+ LOGICAL, PUBLIC, PROTECTED :: SG = .true.! shifted grid flag
INTEGER, PUBLIC :: ir,iz ! counters
integer(C_INTPTR_T), PUBLIC :: local_nkx, local_nky
integer(C_INTPTR_T), PUBLIC :: local_nkx_offset, local_nky_offset
@@ -57,12 +60,13 @@ MODULE grid
! Grids containing position in fourier space
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kxarray, kxarray_full
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kyarray, kyarray_full
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, PUBLIC :: kparray
+ ! Kperp array depends on kx, ky, z (geometry), eo (even or odd zgrid)
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC :: kparray
REAL(dp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max!, kp_max
REAL(dp), PUBLIC, PROTECTED :: local_kxmax, local_kymax
INTEGER, PUBLIC, PROTECTED :: ikxs, ikxe, ikys, ikye!, ikps, ikpe
INTEGER, PUBLIC, PROTECTED :: ikx_0, iky_0, ikx_max, iky_max ! Indices of k-grid origin and max
- INTEGER, PUBLIC :: ikx, iky, ip, ij, ikp, pp2 ! counters
+ INTEGER, PUBLIC :: ikx, iky, ip, ij, ikp, pp2, eo ! counters
LOGICAL, PUBLIC, PROTECTED :: contains_kx0 = .false. ! flag if the proc contains kx=0 index
LOGICAL, PUBLIC, PROTECTED :: contains_ky0 = .false. ! flag if the proc contains ky=0 index
LOGICAL, PUBLIC, PROTECTED :: contains_kxmax = .false. ! flag if the proc contains kx=kxmax index
@@ -102,7 +106,7 @@ CONTAINS
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- Nx, Lx, Ny, Ly, Nz, q0, shear, eps
+ Nx, Lx, Ny, Ly, Nz, q0, shear, eps, SG
READ(lu_in,grid)
!! Compute the maximal degree of full GF moments set
@@ -356,19 +360,35 @@ CONTAINS
USE prec_const
IMPLICIT NONE
INTEGER :: i_, ngz
+ REAL :: grids_shift
! Start and END indices of grid
- izs = 1
- ize = Nz
- ALLOCATE(zarray(izs:ize))
+ izs = 1
+ ize = Nz
+ ALLOCATE(zarray(izs:ize,0:1))
IF (Nz .EQ. 1) THEN ! full perp case
- deltaz = 1._dp
- zarray(1) = 0
- ELSE
- deltaz = 2._dp*PI/REAL(Nz,dp)
+ deltaz = 1._dp
+ zarray(1,0) = 0._dp
+ zarray(1,1) = 0._dp
+ SG = .false. ! unique perp plane at z=0
+ izgs = izs
+ izge = ize
+ ELSE ! Parallel dimension exists
+ deltaz = 2._dp*PI/REAL(Nz,dp)
inv_deltaz = 1._dp/deltaz
+ IF(SG) THEN ! Shifted grids option
+ grids_shift = deltaz/2._dp ! we shift both z grid
+ ELSE
+ grids_shift = 0._dp
+ ENDIF
DO iz = izs,ize
- zarray(iz) = REAL((iz-1),dp)*deltaz - PI
+ ! Even z grid (Napj with p even and phi)
+ zarray(iz,0) = REAL((iz-1),dp)*deltaz - PI
+ ! Odd z grid (Napj with p odd)
+ zarray(iz,1) = REAL((iz-1),dp)*deltaz - PI + grids_shift
ENDDO
+ ! 2 ghosts cell for four point stencil
+ izgs = izs-2
+ izge = ize+2
ENDIF
if(my_id.EQ.0) write(*,*) '#parallel planes = ', Nz
! Build the full grids on process 0 to diagnose it without comm
@@ -421,6 +441,7 @@ CONTAINS
CALL attach(fidres, TRIM(str), "Lkx", Lkx)
CALL attach(fidres, TRIM(str), "Nky", Nky)
CALL attach(fidres, TRIM(str), "Lky", Lky)
+ CALL attach(fidres, TRIM(str), "SG", SG)
END SUBROUTINE grid_outputinputs
FUNCTION bare(p_,j_)
diff --git a/src/inital.F90 b/src/inital.F90
index 0627b10e52124356cc6bf32e09299a9548ddd6b7..c7e1db9de7d71e996282e7ac7758bb56880486bd 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -286,9 +286,8 @@ SUBROUTINE init_phi
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz=izs,ize
- kp = kparray(ikx,iky,iz)
CALL RANDOM_NUMBER(noise)
- phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))*EXP(-0.1*kp**2)!*AA_x(ikx)*AA_y(iky)
+ phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
ENDDO
END DO
END DO
diff --git a/src/memory.F90 b/src/memory.F90
index 0cfc3203f5d87ed3d4fc474e0ca2afd8a09509d9..79249e7572f80faf763f5fe76210a569a27836a7 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -20,7 +20,7 @@ SUBROUTINE memory
CALL allocate_array( Ne00, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( dens_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( temp_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
@@ -46,7 +46,7 @@ SUBROUTINE memory
CALL allocate_array(Ni00, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array(dens_i, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array(temp_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
@@ -76,22 +76,19 @@ SUBROUTINE memory
CALL allocate_array( xphijm1, ips_i,ipe_i, ijs_i,ije_i)
! Curvature and geometry
- CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(Jacobian,izs,ize)
- CALL allocate_array(gxx, izs,ize)
- CALL allocate_array(gxy, izs,ize)
- CALL allocate_array(gyy, izs,ize)
- CALL allocate_array(gyz, izs,ize)
- CALL allocate_array(gxz, izs,ize)
- CALL allocate_array(gradzB,izs,ize)
- CALL allocate_array(gradxB,izs,ize)
- CALL allocate_array(hatB,izs,ize)
- CALL allocate_array(hatR,izs,ize)
- CALL allocate_array(Gamma1, izs,ize)
- call allocate_array(Gamma2, izs,ize)
- call allocate_array(Gamma3, izs, ize)
- call allocate_array(gradz_coeff, izs, ize)
+ CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
+ CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
+ CALL allocate_array( Jacobian, izs,ize, 0,1)
+ CALL allocate_array( gxx, izs,ize, 0,1)
+ CALL allocate_array( gxy, izs,ize, 0,1)
+ CALL allocate_array( gyy, izs,ize, 0,1)
+ CALL allocate_array( gyz, izs,ize, 0,1)
+ CALL allocate_array( gxz, izs,ize, 0,1)
+ CALL allocate_array( gradzB, izs,ize, 0,1)
+ CALL allocate_array( gradxB, izs,ize, 0,1)
+ CALL allocate_array( hatB, izs,ize, 0,1)
+ CALL allocate_array( hatR, izs,ize, 0,1)
+ call allocate_array(gradz_coeff, izs,ize, 0,1)
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 34a424de5adbfc3776ac3d27693d5de0888c58c1..ba1fc93e68dd211da1af6981d4d31be12edd1f75 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -15,6 +15,7 @@ SUBROUTINE moments_eq_rhs_e
USE prec_const
USE collision
use geometry
+ USE calculus, ONLY : interp_z, grad_z
IMPLICIT NONE
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
@@ -24,36 +25,42 @@ SUBROUTINE moments_eq_rhs_e
COMPLEX(dp) :: UNepm1j, UNepm1jp1, UNepm1jm1 ! Terms from mirror force with adiab moments
COMPLEX(dp) :: TColl ! terms of the rhs
COMPLEX(dp) :: i_ky
- REAL(dp) :: delta_p0, delta_p1, delta_p2
INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
+ ! To store derivatives and odd-even z grid interpolations
+ COMPLEX(dp), DIMENSION(izs:ize) :: ddznepp1j, ddznepm1j, &
+ nepp1j, nepp1jm1, nepm1j, nepm1jm1, nepm1jp1
! Measuring execution time
CALL cpu_time(t0_rhs)
+ ! Kinetic loops
ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices
- p_int = parray_e(ip) ! Hermite polynom. degree
- delta_p0 = 0._dp; delta_p1 = 0._dp; delta_p2 = 0._dp
- IF(p_int .EQ. 0) delta_p0 = 1._dp
- IF(p_int .EQ. 1) delta_p1 = 1._dp
- IF(p_int .EQ. 2) delta_p2 = 1._dp
-
+ p_int = parray_e(ip) ! Hermite polynom. degree
+ eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid
jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices
j_int = jarray_e(ij)
- ! Loop on kspace
- zloope : DO iz = izs,ize
- ! Obtain the index with an array that accounts for boundary conditions
- ! e.g. : 4 stencil with periodic BC, izarray(Nz+2) = 2, izarray(-1) = Nz-1
- izp1 = izarray(iz+1); izp2 = izarray(iz+2);
- izm1 = izarray(iz-1); izm2 = izarray(iz-2);
- !
- kxloope : DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! radial wavevector
- kyloope : DO iky = ikys,ikye
- ky = kyarray(iky) ! toroidal wavevector
- i_ky = imagu * ky ! toroidal derivative
- IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
- ! kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
- kperp2= kparray(ikx,iky,iz)**2
+ ! Spatial loops
+ kxloope : DO ikx = ikxs,ikxe
+ kx = kxarray(ikx) ! radial wavevector
+ kyloope : DO iky = ikys,ikye
+ ky = kyarray(iky) ! toroidal wavevector
+ i_ky = imagu * ky ! toroidal derivative
+ IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
+ ! Compute z derivatives and odd-even z interpolations
+ CALL grad_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,:),ddznepp1j(:))
+ CALL grad_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,:),ddznepm1j(:))
+ CALL interp_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,:), nepp1j (:))
+ CALL interp_z(eo,nadiab_moments_e(ip+1,ij-1,ikx,iky,:), nepp1jm1(:))
+ CALL interp_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,:), nepm1j (:))
+ CALL interp_z(eo,nadiab_moments_e(ip-1,ij+1,ikx,iky,:), nepm1jp1(:))
+ CALL interp_z(eo,nadiab_moments_e(ip-1,ij-1,ikx,iky,:), nepm1jm1(:))
+ zloope : DO iz = izs,ize
+ ! Obtain the index with an array that accounts for boundary conditions
+ ! e.g. : 4 stencil with periodic BC, izarray(Nz+2) = 2, izarray(-1) = Nz-1
+ izp1 = izarray(iz+1); izp2 = izarray(iz+2);
+ izm1 = izarray(iz-1); izm2 = izarray(iz-2);
+ !
+ kperp2= kparray(ikx,iky,iz,eo)**2
!! Compute moments mixing terms
! Perpendicular dynamic
@@ -71,33 +78,24 @@ SUBROUTINE moments_eq_rhs_e
Tpare = 0._dp; Tmir = 0._dp
IF(Nz .GT. 1) THEN
! ddz derivative for Landau damping term
- Tpare = xnepp1j(ip) * &
- ( onetwelfth*nadiab_moments_e(ip+1,ij,ikx,iky,izm2)&
- - twothird*nadiab_moments_e(ip+1,ij,ikx,iky,izm1)&
- + twothird*nadiab_moments_e(ip+1,ij,ikx,iky,izp1)&
- -onetwelfth*nadiab_moments_e(ip+1,ij,ikx,iky,izp2))&
- +xnepm1j(ip) * &
- ( onetwelfth*nadiab_moments_e(ip-1,ij,ikx,iky,izm2)&
- - twothird*nadiab_moments_e(ip-1,ij,ikx,iky,izm1)&
- + twothird*nadiab_moments_e(ip-1,ij,ikx,iky,izp1)&
- -onetwelfth*nadiab_moments_e(ip-1,ij,ikx,iky,izp2))
+ Tpare = xnepp1j(ip) * ddznepp1j(iz) + xnepm1j(ip) * ddznepm1j(iz)
! Mirror terms
- Tnepp1j = ynepp1j(ip,ij) * nadiab_moments_e(ip+1,ij ,ikx,iky,iz)
- Tnepp1jm1 = ynepp1jm1(ip,ij) * nadiab_moments_e(ip+1,ij-1,ikx,iky,iz)
- Tnepm1j = ynepm1j(ip,ij) * nadiab_moments_e(ip-1,ij ,ikx,iky,iz)
- Tnepm1jm1 = ynepm1jm1(ip,ij) * nadiab_moments_e(ip-1,ij-1,ikx,iky,iz)
+ Tnepp1j = ynepp1j (ip,ij) * nepp1j (iz)
+ Tnepp1jm1 = ynepp1jm1(ip,ij) * nepp1jm1(iz)
+ Tnepm1j = ynepm1j (ip,ij) * nepm1j (iz)
+ Tnepm1jm1 = ynepm1jm1(ip,ij) * nepm1jm1(iz)
! Trapping terms
- UNepm1j = zNepm1j(ip,ij) * nadiab_moments_e(ip-1,ij ,ikx,iky,iz)
- UNepm1jp1 = zNepm1jp1(ip,ij) * nadiab_moments_e(ip-1,ij+1,ikx,iky,iz)
- UNepm1jm1 = zNepm1jm1(ip,ij) * nadiab_moments_e(ip-1,ij-1,ikx,iky,iz)
+ UNepm1j = znepm1j (ip,ij) * nepm1j (iz)
+ UNepm1jp1 = znepm1jp1(ip,ij) * nepm1jp1(iz)
+ UNepm1jm1 = znepm1jm1(ip,ij) * nepm1jm1(iz)
Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + UNepm1j + UNepm1jp1 + UNepm1jm1
ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_e(ij,ikx,iky,iz) &
- + xphijp1(ip,ij)*kernel_e(ij+1,ikx,iky,iz) &
- + xphijm1(ip,ij)*kernel_e(ij-1,ikx,iky,iz) )
+ Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_e(ij,ikx,iky,iz,eo) &
+ + xphijp1(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) &
+ + xphijm1(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo) )
ELSE
Tphi = 0._dp
ENDIF
@@ -114,11 +112,11 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz)*hatR(iz)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
+ - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
! Parallel coupling (Landau Damping)
- - Tpare*inv_deltaz*gradz_coeff(iz) &
+ - Tpare*inv_deltaz*gradz_coeff(iz,eo) &
! Mirror term (parallel magnetic gradient)
- - gradzB(iz)* Tmir *gradz_coeff(iz) &
+ - gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) &
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Electrostatic background gradients
@@ -133,10 +131,9 @@ SUBROUTINE moments_eq_rhs_e
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) - Sepj(ip,ij,ikx,iky,iz)
ENDIF
-
- END DO kyloope
- END DO kxloope
- END DO zloope
+ END DO zloope
+ END DO kyloope
+ END DO kxloope
END DO jloope
END DO ploope
@@ -161,6 +158,7 @@ SUBROUTINE moments_eq_rhs_i
USE model
USE prec_const
USE collision
+ USE calculus, ONLY : interp_z, grad_z
IMPLICIT NONE
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
@@ -170,81 +168,73 @@ SUBROUTINE moments_eq_rhs_i
COMPLEX(dp) :: UNipm1j, UNipm1jp1, UNipm1jm1 ! Terms from mirror force with adiab moments
COMPLEX(dp) :: TColl ! terms of the rhs
COMPLEX(dp) :: i_ky
- REAL(dp) :: delta_p0, delta_p1, delta_p2
INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
-
+ ! To store derivatives and odd-even z grid interpolations
+ COMPLEX(dp), DIMENSION(izs:ize) :: ddznipp1j, ddznipm1j, &
+ nipp1j, nipp1jm1, nipm1j, nipm1jm1, nipm1jp1
! Measuring execution time
CALL cpu_time(t0_rhs)
+ ! Kinetic loops
ploopi : DO ip = ips_i, ipe_i ! Hermite loop
- p_int= parray_i(ip) ! Hermite degree
-
- delta_p0 = 0._dp; delta_p1 = 0._dp; delta_p2 = 0._dp
- IF(p_int .EQ. 0) delta_p0 = 1._dp
- IF(p_int .EQ. 1) delta_p1 = 1._dp
- IF(p_int .EQ. 2) delta_p2 = 1._dp
-
+ p_int = parray_i(ip) ! Hermite degree
+ eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
j_int = jarray_i(ij)
- ! Loop on kspace
- zloopi : DO iz = izs,ize
- ! Obtain the index with an array that accounts for boundary conditions
- ! e.g. : 4 stencil with periodic BC, izarray(Nz+2) = 2, izarray(-1) = Nz-1
- izp1 = izarray(iz+1); izp2 = izarray(iz+2);
- izm1 = izarray(iz-1); izm2 = izarray(iz-2);
- !
- kxloopi : DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! radial wavevector
- kyloopi : DO iky = ikys,ikye
- ky = kyarray(iky) ! toroidal wavevector
- i_ky = imagu * ky ! toroidal derivative
- IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
+ ! Spatial loops
+ kxloopi : DO ikx = ikxs,ikxe
+ kx = kxarray(ikx) ! radial wavevector
+ kyloopi : DO iky = ikys,ikye
+ ky = kyarray(iky) ! toroidal wavevector
+ i_ky = imagu * ky ! toroidal derivative
+ IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
+ ! Compute z derivatives and odd-even z interpolations
+ CALL grad_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,:),ddznipp1j(:))
+ CALL grad_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,:),ddznipm1j(:))
+ CALL interp_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,:), nipp1j (:))
+ CALL interp_z(eo,nadiab_moments_i(ip+1,ij-1,ikx,iky,:), nipp1jm1(:))
+ CALL interp_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,:), nipm1j (:))
+ CALL interp_z(eo,nadiab_moments_i(ip-1,ij+1,ikx,iky,:), nipm1jp1(:))
+ CALL interp_z(eo,nadiab_moments_i(ip-1,ij-1,ikx,iky,:), nipm1jm1(:))
+ zloopi : DO iz = izs,ize
! kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
- kperp2= kparray(ikx,iky,iz)**2
+ kperp2= kparray(ikx,iky,iz,eo)**2
!! Compute moments mixing terms
! Perpendicular dynamic
! term propto n_i^{p,j}
- Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip,ij,ikx,iky,iz)
+ Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,ikx,iky,iz)
! term propto n_i^{p+2,j}
- Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij,ikx,iky,iz)
+ Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,ikx,iky,iz)
! term propto n_i^{p-2,j}
- Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij,ikx,iky,iz)
+ Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,ikx,iky,iz)
! term propto n_e^{p,j+1}
- Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip,ij+1,ikx,iky,iz)
+ Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,ikx,iky,iz)
! term propto n_e^{p,j-1}
- Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip,ij-1,ikx,iky,iz)
+ Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,ikx,iky,iz)
! Parallel dynamic
Tpari = 0._dp; Tmir = 0._dp
IF(Nz .GT. 1) THEN
! term propto N_i^{p,j+1}, centered FDF
- Tpari = xnipp1j(ip) * &
- ( onetwelfth*nadiab_moments_i(ip+1,ij,ikx,iky,izm2)&
- - twothird*nadiab_moments_i(ip+1,ij,ikx,iky,izm1)&
- + twothird*nadiab_moments_i(ip+1,ij,ikx,iky,izp1)&
- -onetwelfth*nadiab_moments_i(ip+1,ij,ikx,iky,izp2))&
- +xnipm1j(ip) * &
- ( onetwelfth*nadiab_moments_i(ip-1,ij,ikx,iky,izm2)&
- - twothird*nadiab_moments_i(ip-1,ij,ikx,iky,izm1)&
- + twothird*nadiab_moments_i(ip-1,ij,ikx,iky,izp1)&
- -onetwelfth*nadiab_moments_i(ip-1,ij,ikx,iky,izp2))
+ Tpari = xnipp1j(ip) * ddznipp1j(iz) + xnipm1j(ip) * ddznipm1j(iz)
+
! Mirror terms
- Tnipp1j = ynipp1j(ip,ij) * nadiab_moments_i(ip+1,ij ,ikx,iky,iz)
- Tnipp1jm1 = ynipp1jm1(ip,ij) * nadiab_moments_i(ip+1,ij-1,ikx,iky,iz)
- Tnipm1j = ynipm1j(ip,ij) * nadiab_moments_i(ip-1,ij ,ikx,iky,iz)
- Tnipm1jm1 = ynipm1jm1(ip,ij) * nadiab_moments_i(ip-1,ij-1,ikx,iky,iz)
+ Tnipp1j = ynipp1j (ip,ij) * nipp1j (iz)
+ Tnipp1jm1 = ynipp1jm1(ip,ij) * nipp1jm1(iz)
+ Tnipm1j = ynipm1j (ip,ij) * nipm1j (iz)
+ Tnipm1jm1 = ynipm1jm1(ip,ij) * nipm1jm1(iz)
! Trapping terms
- Unipm1j = znipm1j(ip,ij) * nadiab_moments_i(ip-1,ij ,ikx,iky,iz)
- Unipm1jp1 = znipm1jp1(ip,ij) * nadiab_moments_i(ip-1,ij+1,ikx,iky,iz)
- Unipm1jm1 = znipm1jm1(ip,ij) * nadiab_moments_i(ip-1,ij-1,ikx,iky,iz)
+ UNipm1j = znipm1j (ip,ij) * nipm1j (iz)
+ UNipm1jp1 = znipm1jp1(ip,ij) * nipm1jp1(iz)
+ UNipm1jm1 = znipm1jm1(ip,ij) * nipm1jm1(iz)
Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + UNipm1j + UNipm1jp1 + UNipm1jm1
ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_i(ij,ikx,iky,iz) &
- + xphijp1(ip,ij)*kernel_i(ij+1,ikx,iky,iz) &
- + xphijm1(ip,ij)*kernel_i(ij-1,ikx,iky,iz) )
+ Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_i(ij,ikx,iky,iz,eo) &
+ + xphijp1(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) &
+ + xphijm1(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo) )
ELSE
Tphi = 0._dp
ENDIF
@@ -261,11 +251,11 @@ SUBROUTINE moments_eq_rhs_i
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz)*hatR(iz)*(Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
+ - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo)*(Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
! Parallel coupling (Landau Damping)
- - Tpari*inv_deltaz*gradz_coeff(iz) &
+ - Tpari*inv_deltaz*gradz_coeff(iz,eo) &
! Mirror term (parallel magnetic gradient)
- - gradzB(iz)*Tmir*gradz_coeff(iz) &
+ - gradzB(iz,eo)*Tmir*gradz_coeff(iz,eo) &
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Electrostatic background gradients
@@ -280,10 +270,9 @@ SUBROUTINE moments_eq_rhs_i
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) - Sipj(ip,ij,ikx,iky,iz)
ENDIF
-
- END DO kyloopi
- END DO kxloopi
- END DO zloopi
+ END DO zloopi
+ END DO kyloopi
+ END DO kxloopi
END DO jloopi
END DO ploopi
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 06235c2fae37f67ebcef16cc6af8f5708e6655f0..da59158ef16b5bd1d948b78ba0c7676fa11a1af4 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -57,6 +57,7 @@ SUBROUTINE evaluate_kernels
INTEGER :: j_int
REAL(dp) :: j_dp, y_, kp2_, kx_, ky_
+DO eo = 0,1
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
DO iz = izs,ize
@@ -65,20 +66,21 @@ DO iz = izs,ize
DO ij = ijsg_e, ijeg_e
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp)
- y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz)**2
- kernel_e(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
+ y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz,eo)**2
+ kernel_e(ij,ikx,iky,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
ENDIF
!!!!! Ion kernels !!!!!
DO ij = ijsg_i, ijeg_i
j_int = jarray_i(ij)
j_dp = REAL(j_int,dp)
- y_ = sigmai2_taui_o2 * kparray(ikx,iky,iz)**2
- kernel_i(ij,ikx,iky,iz) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
+ y_ = sigmai2_taui_o2 * kparray(ikx,iky,iz,eo)**2
+ kernel_i(ij,ikx,iky,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
ENDDO
ENDDO
ENDDO
+ENDDO
END SUBROUTINE evaluate_kernels
!******************************************************************************!
@@ -98,6 +100,7 @@ SUBROUTINE evaluate_poisson_op
kxloop: DO ikx = ikxs,ikxe
kyloop: DO iky = ikys,ikye
zloop: DO iz = izs,ize
+ ! This term is evalued on the even z grid since poisson uses only p=0 and phi
IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
inv_poisson_op(ikx, iky, iz) = 0._dp
ELSE
@@ -105,7 +108,7 @@ SUBROUTINE evaluate_poisson_op
! loop over n only if the max polynomial degree
pol_i = 0._dp
DO ini=1,jmaxi+1
- pol_i = pol_i + qi2_taui*kernel_i(ini,ikx,iky,iz)**2 ! ... sum recursively ...
+ pol_i = pol_i + qi2_taui*kernel_i(ini,ikx,iky,iz,0)**2 ! ... sum recursively ...
END DO
!!!!!!!!!!!!! Electron contribution\
pol_e = 0._dp
@@ -113,7 +116,7 @@ SUBROUTINE evaluate_poisson_op
IF (KIN_E) THEN
! loop over n only if the max polynomial degree
DO ine=1,jmaxe+1 ! ine = n+1
- pol_e = pol_e + qe2_taue*kernel_e(ine,ikx,iky,iz)**2 ! ... sum recursively ...
+ pol_e = pol_e + qe2_taue*kernel_e(ine,ikx,iky,iz,0)**2 ! ... sum recursively ...
END DO
! Adiabatic model
ELSE
diff --git a/src/poisson.F90 b/src/poisson.F90
index 8cc242ca37e3e6b4d6f0bc2bf650436554b0bb81..95b5342a1b1c2fedb165055973321bd9e0dd08b4 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -6,7 +6,8 @@ SUBROUTINE poisson
USE array
USE fields
USE grid
- USE utility
+ USE calculus, ONLY : simpson_rule_z
+ USE utility, ONLY : manual_3D_bcast
use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD, KIN_E
USE processing, ONLY : compute_density
USE prec_const
@@ -31,7 +32,7 @@ SUBROUTINE poisson
rho_i = 0._dp
DO ini=1,jmaxi+1
rho_i = rho_i &
- +q_i*kernel_i(ini,ikx,iky,:)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)
+ +q_i*kernel_i(ini,ikx,iky,:,0)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)
END DO
!!!! Compute electron gyro charge density
@@ -39,15 +40,15 @@ SUBROUTINE poisson
IF (KIN_E) THEN ! Kinetic electrons
DO ine=1,jmaxe+1
rho_e = rho_e &
- +q_e*kernel_e(ine,ikx,iky,:)*moments_e(ip0_e,ine, ikx,iky,:,updatetlevel)
+ +q_e*kernel_e(ine,ikx,iky,:,0)*moments_e(ip0_e,ine, ikx,iky,:,updatetlevel)
END DO
ELSE ! Adiabatic electrons
! Adiabatic charge density (linked to flux averaged phi)
fa_phi = 0._dp
IF(kyarray(iky).EQ.0._dp) THEN
DO ini=1,jmaxi+1
- rho_e(:) = Jacobian(:)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)&
- *kernel_i(ini,ikx,iky,:)*(inv_poisson_op(ikx,iky,:)-1._dp)
+ rho_e(:) = Jacobian(:,0)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)&
+ *kernel_i(ini,ikx,iky,:,0)*(inv_poisson_op(ikx,iky,:)-1._dp)
call simpson_rule_z(rho_e,intf_)
fa_phi = fa_phi + intf_
ENDDO
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 2866a44d0fad2ba3d056b1c836f2e59fa223a0b6..5419e3b7cc1ae80527f919ec4d71a3175874ce63 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -35,7 +35,7 @@ SUBROUTINE compute_radial_ion_transport
imagu * ky_ * moments_i(ip0_i,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,iz) * moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
+ imagu * ky_ * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
ENDDO
ENDDO
ENDDO
@@ -87,20 +87,20 @@ SUBROUTINE compute_radial_heatflux
DO ij = ijs_i, ije_i
j_dp = REAL(ij-1,dp)
hflux_local = hflux_local - imagu*ky_*CONJG(phi(ikx,iky,iz))&
- *(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(ip0_i,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(ip2_i,ij,ikx,iky,iz,updatetlevel))
+ *(2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz,0) - j_dp*kernel_i(ij-1,ikx,iky,iz,0))&
+ * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+ + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel))
ENDDO
IF(KIN_E) THEN
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
hflux_local = hflux_local - imagu*ky_*CONJG(phi(ikx,iky,iz))&
- *(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))&
+ *(2._dp/3._dp *(2._dp*j_dp * kernel_e(ij ,ikx,iky,iz,0) &
+ -(j_dp+1) * kernel_e(ij+1,ikx,iky,iz,0) &
+ -j_dp * kernel_e(ij-1,ikx,iky,iz,0))&
*(moments_e(ip0_e,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(ip2_e,ij,ikx,iky,iz,updatetlevel))
+ +q_e/tau_e * kernel_e(ij ,ikx,iky,iz,0) * phi(ikx,iky,iz)) &
+ +SQRT2/3._dp * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel))
ENDDO
ENDIF
ENDDO
@@ -146,7 +146,7 @@ SUBROUTINE compute_nadiab_moments
DO ij=ijsg_e,ijeg_e
nadiab_moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
= moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
- + qe_taue*kernel_e(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
+ + qe_taue*kernel_e(ij,ikxs:ikxe,ikys:ikye,izs:ize,0)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
ENDDO
ELSE
nadiab_moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize) &
@@ -160,7 +160,7 @@ SUBROUTINE compute_nadiab_moments
DO ij=ijsg_i,ijeg_i
nadiab_moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
= moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
- + qi_taui*kernel_i(ij,ikxs:ikxe,ikys:ikye,izs:ize)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
+ + qi_taui*kernel_i(ij,ikxs:ikxe,ikys:ikye,izs:ize,0)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
ENDDO
ELSE
nadiab_moments_i(ip,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize) &
@@ -187,15 +187,15 @@ 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,iz) * &
- (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz))
+ dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky,iz,0) * &
+ (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz))
ENDDO
ENDIF
! 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,iz) * &
- (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iz))
+ dens_i(ikx,iky,iz) = dens_i(ikx,iky,iz) + kernel_i(ij,ikx,iky,iz,0) * &
+ (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz))
ENDDO
ENDDO
ENDDO
@@ -227,9 +227,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,iz) - (j_dp+1)*kernel_e(ij+1,ikx,iky,iz) - j_dp*kernel_e(ij-1,ikx,iky,iz))&
- * (moments_e(ip0_e,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(ip2_e,ij,ikx,iky,iz,updatetlevel)
+ 2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_e(ij+1,ikx,iky,iz,0) - j_dp*kernel_e(ij-1,ikx,iky,iz,0))&
+ * (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+ + SQRT2/3._dp * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel)
ENDDO
ENDIF
! ion temperature
@@ -237,9 +237,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,iz) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz) - j_dp*kernel_i(ij-1,ikx,iky,iz))&
- * (moments_i(ip0_i,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(ip2_i,ij,ikx,iky,iz,updatetlevel)
+ 2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz,0) - j_dp*kernel_i(ij-1,ikx,iky,iz,0))&
+ * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+ + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel)
ENDDO
ENDDO
ENDDO
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index def92f7a9f0585a36c2c428461e443089e5fde1d..0bb98587c82da2e93cb14d1c7e1ecdd797d277f5 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -5,7 +5,7 @@ USE grid
USE fields
USE diagnostics_par
USE time_integration
-
+USE model, ONLY: KIN_E
IMPLICIT NONE
INTEGER :: rank, sz_, n_
@@ -32,15 +32,17 @@ CONTAINS
! Open file
CALL openf(rstfile, fidrst,mpicomm=comm0)
! Get the checkpoint moments degrees to allocate memory
+ IF (KIN_E) THEN
CALL getatt(fidrst,"/data/input/" , "pmaxe", pmaxe_cp)
CALL getatt(fidrst,"/data/input/" , "jmaxe", jmaxe_cp)
+ ENDIF
CALL getatt(fidrst,"/data/input/" , "pmaxi", pmaxi_cp)
CALL getatt(fidrst,"/data/input/" , "jmaxi", jmaxi_cp)
- IF (my_id .EQ. 0) WRITE(*,*) "Pe_cp = ", pmaxe_cp
- IF (my_id .EQ. 0) WRITE(*,*) "Je_cp = ", jmaxe_cp
+ IF (my_id .EQ. 0) WRITE(*,*) "Pi_cp = ", pmaxi_cp
+ IF (my_id .EQ. 0) WRITE(*,*) "Ji_cp = ", jmaxi_cp
CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
- IF ((pmaxe_cp .NE. pmaxe) .OR. (jmaxe_cp .NE. jmaxe) .OR.&
+ IF ((KIN_E .AND. ((pmaxe_cp .NE. pmaxe) .OR. (jmaxe_cp .NE. jmaxe))) .OR.&
(pmaxi_cp .NE. pmaxi) .OR. (jmaxi_cp .NE. jmaxi)) THEN
IF(my_id.EQ.0)WRITE(*,*) '! Extending the polynomials basis !'
CALL load_output_adapt_pj
@@ -48,15 +50,15 @@ CONTAINS
! Find the last results of the checkpoint file by iteration
n_ = n0+1
- WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_ ! start with moments_e/000001
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! start with moments_e/000001
DO WHILE (isdataset(fidrst, dset_name)) ! If n_ is not a file we stop the loop
n_ = n_ + 1
- WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_ ! updtate file number
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! updtate file number
ENDDO
n_ = n_ - 1 ! n_ is not a file so take the previous one n_-1
! Read time dependent attributes to continue simulation
- WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
CALL getatt(fidrst, dset_name, 'cstep', cstep)
CALL getatt(fidrst, dset_name, 'time', time)
CALL getatt(fidrst, dset_name, 'jobnum', jobnum)
@@ -67,8 +69,10 @@ CONTAINS
IF(my_id.EQ.0) WRITE(*,*) '.. restart from t = ', time
! Read state of system from checkpoint file
+ IF (KIN_E) THEN
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3/))
+ ENDIF
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3/))
@@ -93,18 +97,16 @@ CONTAINS
IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Resume from ", rstfile
! Open file
CALL openf(rstfile, fidrst,mpicomm=comm0)
+ !!!!!!!!! Load electron moments
+ IF (KIN_E) THEN
! Get the checkpoint moments degrees to allocate memory
CALL getatt(fidrst,"/data/input/" , "pmaxe", pmaxe_cp)
CALL getatt(fidrst,"/data/input/" , "jmaxe", jmaxe_cp)
- CALL getatt(fidrst,"/data/input/" , "pmaxi", pmaxi_cp)
- CALL getatt(fidrst,"/data/input/" , "jmaxi", jmaxi_cp)
IF (my_id .EQ. 0) WRITE(*,*) "Pe_cp = ", pmaxe_cp
IF (my_id .EQ. 0) WRITE(*,*) "Je_cp = ", jmaxe_cp
CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
-
! Allocate the required size to load checkpoints moments
CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_ ! start with moments_e/000001
@@ -113,16 +115,12 @@ CONTAINS
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_ ! updtate file number
ENDDO
n_ = n_ - 1 ! n_ is not a file so take the previous one n_-1
-
! Read state of system from checkpoint file and load every moment to change the distribution
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
CALL getarrnd(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
- WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
- CALL getarrnd(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
-
! Initialize simulation moments array with checkpoints ones
! (they may have a larger number of polynomials, set to 0 at the begining)
- moments_e = 0._dp; moments_i = 0._dp
+ moments_e = 0._dp;
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
DO ikx=ikxs,ikxe
@@ -134,7 +132,33 @@ CONTAINS
ENDDO
ENDDO
ENDDO
+ ! Deallocate checkpoint arrays
+ DEALLOCATE(moments_e_cp)
+ ENDIF
+ !!!!!!! Load ion moments
+ ! Get the checkpoint moments degrees to allocate memory
+ CALL getatt(fidrst,"/data/input/" , "pmaxi", pmaxi_cp)
+ CALL getatt(fidrst,"/data/input/" , "jmaxi", jmaxi_cp)
+ IF (my_id .EQ. 0) WRITE(*,*) "Pi_cp = ", pmaxi_cp
+ IF (my_id .EQ. 0) WRITE(*,*) "Ji_cp = ", jmaxi_cp
+ CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
+ ! Allocate the required size to load checkpoints moments
+ CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
+ ! Find the last results of the checkpoint file by iteration
+ n_ = n0+1
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! start with moments_e/000001
+ DO WHILE (isdataset(fidrst, dset_name)) ! If n_ is not a file we stop the loop
+ n_ = n_ + 1
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! updtate file number
+ ENDDO
+ n_ = n_ - 1 ! n_ is not a file so take the previous one n_-1
+ ! Read state of system from checkpoint file and load every moment to change the distribution
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
+ CALL getarrnd(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
+ ! Initialize simulation moments array with checkpoints ones
+ ! (they may have a larger number of polynomials, set to 0 at the begining)
+ moments_i = 0._dp
DO ip=1,pmaxi_cp+1
DO ij=1,jmaxi_cp+1
DO ikx=ikxs,ikxe
@@ -147,7 +171,6 @@ CONTAINS
ENDDO
ENDDO
! Deallocate checkpoint arrays
- DEALLOCATE(moments_e_cp)
DEALLOCATE(moments_i_cp)
! Read time dependent attributes to continue simulation
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 85bcb0f4f0277904b7ad93de01a72c066912e435..e920419f6b3b8dfaacda6b07fb54bd74c4eacf9f 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='d29af00-dirty')
+parameter (VERSION='662e3f7-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Fri Oct 29 18:03:02 CEST 2021')
+parameter (EXECDATE='Mon Nov 1 11:15:33 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 85bcb0f4f0277904b7ad93de01a72c066912e435..e920419f6b3b8dfaacda6b07fb54bd74c4eacf9f 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='d29af00-dirty')
+parameter (VERSION='662e3f7-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Fri Oct 29 18:03:02 CEST 2021')
+parameter (EXECDATE='Mon Nov 1 11:15:33 CET 2021')
parameter (HOST ='spcpc606')
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index 556ed5c21220a822da087c70aa4bed14d051f5f9..a295575e5464fdf0efdf415e2f9f08f84fd7dc0c 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -4,8 +4,7 @@ MODULE utility
use prec_const
IMPLICIT NONE
- PUBLIC :: manual_2D_bcast, manual_3D_bcast,&
- simpson_rule_z, o2e_z, e2o_z
+ PUBLIC :: manual_2D_bcast, manual_3D_bcast
CONTAINS
@@ -150,74 +149,4 @@ SUBROUTINE manual_3D_bcast(field_)
ENDIF
END SUBROUTINE manual_3D_bcast
-SUBROUTINE simpson_rule_z(f,intf)
- ! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
- !from molix BJ Frei
- use prec_const
- use grid
- !
- implicit none
- !
- complex(dp),dimension(izs:ize), intent(in) :: f
- COMPLEX(dp), intent(out) :: intf
- !
- COMPLEX(dp) :: buff_
- !
- IF(Nz .GT. 1) THEN
- IF(mod(Nz,2) .ne. 0 ) THEN
- ERROR STOP 'Simpson rule: Nz must be an even number !!!!'
- ENDIF
- !
- buff_ = 0._dp
- !
- DO iz = izs, Nz/2
- IF(iz .eq. Nz/2) THEN ! ... iz = ize
- buff_ = buff_ + (f(izs) + 4._dp*f(ize) + f(ize-1 ))
- ELSE
- buff_ = buff_ + (f(2*iz+1) + 4._dp*f(2*iz) + f(2*iz-1 ))
- ENDIF
- ENDDO
- !
- !
- intf = buff_*deltaz/3._dp
- !
- ELSE
- intf = f(izs)
- ENDIF
-END SUBROUTINE simpson_rule_z
-
-SUBROUTINE o2e_z(fo,fe)
- ! gives the value of a field from the odd grid to the even one
- use prec_const
- use grid
- !
- implicit none
- !
- COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
- COMPLEX(dp),dimension(1:Nz), intent(out) :: fe !
- !
- DO iz = 2,Nz
- fe(iz) = 0.5_dp*(fo(iz)+fo(iz-1))
- ENDDO
- ! Periodic boundary conditions
- fe(1) = 0.5_dp*(fo(1) + fo(Nz))
-END SUBROUTINE o2e_z
-
-SUBROUTINE e2o_z(fe,fo)
- ! gives the value of a field from the even grid to the odd one
- use prec_const
- use grid
- !
- implicit none
- !
- COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
- COMPLEX(dp),dimension(1:Nz), intent(out) :: fo
- !
- DO iz = 1,Nz-1
- fo(iz) = 0.5_dp*(fe(iz+1)+fe(iz))
- ENDDO
- ! Periodic boundary conditions
- fo(Nz) = 0.5_dp*(fe(1) + fe(Nz))
-END SUBROUTINE e2o_z
-
END MODULE utility