From e581e2d06aa83c1d67b9845f409a799adfc77049 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 26 Jan 2023 13:39:51 +0100
Subject: [PATCH] Major change of the module. The e-i routines have been merged
in one with input parameters.
---
src/moments_eq_rhs_mod.F90 | 430 ++++++++++++++-----------------------
1 file changed, 167 insertions(+), 263 deletions(-)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index eaefef5d..4ab8c91d 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -4,313 +4,217 @@ MODULE moments_eq_rhs
CONTAINS
SUBROUTINE compute_moments_eq_rhs
- USE model, only: KIN_E
- IMPLICIT NONE
-
- IF(KIN_E) CALL moments_eq_rhs_e
- CALL moments_eq_rhs_i
- !! BETA TESTING !!
- IF(.false.) CALL add_Maxwellian_background_terms
-END SUBROUTINE compute_moments_eq_rhs
-!_____________________________________________________________________________!
-!_____________________________________________________________________________!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!! Electrons moments RHS !!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!_____________________________________________________________________________!
-SUBROUTINE moments_eq_rhs_e
-
- USE basic
- USE time_integration
- USE array, ONLY: xnepj, xnepp2j, xnepm2j, xnepjp1, xnepjm1, xnepp1j, xnepm1j,&
- ynepp1j, ynepp1jm1, ynepm1j, ynepm1jm1, &
- znepm1j, znepm1jp1, znepm1jm1, &
- xphij_e, xphijp1_e, xphijm1_e, xpsij_e, xpsijp1_e, xpsijm1_e,&
- nadiab_moments_e, ddz_nepj, interp_nepj, moments_rhs_e, Sepj,&
- TColl_e, ddzND_nepj, kernel_e
- USE fields, ONLY: phi, psi, moments_e
- USE grid
USE model
- USE prec_const
- USE collision
- USE geometry, ONLY: gradz_coeff, dBdz, Ckxky, hatB_NL, hatB
- USE calculus, ONLY : interp_z, grad_z, grad_z2
- IMPLICIT NONE
-
- INTEGER :: p_int, j_int ! loops indices and polynom. degrees
- REAL(dp) :: kx, ky, kperp2
- COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives
- COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
- COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi
- COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: i_kx,i_ky,phikykxz, psikykxz
-
- ! Measuring execution time
- CALL cpu_time(t0_rhs)
-
- ! Spatial loops
- zloope : DO iz = izs,ize
- kxloope : DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! radial wavevector
- i_kx = imagu * kx ! toroidal derivative
-
- kyloope : DO iky = ikys,ikye
- ky = kyarray(iky) ! toroidal wavevector
- i_ky = imagu * ky ! toroidal derivative
- phikykxz = phi(iky,ikx,iz)! tmp phi value
- psikykxz = psi(iky,ikx,iz)! tmp psi value
-
- ! Kinetic loops
- jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1
- j_int = jarray_e(ij)
-
- ploope : DO ip = ips_e, ipe_e ! Hermite loop
- p_int = parray_e(ip) ! Hermite degree
- eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
- kperp2= kparray(iky,ikx,iz,eo)**2
-
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN
- !! Compute moments mixing terms
- Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
- ! Perpendicular dynamic
- ! term propto n_e^{p,j}
- Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz)
- ! term propto n_e^{p+2,j}
- Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz)
- ! term propto n_e^{p-2,j}
- Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz)
- ! term propto n_e^{p,j+1}
- Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz)
- ! term propto n_e^{p,j-1}
- Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz)
- ! Tperp
- Tperp = Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1
- ! Parallel dynamic
- ! ddz derivative for Landau damping term
- Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) &
- + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz)
- ! Mirror terms
- Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz)
- Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz)
- Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz)
- Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz)
- ! Trapping terms
- Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz)
- Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz)
- Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz)
-
- Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1
- !! Electrical potential term
- IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = (xphij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
- + xphijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) &
- + xphijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz
- ELSE
- Tphi = 0._dp
- ENDIF
-
- !! Vector potential term
- IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3
- Tpsi = (xpsij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
- + xpsijp1_e(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) &
- + xpsijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*psikykxz
- ELSE
- Tpsi = 0._dp
- ENDIF
-
- !! Sum of all RHS terms
- moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
- ! Perpendicular magnetic gradient/curvature effects
- -imagu*Ckxky(iky,ikx,iz,eo) * Tperp&
- ! Perpendicular pressure effects
- -i_ky*beta*dpdx * Tperp&
- ! Parallel coupling (Landau Damping)
- -gradz_coeff(iz,eo) * Tpar &
- ! Mirror term (parallel magnetic gradient)
- -dBdz(iz,eo)*gradz_coeff(iz,eo) * Tmir&
- ! Drives (density + temperature gradients)
- -i_ky * (Tphi - Tpsi) &
- ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- -mu_x*diff_kx_coeff*kx**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
- -mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
- ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
- -mu_z*diff_dz_coeff*ddzND_nepj(ip,ij,iky,ikx,iz) &
- ! Collision term
- +TColl_e(ip,ij,iky,ikx,iz) &
- ! Nonlinear term
- -hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz)
-
- ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
- ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! ! (not used often so not parallelized)
- ! moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
- ! moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) &
- ! + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel)
-
- ELSE
- moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
- ENDIF
- END DO ploope
- END DO jloope
- END DO kyloope
- END DO kxloope
- END DO zloope
-
- ! Execution time end
- CALL cpu_time(t1_rhs)
- tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
-
-END SUBROUTINE moments_eq_rhs_e
-!_____________________________________________________________________________!
-!_____________________________________________________________________________!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!! Ions moments RHS !!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!_____________________________________________________________________________!
-SUBROUTINE moments_eq_rhs_i
-
- USE basic
- USE time_integration, ONLY: updatetlevel
- USE array, ONLY: xnipj, xnipp2j, xnipm2j, xnipjp1, xnipjm1, xnipp1j, xnipm1j,&
- ynipp1j, ynipp1jm1, ynipm1j, ynipm1jm1, &
- znipm1j, znipm1jp1, znipm1jm1, &
- xphij_i, xphijp1_i, xphijm1_i, xpsij_i, xpsijp1_i, xpsijm1_i,&
- nadiab_moments_i, ddz_nipj, interp_nipj, moments_rhs_i, Sipj,&
- TColl_i, ddzND_nipj, kernel_i
- USE fields, ONLY: phi, psi, moments_i
+ USE array
+ USE fields
USE grid
- USE model
+ USE basic
USE prec_const
USE collision
+ USE time_integration
USE geometry, ONLY: gradz_coeff, dBdz, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
- INTEGER :: p_int, j_int ! loops indices and polynom. degrees
- REAL(dp) :: kx, ky, kperp2
- COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1
- COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
- COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi
- COMPLEX(dp) :: i_kx, i_ky, phikykxz, psikykxz
-
- ! Measuring execution time
- CALL cpu_time(t0_rhs)
-
-
- ! Spatial loops
- zloopi : DO iz = izs,ize
- kxloopi : DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! radial wavevector
- i_kx = imagu * kx ! toroidal derivative
- kyloopi : DO iky = ikys,ikye
- ky = kyarray(iky) ! toroidal wavevector
- i_ky = imagu * ky ! toroidal derivative
+ !compute ion moments_eq_rhs
+ CALL moments_eq_rhs(ips_i,ipe_i,ipgs_i,ipge_i,ijs_i,ije_i,ijgs_i,ijge_i,jarray_i,parray_i,&
+ xnipj, xnipp2j, xnipm2j, xnipjp1, xnipjm1, xnipp1j, xnipm1j,&
+ ynipp1j, ynipp1jm1, ynipm1j, ynipm1jm1, &
+ znipm1j, znipm1jp1, znipm1jm1, &
+ xphij_i, xphijp1_i, xphijm1_i, xpsij_i, xpsijp1_i, xpsijm1_i,&
+ kernel_i, nadiab_moments_i, ddz_nipj, interp_nipj, Sipj,&
+ moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel),&
+ TColl_i, ddzND_nipj, &
+ moments_rhs_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel))
+
+ !compute ion moments_eq_rhs
+ IF(KIN_E) &
+ CALL moments_eq_rhs(ips_e,ipe_e,ipgs_e,ipge_e,ijs_e,ije_e,ijgs_e,ijge_e,jarray_e,parray_e,&
+ xnepj, xnepp2j, xnepm2j, xnepjp1, xnepjm1, xnepp1j, xnepm1j,&
+ ynepp1j, ynepp1jm1, ynepm1j, ynepm1jm1, &
+ znepm1j, znepm1jp1, znepm1jm1, &
+ xphij_e, xphijp1_e, xphijm1_e, xpsij_e, xpsijp1_e, xpsijm1_e,&
+ kernel_e, nadiab_moments_e, ddz_nepj, interp_nepj, Sepj,&
+ moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel),&
+ TColl_e, ddzND_nepj,&
+ moments_rhs_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel))
+
+ CONTAINS
+ !_____________________________________________________________________________!
+ !_____________________________________________________________________________!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!! moments_ RHS computation !!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This routine assemble the RHS of the moment hierarchy equations. It uses
+ ! linear coefficients that are stored in arrays (xn*, yn* and zn*) computed in
+ ! numerics_mod.F90. Otherwise it simply adds the collision term TColl that is
+ ! computed in collision_mod.F90 and the nonlinear term Sapj computed in
+ ! nonlinear_mod.F90.
+ ! All arguments of the subroutines are inputs only except the last one,
+ ! moments_rhs_ that will contain the sum of every terms in the RHS.
+ !_____________________________________________________________________________!
+ SUBROUTINE moments_eq_rhs(ips,ipe,ipgs,ipge,ijs,ije,ijgs,ijge,jarray,parray,&
+ xnapj, xnapp2j, xnapm2j, xnapjp1, xnapjm1, xnapp1j, xnapm1j,&
+ ynapp1j, ynapp1jm1, ynapm1j, ynapm1jm1, &
+ znapm1j, znapm1jp1, znapm1jm1, &
+ xphij, xphijp1, xphijm1, xpsij, xpsijp1, xpsijm1,&
+ kernel, nadiab_moments, ddz_napj, interp_napj, Sapj,&
+ moments_, TColl_, ddzND_napj, moments_rhs_)
+
+ IMPLICIT NONE
+ !! INPUTS
+ INTEGER, INTENT(IN) :: ips, ipe, ipgs, ipge, ijs, ije, ijgs, ijge
+ INTEGER, DIMENSION(ips:ipe), INTENT(IN) :: parray
+ INTEGER, DIMENSION(ijs:ije), INTENT(IN) :: jarray
+ REAL(dp), DIMENSION(ips:ipe,ijs:ije), INTENT(IN) :: xnapj,&
+ ynapp1j, ynapm1j, ynapp1jm1, ynapm1jm1,&
+ znapm1j, znapm1jp1, znapm1jm1,&
+ xphij, xphijp1, xphijm1,&
+ xpsij, xpsijp1, xpsijm1
+ REAL(dp), DIMENSION(ips:ipe), INTENT(IN) :: &
+ xnapp1j, xnapm1j, xnapp2j, xnapm2j, xnapjp1, xnapjm1
+
+ REAL(dp), DIMENSION(ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge,0:1),INTENT(IN) :: kernel
+
+ COMPLEX(dp), DIMENSION(ipgs:ipge,ijgs:ijge,ikys:ikye,ikxs:ikxe,izgs:izge),INTENT(IN) ::&
+ moments_,nadiab_moments, ddz_napj, interp_napj, ddzND_napj
+ COMPLEX(dp), DIMENSION(ips:ipe,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize),INTENT(IN) :: Sapj, TColl_
+
+ !! OUTPUT
+ COMPLEX(dp), DIMENSION(ips:ipe,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize),INTENT(OUT) :: moments_rhs_
+
+ INTEGER :: p_int, j_int ! loops indices and polynom. degrees
+ REAL(dp) :: kx, ky, kperp2
+ COMPLEX(dp) :: Tnapj, Tnapp2j, Tnapm2j, Tnapjp1, Tnapjm1 ! Terms from b x gradB and drives
+ COMPLEX(dp) :: Tnapp1j, Tnapm1j, Tnapp1jm1, Tnapm1jm1 ! Terms from mirror force with non adiab moments_
+ COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi, Tpsi
+ COMPLEX(dp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
+ COMPLEX(dp) :: i_kx,i_ky,phikykxz, psikykxz
+
+ ! Measuring execution time
+ CALL cpu_time(t0_rhs)
+
+ ! Spatial loops
+ zloop : DO iz = izs,ize
+ kxloop : DO ikx = ikxs,ikxe
+ kx = kxarray(ikx) ! radial wavevector
+ i_kx = imagu * kx ! radial derivative
+
+ kyloop : DO iky = ikys,ikye
+ ky = kyarray(iky) ! binormal wavevector
+ i_ky = imagu * ky ! binormal derivative
phikykxz = phi(iky,ikx,iz)! tmp phi value
- psikykxz = psi(iky,ikx,iz)! tmp phi value
+ psikykxz = psi(iky,ikx,iz)! tmp psi value
- ! Kinetic loops
- jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
- j_int = jarray_i(ij)
+ ! Kinetic loops
+ jloop : DO ij = ijs, ije ! This loop is from 1 to jmaxi+1
+ j_int = jarray(ij)
- ploopi : DO ip = ips_i, ipe_i ! Hermite loop
- p_int = parray_i(ip) ! Hermite degree
- eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
- kperp2= kparray(iky,ikx,iz,eo)**2
+ ploop : DO ip = ips, ipe ! Hermite loop
+ p_int = parray(ip) ! Hermite degree
+ eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
+ kperp2= kparray(iky,ikx,iz,eo)**2
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN
- !! Compute moments mixing terms
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN
+ !! Compute moments_ mixing terms
Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
- ! term propto n_i^{p,j}
- Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz)
- ! term propto n_i^{p+2,j}
- Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz)
- ! term propto n_i^{p-2,j}
- Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz)
- ! term propto n_i^{p,j+1}
- Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz)
- ! term propto n_i^{p,j-1}
- Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz)
+ ! term propto n^{p,j}
+ Tnapj = xnapj(ip,ij)* nadiab_moments(ip,ij,iky,ikx,iz)
+ ! term propto n^{p+2,j}
+ Tnapp2j = xnapp2j(ip) * nadiab_moments(ip+pp2,ij,iky,ikx,iz)
+ ! term propto n^{p-2,j}
+ Tnapm2j = xnapm2j(ip) * nadiab_moments(ip-pp2,ij,iky,ikx,iz)
+ ! term propto n^{p,j+1}
+ Tnapjp1 = xnapjp1(ij) * nadiab_moments(ip,ij+1,iky,ikx,iz)
+ ! term propto n^{p,j-1}
+ Tnapjm1 = xnapjm1(ij) * nadiab_moments(ip,ij-1,iky,ikx,iz)
! Tperp
- Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1
+ Tperp = Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1
! Parallel dynamic
! ddz derivative for Landau damping term
- Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) &
- + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz)
+ Tpar = xnapp1j(ip) * ddz_napj(ip+1,ij,iky,ikx,iz) &
+ + xnapm1j(ip) * ddz_napj(ip-1,ij,iky,ikx,iz)
! Mirror terms
- Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz)
- Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz)
- Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz)
- Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz)
+ Tnapp1j = ynapp1j (ip,ij) * interp_napj(ip+1,ij ,iky,ikx,iz)
+ Tnapp1jm1 = ynapp1jm1(ip,ij) * interp_napj(ip+1,ij-1,iky,ikx,iz)
+ Tnapm1j = ynapm1j (ip,ij) * interp_napj(ip-1,ij ,iky,ikx,iz)
+ Tnapm1jm1 = ynapm1jm1(ip,ij) * interp_napj(ip-1,ij-1,iky,ikx,iz)
! Trapping terms
- Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz)
- Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz)
- Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz)
-
- Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1
+ Unapm1j = znapm1j (ip,ij) * interp_napj(ip-1,ij ,iky,ikx,iz)
+ Unapm1jp1 = znapm1jp1(ip,ij) * interp_napj(ip-1,ij+1,iky,ikx,iz)
+ Unapm1jm1 = znapm1jm1(ip,ij) * interp_napj(ip-1,ij-1,iky,ikx,iz)
+ Tmir = Tnapp1j + Tnapp1jm1 + Tnapm1j + Tnapm1jm1 + Unapm1j + Unapm1jp1 + Unapm1jm1
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) &
- + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) &
- + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz
+ Tphi = (xphij (ip,ij)*kernel(ij ,iky,ikx,iz,eo) &
+ + xphijp1(ip,ij)*kernel(ij+1,iky,ikx,iz,eo) &
+ + xphijm1(ip,ij)*kernel(ij-1,iky,ikx,iz,eo))*phikykxz
ELSE
Tphi = 0._dp
ENDIF
!! Vector potential term
IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3
- Tpsi = (xpsij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) &
- + xpsijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) &
- + xpsijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*psikykxz
+ Tpsi = (xpsij (ip,ij)*kernel(ij ,iky,ikx,iz,eo) &
+ + xpsijp1(ip,ij)*kernel(ij+1,iky,ikx,iz,eo) &
+ + xpsijm1(ip,ij)*kernel(ij-1,iky,ikx,iz,eo))*psikykxz
ELSE
Tpsi = 0._dp
ENDIF
-
!! Sum of all RHS terms
- moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
+ moments_rhs_(ip,ij,iky,ikx,iz) = &
! Perpendicular magnetic gradient/curvature effects
- -imagu*Ckxky(iky,ikx,iz,eo) * Tperp &
+ -imagu*Ckxky(iky,ikx,iz,eo) * Tperp&
! Perpendicular pressure effects
-i_ky*beta*dpdx * Tperp&
- ! Parallel coupling (Landau damping)
+ ! Parallel coupling (Landau Damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
- -dBdz(iz,eo)*gradz_coeff(iz,eo) * Tmir &
+ -dBdz(iz,eo)*gradz_coeff(iz,eo) * Tmir&
! Drives (density + temperature gradients)
-i_ky * (Tphi - Tpsi) &
- ! Numerical hyperdiffusion (totally artificial, for stability purpose)
- -mu_x*diff_kx_coeff*kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
- -mu_y*diff_ky_coeff*ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
- ! Numerical parallel hyperdiffusion "mu_z*ddz**4"
- -mu_z*diff_dz_coeff*ddzND_nipj(ip,ij,iky,ikx,iz) &
+ ! Numerical Hermite hyperdiffusion (GX version)
+ -mu_p*diff_pe_coeff*p_int**4*moments_(ip,ij,iky,ikx,iz)&
+ ! Numerical Laguerre hyperdiffusion (GX version)
+ -mu_j*diff_je_coeff*j_int**4*moments_(ip,ij,iky,ikx,iz)&
+ ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
+ -mu_x*diff_kx_coeff*kx**N_HD*moments_(ip,ij,iky,ikx,iz) &
+ -mu_y*diff_ky_coeff*ky**N_HD*moments_(ip,ij,iky,ikx,iz) &
+ ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
+ -mu_z*diff_dz_coeff*ddzND_napj(ip,ij,iky,ikx,iz) &
! Collision term
- +TColl_i(ip,ij,iky,ikx,iz)&
- ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor
- -hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz)
+ +TColl_(ip,ij,iky,ikx,iz) &
+ ! Nonlinear term
+ -hatB_NL(iz,eo) * Sapj(ip,ij,iky,ikx,iz)
- ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
- ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! ! (not used often so not parallelized)
- ! moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
- ! moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) &
- ! + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel)
- ELSE
- moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
- ENDIF
- END DO ploopi
- END DO jloopi
- END DO kyloopi
- END DO kxloopi
- END DO zloopi
+ ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
+ ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
+ ! ! (not used often so not parallelized)
+ ! moments_rhs_(ip,ij,iky,ikx,iz) = &
+ ! moments_rhs_(ip,ij,iky,ikx,iz) &
+ ! + mu_p * moments_(ip-4,ij,iky,ikx,iz)
- ! Execution time end
- CALL cpu_time(t1_rhs)
- tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
+ ELSE
+ moments_rhs_(ip,ij,iky,ikx,iz) = 0._dp
+ ENDIF
+ END DO ploop
+ END DO jloop
+ END DO kyloop
+ END DO kxloop
+ END DO zloop
+
+ ! Execution time end
+ CALL cpu_time(t1_rhs)
+ tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
-END SUBROUTINE moments_eq_rhs_i
+ END SUBROUTINE moments_eq_rhs
+ !_____________________________________________________________________________!
+ !_____________________________________________________________________________!
+
+END SUBROUTINE compute_moments_eq_rhs
SUBROUTINE add_Maxwellian_background_terms
! This routine is meant to add the terms rising from the magnetic operator,
--
GitLab