Skip to content
Snippets Groups Projects
Commit 5a4c7c3c authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

Split moment_eq_rhs in two routines (i,e)

parent bb9d5767
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/moments_eq_rhs.F90
+ 80
− 144
View file @ 5a4c7c3c
SUBROUTINE moments_eq_rhs
!_____________________________________________________________________________!
!_____________________________________________________________________________!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!! Electrons moments RHS !!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!_____________________________________________________________________________!
SUBROUTINE moments_eq_rhs_e
USE basic
USE time_integration
......@@ -20,14 +26,11 @@ SUBROUTINE moments_eq_rhs
REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
COMPLEX(dp) :: test_nan, i_kz
COMPLEX(dp) :: i_kz
! Measuring execution time
CALL cpu_time(t0_rhs)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! Electrons moments RHS !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices
p_int= parray_e(ip) ! Hermite polynom. degree
p_dp = REAL(p_int,dp) ! REAL of the Hermite degree
......@@ -42,8 +45,6 @@ SUBROUTINE moments_eq_rhs
! N_e^{p-2,j} coeff
xNapm2j = taue_qe_etaB * SQRT(p_dp * (p_dp - 1._dp))
! factj = 1.0 ! Start of the recursive factorial
jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices
j_int= jarray_e(ij) ! Laguerre polynom. degree
j_dp = REAL(j_int,dp) ! REAL of degree
......@@ -92,13 +93,6 @@ SUBROUTINE moments_eq_rhs
xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp
ENDIF
! Recursive factorial
! IF (j_dp .GT. 0) THEN
! factj = factj * j_dp
! ELSE
! factj = 1._dp
! ENDIF
! Loop on kspace
krloope : DO ikr = ikrs,ikre
kzloope : DO ikz = ikzs,ikze
......@@ -113,69 +107,36 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_e^{p,j}
TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
TNapp1j = 0._dp; TNapm1j = 0._dp
TNapp2j = 0._dp; TNapm2j = 0._dp
TNapjp1 = 0._dp; TNapjm1 = 0._dp
! term propto N_e^{p+1,j} and kparallel
IF ( (p_int+1 .LE. pmaxe) ) THEN ! OoB check
TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
ELSE
TNapp1j = 0._dp
ENDIF
IF ( p_int+1 .LE. pmaxe ) TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p-1,j} and kparallel
IF ( (p_int-1 .GE. 0) ) THEN ! OoB check
TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
ELSE
TNapm1j = 0._dp
ENDIF
IF ( p_int-1 .GE. 0 ) TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p+2,j}
IF (p_int+2 .LE. pmaxe) THEN ! OoB check
TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
TNapp2j = 0._dp
ENDIF
IF ( p_int+2 .LE. pmaxe ) TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
! term propto N_e^{p-2,j}
IF (p_int-2 .GE. 0) THEN ! OoB check
TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
ELSE
TNapm2j = 0._dp
ENDIF
IF ( p_int-2 .GE. 0 ) TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
! xterm propto N_e^{p,j+1}
IF (j_int+1 .LE. jmaxe) THEN ! OoB check
TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
ELSE
TNapjp1 = 0._dp
ENDIF
IF ( j_int+1 .LE. jmaxe ) TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
! term propto N_e^{p,j-1}
IF (j_int-1 .GE. 0) THEN ! OoB check
TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
ELSE
TNapjm1 = 0._dp
ENDIF
IF ( j_int-1 .GE. 0 ) TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
!! Collision
IF (CO .EQ. -3) THEN ! GK Dougherty
CALL DoughertyGK_e(ip,ij,ikr,ikz,TColl)
ELSEIF (CO .EQ. -2) THEN ! DK Dougherty
IF ( (pmaxe .GE. 2) ) THEN ! OoB check
TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
ELSE
TColl20 = 0._dp
ENDIF
IF ( (jmaxe .GE. 1) ) THEN ! OoB check
TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
ELSE
TColl01 = 0._dp
ENDIF
IF ( (pmaxe .GE. 1) ) THEN ! OoB + odd number for Hermite degrees check
TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
ELSE
TColl10 = 0._dp
ENDIF
TColl20 = 0._dp; TColl01 = 0._dp; TColl10 = 0._dp
IF ( (pmaxe .GE. 2) ) TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
IF ( (jmaxe .GE. 1) ) TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
IF ( (pmaxe .GE. 1) ) TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
! Total collisional term
TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
ELSEIF (CO .EQ. -1) THEN ! Full Coulomb for electrons (COSOlver matrix)
TColl = 0._dp ! Initialization
ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms
......@@ -198,23 +159,15 @@ SUBROUTINE moments_eq_rhs
ELSEIF (CO .EQ. 0) THEN ! Lenard Bernstein
TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
kernelj = kernel_e(ij, ikr, ikz)
IF ( j_int+1 .LE. jmaxe ) THEN
kerneljp1 = kernel_e(ij+1, ikr, ikz)
ELSE
kerneljp1 = 0._dp
ENDIF
IF ( j_int-1 .GE. 0 ) THEN
kerneljm1 = kernel_e(ij-1, ikr, ikz)
ELSE
kerneljm1 = 0._dp
ENDIF
kerneljp1 = 0._dp; kerneljm1 = 0._dp
IF ( j_int+1 .LE. jmaxe ) kerneljp1 = kernel_e(ij+1, ikr, ikz)
IF ( j_int-1 .GE. 0 ) kerneljm1 = kernel_e(ij-1, ikr, ikz)
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
......@@ -239,12 +192,44 @@ SUBROUTINE moments_eq_rhs
END DO jloope
END DO ploope
!_____________________________________________________________________________!
!_____________________________________________________________________________!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!! Ions moments RHS !!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!_____________________________________________________________________________!
! 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
USE fields
USE grid
USE model
USE prec_const
USE utility, ONLY : is_nan
USE collision
IMPLICIT NONE
INTEGER :: ip2, ij2, p_int, j_int, p2_int, j2_int ! loops indices and polynom. degrees
REAL(dp) :: p_dp, j_dp
REAL(dp) :: kr, kz, kperp2
REAL(dp) :: kernelj, kerneljp1, kerneljm1 ! Kernel functions and variable
REAL(dp) :: xNapj, xNapp1j, xNapm1j, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop
REAL(dp) :: xphij, xphijp1, xphijm1, xphijpar ! ESpot. factors depending on the pj loop
REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
COMPLEX(dp) :: i_kz
! Measuring execution time
CALL cpu_time(t0_rhs)
ploopi : DO ip = ips_i, ipe_i ! Hermite loop
p_int= parray_i(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
......@@ -253,14 +238,11 @@ SUBROUTINE moments_eq_rhs
xNapp1j = sqrtTaui_qi * SQRT(p_dp + 1)
! N_i^{p-1,j} coeff
xNapm1j = sqrtTaui_qi * SQRT(p_dp)
! x N_i^{p+2,j} coeff
xNapp2j = taui_qi_etaB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
! x N_i^{p-2,j} coeff
xNapm2j = taui_qi_etaB * SQRT(p_dp * (p_dp - 1._dp))
! factj = 1._dp ! Start of the recursive factorial
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
j_int= jarray_i(ij) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
......@@ -309,13 +291,6 @@ SUBROUTINE moments_eq_rhs
xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar = 0._dp
ENDIF
! ! Recursive factorial
! IF (j_dp .GT. 0) THEN
! factj = factj * j_dp
! ELSE
! factj = 1._dp
! ENDIF
! Loop on kspace
krloopi : DO ikr = ikrs,ikre
kzloopi : DO ikz = ikzs,ikze
......@@ -328,63 +303,31 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_i^{p,j}
TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
TNapp1j = 0._dp; TNapm1j = 0._dp
TNapp2j = 0._dp; TNapm2j = 0._dp
TNapjp1 = 0._dp; TNapjm1 = 0._dp
! term propto N_i^{p+1,j}
IF ( (p_int+1 .LE. pmaxi) ) THEN ! OoB check
TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
ELSE
TNapp1j = 0._dp
ENDIF
IF ( p_int+1 .LE. pmaxi ) TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p-1,j}
IF ( (p_int-1 .GE. 0) ) THEN ! OoB check
TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
ELSE
TNapm1j = 0._dp
ENDIF
IF ( p_int-1 .GE. 0 ) TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p+2,j}
IF (p_int+2 .LE. pmaxi) THEN ! OoB check
TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
TNapp2j = 0._dp
ENDIF
IF ( p_int+2 .LE. pmaxi ) TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
! term propto N_i^{p-2,j}
IF (p_int-2 .GE. 0) THEN ! OoB check
TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
ELSE
TNapm2j = 0._dp
ENDIF
IF ( p_int-2 .GE. 0 ) TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
! xterm propto N_i^{p,j+1}
IF (j_int+1 .LE. jmaxi) THEN ! OoB check
TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
ELSE
TNapjp1 = 0._dp
ENDIF
IF ( j_int+1 .LE. jmaxi ) TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
! term propto N_i^{p,j-1}
IF (j_int-1 .GE. 0) THEN ! OoB check
TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
ELSE
TNapjm1 = 0._dp
ENDIF
IF ( j_int-1 .GE. 0 ) TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
!! Collision
IF (CO .EQ. -3) THEN ! Gyrokin. Dougherty Collision terms
CALL DoughertyGK_i(ip,ij,ikr,ikz,TColl)
ELSEIF (CO .EQ. -2) THEN ! Dougherty Collision terms
IF ( (pmaxi .GE. 2) ) THEN ! OoB check
TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
ELSE
TColl20 = 0._dp
ENDIF
IF ( (jmaxi .GE. 1) ) THEN ! OoB check
TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
ELSE
TColl01 = 0._dp
ENDIF
IF ( (pmaxi .GE. 1) ) THEN ! OoB check
TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
ELSE
TColl10 = 0._dp
ENDIF
TColl20 = 0._dp; TColl01 = 0._dp; TColl10 = 0._dp
IF ( (pmaxe .GE. 2) ) TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
IF ( (jmaxe .GE. 1) ) TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
IF ( (pmaxe .GE. 1) ) TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
! Total collisional term
TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
......@@ -417,16 +360,9 @@ SUBROUTINE moments_eq_rhs
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
kernelj = kernel_i(ij, ikr, ikz)
IF (j_int+1 .LE. jmaxi) THEN
kerneljp1 = kernel_i(ij+1, ikr, ikz)
ELSE
kerneljp1 = 0.0_dp
ENDIF
IF ( j_int-1 .GE. 0 ) THEN
kerneljm1 = kernel_i(ij-1, ikr, ikz)
ELSE
kerneljm1 = 0._dp
ENDIF
kerneljp1 = 0._dp; kerneljm1 = 0._dp
IF ( j_int+1 .LE. jmaxe ) kerneljp1 = kernel_i(ij+1, ikr, ikz)
IF ( j_int-1 .GE. 0 ) kerneljm1 = kernel_i(ij-1, ikr, ikz)
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
......@@ -455,4 +391,4 @@ SUBROUTINE moments_eq_rhs
CALL cpu_time(t1_rhs)
tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
END SUBROUTINE moments_eq_rhs
END SUBROUTINE moments_eq_rhs_i
src/stepon.F90
+ 2
− 1
View file @ 5a4c7c3c
......@@ -17,7 +17,8 @@ SUBROUTINE stepon
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
! Compute right hand side of moments hierarchy equation
CALL moments_eq_rhs
CALL moments_eq_rhs_e
CALL moments_eq_rhs_i
! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
CALL advance_time_level
......
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