Skip to content
Snippets Groups Projects
Commit ce66464e authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

EM implementation and small rewriting

parent 6796aee9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/moments_eq_rhs_mod.F90
+ 55
− 41
View file @ ce66464e
...@@ -6,6 +6,7 @@ CONTAINS ...@@ -6,6 +6,7 @@ CONTAINS
SUBROUTINE compute_moments_eq_rhs SUBROUTINE compute_moments_eq_rhs
USE model, only: KIN_E USE model, only: KIN_E
IMPLICIT NONE IMPLICIT NONE
IF(KIN_E) CALL moments_eq_rhs_e IF(KIN_E) CALL moments_eq_rhs_e
CALL moments_eq_rhs_i CALL moments_eq_rhs_i
!! BETA TESTING !! !! BETA TESTING !!
...@@ -21,13 +22,18 @@ SUBROUTINE moments_eq_rhs_e ...@@ -21,13 +22,18 @@ SUBROUTINE moments_eq_rhs_e
USE basic USE basic
USE time_integration USE time_integration
USE array USE array, ONLY: xnepj, xnepp2j, xnepm2j, xnepjp1, xnepjm1, xnepp1j, xnepm1j,&
USE fields 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 grid
USE model USE model
USE prec_const USE prec_const
USE collision USE collision
use geometry USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR
USE calculus, ONLY : interp_z, grad_z, grad_z2 USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE IMPLICIT NONE
...@@ -77,6 +83,8 @@ SUBROUTINE moments_eq_rhs_e ...@@ -77,6 +83,8 @@ SUBROUTINE moments_eq_rhs_e
Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz) Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz)
! term propto n_e^{p,j-1} ! term propto n_e^{p,j-1}
Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz) Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz)
! Tperp
Tperp = Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1
! Parallel dynamic ! Parallel dynamic
! ddz derivative for Landau damping term ! ddz derivative for Landau damping term
Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) & Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) &
...@@ -102,7 +110,7 @@ SUBROUTINE moments_eq_rhs_e ...@@ -102,7 +110,7 @@ SUBROUTINE moments_eq_rhs_e
ENDIF ENDIF
!! Vector potential term !! Vector potential term
IF ( (p_int .LE. 3) .AND. (p_int .GT. 1) ) THEN ! Kronecker p1 or p3 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) & Tpsi = (xpsij_e (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
+ xpsijp1_e(ip,ij)*kernel_e(ij+1,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 + xpsijm1_e(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*psikykxz
...@@ -113,22 +121,22 @@ SUBROUTINE moments_eq_rhs_e ...@@ -113,22 +121,22 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all RHS terms !! Sum of all RHS terms
moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = & moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects ! Perpendicular magnetic gradient/curvature effects
- imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)& -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp&
! Parallel coupling (Landau Damping) ! Parallel coupling (Landau Damping)
- Tpar*gradz_coeff(iz,eo) & -gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient) ! Mirror term (parallel magnetic gradient)
- gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) & -gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir&
! Drives (density + temperature gradients) ! Drives (density + temperature gradients)
- i_ky * (Tphi - Tpsi) & -i_ky * (Tphi - Tpsi) &
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) ! 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_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) & -mu_y*diff_ky_coeff*ky**N_HD*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" see Pueschel 2010 (eq 25) ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
- mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) & -mu_z*diff_dz_coeff*ddzND_nepj(ip,ij,iky,ikx,iz) &
! Collision term ! Collision term
+ TColl_e(ip,ij,iky,ikx,iz) & +TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term ! Nonlinear term
- hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz) -hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) & IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
...@@ -160,13 +168,18 @@ SUBROUTINE moments_eq_rhs_i ...@@ -160,13 +168,18 @@ SUBROUTINE moments_eq_rhs_i
USE basic USE basic
USE time_integration, ONLY: updatetlevel USE time_integration, ONLY: updatetlevel
USE array USE array, ONLY: xnipj, xnipp2j, xnipm2j, xnipjp1, xnipjm1, xnipp1j, xnipm1j,&
USE fields 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 grid USE grid
USE model USE model
USE prec_const USE prec_const
USE collision USE collision
USE geometry USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR
USE calculus, ONLY : interp_z, grad_z, grad_z2 USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE IMPLICIT NONE
...@@ -201,6 +214,7 @@ SUBROUTINE moments_eq_rhs_i ...@@ -201,6 +214,7 @@ SUBROUTINE moments_eq_rhs_i
p_int = parray_i(ip) ! Hermite degree p_int = parray_i(ip) ! Hermite degree
eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
kperp2= kparray(iky,ikx,iz,eo)**2 kperp2= kparray(iky,ikx,iz,eo)**2
IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN
!! Compute moments mixing terms !! Compute moments mixing terms
Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
...@@ -243,7 +257,7 @@ SUBROUTINE moments_eq_rhs_i ...@@ -243,7 +257,7 @@ SUBROUTINE moments_eq_rhs_i
ENDIF ENDIF
!! Vector potential term !! Vector potential term
IF ( (p_int .LE. 3) .AND. (p_int .GT. 1) ) THEN ! Kronecker p1 or p3 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) & Tpsi = (xpsij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) &
+ xpsijp1_i(ip,ij)*kernel_i(ij+1,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 + xpsijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*psikykxz
...@@ -255,22 +269,22 @@ SUBROUTINE moments_eq_rhs_i ...@@ -255,22 +269,22 @@ SUBROUTINE moments_eq_rhs_i
!! Sum of all RHS terms !! Sum of all RHS terms
moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = & moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects ! Perpendicular magnetic gradient/curvature effects
- imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp & -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp &
! Parallel coupling (Landau damping) ! Parallel coupling (Landau damping)
- gradz_coeff(iz,eo) * Tpar & -gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient) ! Mirror term (parallel magnetic gradient)
- gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir & -gradzB(iz,eo)*gradz_coeff(iz,eo) * Tmir &
! Drives (density + temperature gradients) ! Drives (density + temperature gradients)
- i_ky * (Tphi - Tpsi) & -i_ky * (Tphi - Tpsi) &
! Numerical hyperdiffusion (totally artificial, for stability purpose) ! Numerical hyperdiffusion (totally artificial, for stability purpose)
- mu_x*diff_kx_coeff*kx**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) & -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) & -mu_y*diff_ky_coeff*ky**N_HD*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" ! Numerical parallel hyperdiffusion "mu_z*ddz**4"
+ mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) & -mu_z*diff_dz_coeff*ddzND_nipj(ip,ij,iky,ikx,iz) &
! Collision term ! Collision term
+ TColl_i(ip,ij,iky,ikx,iz)& +TColl_i(ip,ij,iky,ikx,iz)&
! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor ! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor
- hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz) -hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) & IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33) ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
...@@ -300,7 +314,7 @@ SUBROUTINE add_Maxwellian_background_terms ...@@ -300,7 +314,7 @@ SUBROUTINE add_Maxwellian_background_terms
! 40, 01,02, 21 with background gradient dependences. ! 40, 01,02, 21 with background gradient dependences.
USE prec_const USE prec_const
USE time_integration, ONLY : updatetlevel USE time_integration, ONLY : updatetlevel
USE model, ONLY: taue_qe, taui_qi, k_N, eta_N, k_T, eta_T, KIN_E USE model, ONLY: taue_qe, taui_qi, k_Ni, k_Ne, k_Ti, k_Te, KIN_E
USE array, ONLY: moments_rhs_e, moments_rhs_i USE array, ONLY: moments_rhs_e, moments_rhs_i
USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,& USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,&
ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,& ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,&
...@@ -320,28 +334,28 @@ SUBROUTINE add_Maxwellian_background_terms ...@@ -320,28 +334,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1) SELECT CASE (ip-1)
CASE(0) ! Na00 term CASE(0) ! Na00 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taue_qe * sinz(izs:ize) * (1.5_dp*eta_N*k_N - 1.125_dp*eta_N*k_T) +taue_qe * sinz(izs:ize) * (1.5_dp*k_Ne - 1.125_dp*k_Te)
CASE(2) ! Na20 term CASE(2) ! Na20 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*eta_N*k_N - 2.75_dp*eta_T*k_T) +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*k_Ne - 2.75_dp*k_Te)
CASE(4) ! Na40 term CASE(4) ! Na40 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*eta_T*k_T +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*k_Te
END SELECT END SELECT
CASE(1) ! j = 1 CASE(1) ! j = 1
SELECT CASE (ip-1) SELECT CASE (ip-1)
CASE(0) ! Na01 term CASE(0) ! Na01 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
-taue_qe * sinz(izs:ize) * (eta_N*k_N + 3.5_dp*eta_T*k_T) -taue_qe * sinz(izs:ize) * (k_Ne + 3.5_dp*k_Te)
CASE(2) ! Na21 term CASE(2) ! Na21 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
-taue_qe * sinz(izs:ize) * SQRT2*eta_T*k_T -taue_qe * sinz(izs:ize) * SQRT2*k_Te
END SELECT END SELECT
CASE(2) ! j = 2 CASE(2) ! j = 2
SELECT CASE (ip-1) SELECT CASE (ip-1)
CASE(0) ! Na02 term CASE(0) ! Na02 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taue_qe * sinz(izs:ize) * 2._dp*eta_T*k_T +taue_qe * sinz(izs:ize) * 2._dp*k_Te
END SELECT END SELECT
END SELECT END SELECT
ENDDO ENDDO
...@@ -355,28 +369,28 @@ SUBROUTINE add_Maxwellian_background_terms ...@@ -355,28 +369,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1) SELECT CASE (ip-1)
CASE(0) ! Na00 term CASE(0) ! Na00 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taui_qi * sinz(izs:ize) * (1.5_dp*k_N - 1.125_dp*k_T) +taui_qi * sinz(izs:ize) * (1.5_dp*k_Ni - 1.125_dp*k_Ti)
CASE(2) ! Na20 term CASE(2) ! Na20 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_N - 2.75_dp*k_T) +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_Ni - 2.75_dp*k_Ti)
CASE(4) ! Na40 term CASE(4) ! Na40 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_T +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_Ti
END SELECT END SELECT
CASE(1) ! j = 1 CASE(1) ! j = 1
SELECT CASE (ip-1) SELECT CASE (ip-1)
CASE(0) ! Na01 term CASE(0) ! Na01 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
-taui_qi * sinz(izs:ize) * (k_N + 3.5_dp*k_T) -taui_qi * sinz(izs:ize) * (k_Ni + 3.5_dp*k_Ti)
CASE(2) ! Na21 term CASE(2) ! Na21 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
-taui_qi * sinz(izs:ize) * SQRT2*k_T -taui_qi * sinz(izs:ize) * SQRT2*k_Ti
END SELECT END SELECT
CASE(2) ! j = 2 CASE(2) ! j = 2
SELECT CASE (ip-1) SELECT CASE (ip-1)
CASE(0) ! Na02 term CASE(0) ! Na02 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)& moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
+taui_qi * sinz(izs:ize) * 2._dp*k_T +taui_qi * sinz(izs:ize) * 2._dp*k_Ti
END SELECT END SELECT
END SELECT END SELECT
ENDDO ENDDO
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment