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

added maxwellian background term for test purpose

parent 37685bde
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
+ 110
− 1
View file @ 46c46efc
......@@ -8,6 +8,8 @@ SUBROUTINE compute_moments_eq_rhs
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
!_____________________________________________________________________________!
!_____________________________________________________________________________!
......@@ -109,12 +111,19 @@ SUBROUTINE moments_eq_rhs_e
- i_ky * Tphi &
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" see Pueschel 2010 (eq 25)
+ mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) &
! Collision term
+ TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term
- Sepj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
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
......@@ -238,6 +247,12 @@ SUBROUTINE moments_eq_rhs_i
+ TColl_i(ip,ij,iky,ikx,iz)&
! Nonlinear term
- Sipj(ip,ij,iky,ikx,iz)
IF(ip-4 .GT. 0) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
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
......@@ -253,4 +268,98 @@ SUBROUTINE moments_eq_rhs_i
END SUBROUTINE moments_eq_rhs_i
SUBROUTINE add_Maxwellian_background_terms
! This routine is meant to add the terms rising from the magnetic operator,
! i.e. (B x GradB) Grad, applied on the background Maxwellian distribution
! (x_a + spar^2)(b x GradB) GradFaM
! It gives birth to kx=ky=0 sources terms (averages) that hit moments 00, 20,
! 40, 01,02, 21 with background gradient dependences.
USE prec_const
USE time_integration, ONLY : updatetlevel
USE model, ONLY: taue_qe, taui_qi, K_N, K_T, KIN_E
USE array, ONLY: moments_rhs_e, moments_rhs_i
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,&
jarray_e,parray_e,jarray_i,parray_i, zarray, izs,ize,&
ip,ij
IMPLICIT NONE
real(dp), DIMENSION(izs:ize) :: sinz
sinz(izs:ize) = SIN(zarray(izs:ize,0))
IF(contains_kx0 .AND. contains_ky0) THEN
IF(KIN_E) THEN
DO ip = ips_e,ipe_e
DO ij = ijs_e,ije_e
SELECT CASE(ij-1)
CASE(0) ! j = 0
SELECT CASE (ip-1)
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)&
+taue_qe * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T)
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)&
+taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T)
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)&
+taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*K_T
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
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)&
-taue_qe * sinz(izs:ize) * (K_N + 3.5_dp*K_T)
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)&
-taue_qe * sinz(izs:ize) * SQRT2*K_T
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
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)&
+taue_qe * sinz(izs:ize) * 2._dp*K_T
END SELECT
END SELECT
ENDDO
ENDDO
ENDIF
DO ip = ips_i,ipe_i
DO ij = ijs_i,ije_i
SELECT CASE(ij-1)
CASE(0) ! j = 0
SELECT CASE (ip-1)
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)&
+taui_qi * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T)
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)&
+taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T)
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)&
+taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*K_T
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
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)&
-taui_qi * sinz(izs:ize) * (K_N + 3.5_dp*K_T)
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)&
-taui_qi * sinz(izs:ize) * SQRT2*K_T
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
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)&
+taui_qi * sinz(izs:ize) * 2._dp*K_T
END SELECT
END SELECT
ENDDO
ENDDO
ENDIF
END SUBROUTINE
END MODULE moments_eq_rhs
src/prec_const_mod.F90
+ 2
− 1
View file @ 46c46efc
......@@ -22,12 +22,13 @@ MODULE prec_const
INTEGER, private :: MPI_MIN_DP !Min reduction operation for DP datatype
! Some useful constants, to avoid recomputing then too often
! Some useful constants, to avoid recomputing them too often
REAL(dp), PARAMETER :: PI=3.141592653589793238462643383279502884197_dp
REAL(dp), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_dp
REAL(dp), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_dp
REAL(dp), PARAMETER :: ONEOPI=0.3183098861837906912164442019275156781077_dp
REAL(dp), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_dp
REAL(dp), PARAMETER :: SQRT6=SQRT(6._dp)
REAL(dp), PARAMETER :: INVSQRT2=0.7071067811865475244008443621048490392848359377_dp
REAL(dp), PARAMETER :: SQRT3=1.73205080756887729352744634150587236694281_dp
REAL(dp), PARAMETER :: onetwelfth=0.08333333333333333333333333333333333333333333333_dp
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment