Newer
Older
MODULE moments_eq_rhs
IMPLICIT NONE
PUBLIC :: compute_moments_eq_rhs
SUBROUTINE compute_moments_eq_rhs
USE model
USE array
USE fields
USE grid, ONLY: local_na, local_np, local_nj, local_nkx, local_nky, local_nz,&
nzgrid,pp2,ngp,ngj,ngz,&
diff_dz_coeff,diff_kx_coeff,diff_ky_coeff,diff_p_coeff,diff_j_coeff,&
parray,jarray,kxarray, kyarray, kp2array
USE basic
USE closure, ONLY: evolve_mom
USE prec_const
USE collision
USE time_integration
USE species, ONLY: Ptot
USE geometry, ONLY: gradz_coeff, dlnBdz, Ckxky!, Gamma_phipar
USE calculus, ONLY: interp_z, grad_z, grad_z2
IMPLICIT NONE
INTEGER :: ia, iz, iky, ikx, ip ,ij, eo ! counters
INTEGER :: izi,ipi,iji ! interior points counters
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
REAL(xp) :: kx, ky, kperp2
COMPLEX(xp) :: Tnapj, Tnapp2j, Tnapm2j, Tnapjp1, Tnapjm1 ! Terms from b x gradB and drives
COMPLEX(xp) :: Tnapp1j, Tnapm1j, Tnapp1jm1, Tnapm1jm1 ! Terms from mirror force with non adiab moments_
COMPLEX(xp) :: Ldamp, Fmir
COMPLEX(xp) :: Mperp, Mpara, Dphi, Dpsi
COMPLEX(xp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
COMPLEX(xp) :: i_kx,i_ky
COMPLEX(xp) :: Napj, RHS, phikykxz, psikykxz
! Spatial loops
z:DO iz = 1,local_nz
izi = iz + ngz/2
x:DO ikx = 1,local_nkx
y:DO iky = 1,local_nky
kx = kxarray(iky,ikx) ! radial wavevector
ky = kyarray(iky) ! binormal wavevector
i_kx = imagu * kx ! radial derivative
i_ky = imagu * ky ! binormal derivative
phikykxz = phi(iky,ikx,izi)
psikykxz = psi(iky,ikx,izi)
! Kinetic loops
j:DO ij = 1,local_nj ! This loop is from 1 to jmaxi+1
iji = ij+ngj/2
j_int = jarray(iji)
p:DO ip = 1,local_np ! Hermite loop
ipi = ip+ngp/2
p_int = parray(ipi) ! Hermite degree
eo = min(nzgrid,MODULO(p_int,2)+1) ! Indicates if we are on odd or even z grid
kperp2= kp2array(iky,ikx,izi,eo)
! Species loop
a:DO ia = 1,local_na
Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
RHS = 0._xp
IF(evolve_mom(ipi,iji)) THEN ! for the closure scheme
!! Compute moments_ mixing terms
! Perpendicular dynamic
! term propto n^{p,j}
Tnapj = xnapj(ia,ip,ij)* nadiab_moments(ia,ipi, iji, iky,ikx,izi)
! term propto n^{p+2,j}
Tnapp2j = xnapp2j(ia,ip) * nadiab_moments(ia,ipi+pp2,iji, iky,ikx,izi)
! term propto n^{p-2,j}
Tnapm2j = xnapm2j(ia,ip) * nadiab_moments(ia,ipi-pp2,iji, iky,ikx,izi)
! term propto n^{p,j+1}
Tnapjp1 = xnapjp1(ia,ij) * nadiab_moments(ia,ipi, iji+1,iky,ikx,izi)
! term propto n^{p,j-1}
Tnapjm1 = xnapjm1(ia,ij) * nadiab_moments(ia,ipi, iji-1,iky,ikx,izi)
! Perpendicular magnetic term (curvature, gradient drifts and alpha MHD pressure drift)
Mperp = imagu*(Ckxky(iky,ikx,izi,eo) - beta*Ptot*ky)&
*(Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)
! Parallel dynamic
! ddz derivative for Landau damping term
Ldamp = xnapp1j(ia,ip) * ddz_napj(ia,ipi+1,iji,iky,ikx,iz) &
+ xnapm1j(ia,ip) * ddz_napj(ia,ipi-1,iji,iky,ikx,iz)
Antoine Cyril David Hoffmann
committed
! neglect Landau damping in temperature equation and higher moments as in Ivanov 2022
IF(RM_LD_T_EQ .AND. ((j_int .GT. 0) .OR. (p_int .GT. 1))) Ldamp = 0._xp
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
! Mirror terms
Tnapp1j = ynapp1j (ia,ip,ij) * interp_napj(ia,ipi+1,iji ,iky,ikx,iz)
Tnapp1jm1 = ynapp1jm1(ia,ip,ij) * interp_napj(ia,ipi+1,iji-1,iky,ikx,iz)
Tnapm1j = ynapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
Tnapm1jm1 = ynapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
! Trapping terms
Unapm1j = znapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
Unapm1jp1 = znapm1jp1(ia,ip,ij) * interp_napj(ia,ipi-1,iji+1,iky,ikx,iz)
Unapm1jm1 = znapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
! sum the parallel forces
Fmir = dlnBdz(izi,eo)*(Tnapp1j + Tnapp1jm1 + Tnapm1j + Tnapm1jm1 +&
Unapm1j + Unapm1jp1 + Unapm1jm1)
! Parallel magnetic term (Landau damping and the mirror force)
Mpara = gradz_coeff(izi,eo)*(Ldamp + Fmir)
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
Dphi =i_ky*( xphij(ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
+xphijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
+xphijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo) )*phikykxz
ELSE
Dphi = 0._xp
ENDIF
!! Vector potential term
IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3
Dpsi =-i_ky*( xpsij (ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
+xpsijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
+xpsijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo))*psikykxz
ELSE
Dpsi = 0._xp
ENDIF
!! Sum of all RHS terms
RHS = &
! Nonlinear term Sapj_ = {phi,f}
- Sapj(ia,ip,ij,iky,ikx,iz) &
! Perpendicular magnetic term
- Mperp &
! Parallel magnetic term
- Mpara &
! Drives (density + temperature gradients)
- (Dphi + Dpsi) &
! Collision term
+ Capj(ia,ip,ij,iky,ikx,iz) &
! Parallel drive term (should be negligible, to test)
! !! -Gamma_phipar(iz,eo)*Tphi*ddz_phi(iky,ikx,iz) &
! Numerical perpendicular hyperdiffusion
-mu_x*diff_kx_coeff*kx**N_HD*Napj &
-mu_y*diff_ky_coeff*ky**N_HD*Napj &
! ! ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
-mu_z*diff_dz_coeff*ddzND_Napj(ia,ipi,iji,iky,ikx,iz)
!! Velocity space dissipation (should be implemented somewhere else)
SELECT CASE(HYP_V)
CASE('hypcoll') ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
IF (p_int .GT. 2) &
RHS = RHS - mu_p*diff_p_coeff*p_int**6*Napj
IF (j_int .GT. 1) &
RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
CASE('dvpar4')
! fourth order numerical diffusion in vpar
IF(p_int .GE. 4) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
! (not used often so not parallelized)
RHS = RHS + mu_p*dv4_Hp_coeff(p_int)*moments(ia,ipi-4,iji,iky,ikx,izi,updatetlevel)
! + dummy Laguerre diff
IF (j_int .GE. 2) &
RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
CASE DEFAULT
END SELECT
RHS = 0._xp
!! Put RHS in the array
moments_rhs(ia,ip,ij,iky,ikx,iz,updatetlevel) = RHS
END DO a
END DO p
END DO j
END DO y
END DO x
END DO z
END SUBROUTINE compute_moments_eq_rhs
END MODULE moments_eq_rhs