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

Typos and presentation

parent 9e2c5064
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
+ 95
− 93
View file @ f94fe388
...@@ -6,153 +6,154 @@ SUBROUTINE moments_eq_rhs ...@@ -6,153 +6,154 @@ SUBROUTINE moments_eq_rhs
USE fields USE fields
USE fourier_grid USE fourier_grid
USE model USE model
USE prec_const
use prec_const
IMPLICIT NONE IMPLICIT NONE
INTEGER :: ip,ij, ikr,ikz ! loops indices INTEGER :: ip,ij, ikr,ikz ! loops indices
REAL(dp) :: ip_dp, ij_dp REAL(dp) :: ip_dp, ij_dp
REAL(dp) :: kr, kz, kperp2 REAL(dp) :: kr, kz, kperp2
REAL(dp) :: taue_qe_etaB, taui_qi_etaB REAL(dp) :: taue_qe_etaB, taui_qi_etaB
REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
REAL(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables
REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
REAL(dp) :: xphij, xphijp1, xphijm1, xColl REAL(dp) :: xphij, xphijp1, xphijm1, xColl
COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
!Precompute species dependant factors !Precompute species dependant factors
taue_qe_etaB = tau_e/q_e * eta_B taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
xb_e2 = sigma_e**2 * tau_e!/2.0 ! species dependant factor of the Kernel, squared taui_qi_etaB = tau_i/q_i * eta_B
taui_qi_etaB = tau_i/q_i * eta_B sigmae2_taue_o2 = sigma_e**2 * tau_e/2.0 ! factor of the Kernel argument
xb_i2 = sigma_i**2 * tau_i!/2.0 ! species dependant factor of the Kernel, squared sigmai2_taui_o2 = sigma_i**2 * tau_i/2.0
!!!!!!!!! Electrons moments RHS !!!!!!!!! !!!!!!!!! Electrons moments RHS !!!!!!!!!
Tphi = 0 ! electrostatic potential term
ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1 ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1
ip_dp = REAL(ip-1.,dp) ! REAL index is one minus the loop index (0 to pmax) ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmax)
! N_e^{p+2,j} multiplicator ! N_e^{p+2,j} multiplicator
IF (ip+2 .LE. ipe_e) THEN IF (ip+2 .LE. pmaxe+1) THEN
xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.)) xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE ELSE
xNapp2j = 0. xNapp2j = 0.
ENDIF ENDIF
! N_e^{p-2,j} multiplicator ! N_e^{p-2,j} multiplicator
IF (ip-2 .GE. ips_e) THEN IF (ip-2 .GE. 1) THEN
xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.)) xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.))
ELSE ELSE
xNapm2j = 0. xNapm2j = 0.
ENDIF ENDIF
factj = 1.0 ! Start of the recursive factorial
jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1 jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1
ij_dp = REAL(ij-1.,dp) ! REAL index is one minus the loop index (0 to jmax) ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmax)
IF (ij_dp .GT. 0) THEN
factj = factj * ij_dp; ! Recursive factorial
ENDIF
! N_e^{p,j+1} multiplicator ! N_e^{p,j+1} multiplicator
IF (ij+1 .LE. ije_e) THEN IF (ij+1 .LE. jmaxe+1) THEN
xNapjp1 = taue_qe_etaB * (ij_dp + 1.) xNapjp1 = taue_qe_etaB * (ij_dp + 1.)
ELSE ELSE
xNapjp1 = 0. xNapjp1 = 0.
ENDIF ENDIF
! N_e^{p,j-1} multiplicator ! N_e^{p,j-1} multiplicator
IF (ij-1 .GE. ijs_e) THEN IF (ij-1 .GE. 1) THEN
xNapjm1 = taue_qe_etaB * ij_dp xNapjm1 = taue_qe_etaB * ij_dp
ELSE ELSE
xNapjm1 = 0. xNapjm1 = 0.
ENDIF ENDIF
! N_e^{pj} multiplicator ! N_e^{pj} multiplicator
xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.) xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)
! Collision operator (DK Lenard-Bernstein basis) ! Collision operator (DK Lenard-Bernstein basis)
xColl = ip_dp + 2.*ij_dp xColl = ip_dp + 2.*ij_dp
! ... adding Dougherty terms
IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) ) THEN ! kronecker p0 * j1
TColl01 = 2.0/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
TColl20 = 0.0; TColl10 = 0.0;
ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) ) THEN ! kronecker p2 * j0
TColl20 = -SQRT2/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
TColl10 = 0.0; TColl01 = 0.0;
ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
TColl10 = moments_e(2,1,ikr,ikz,updatetlevel)
TColl20 = 0.0; TColl01 = 0.0;
ELSE
TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
ENDIF
! phi multiplicator for different kernel numbers ! phi multiplicator for different kernel numbers
IF (ip .EQ. 1) THEN !(kronecker delta_p^0) IF (ip .EQ. 1) THEN !(kronecker delta_p^0)
xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) ) xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
xphijp1 = -(eta_T - eta_B)*(ij_dp+1.) xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
xphijm1 = -(eta_T - eta_B)* ij_dp xphijm1 = -(eta_T - eta_B)* ij_dp
factj = REAL(Factorial(ij-1),dp)
ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2) ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
xphij = (eta_T/SQRT2 - SQRT2*eta_B) xphij = (eta_T/SQRT2 - SQRT2*eta_B)
factj = REAL(Factorial(ij-1),dp) xphijp1 = 0.; xphijm1 = 0.
ELSE ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0. xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
factj = 1
ENDIF ENDIF
!write(*,*) '(ip,ij) = (', ip,',', ij,')'
! Loop on kspace ! Loop on kspace
krloope : DO ikr = ikrs,ikre krloope : DO ikr = ikrs,ikre
kzloope : DO ikz = ikzs,ikze kzloope : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector kz = kzarray(ikz) ! Toroidal wavevector
kperp2 = kr**2 + kz**2 ! perpendicular wavevector kperp2 = kr**2 + kz**2 ! perpendicular wavevector
b_e2 = kperp2 * xb_e2 ! Bessel argument b_e2 = kperp2 * sigmae2_taue_o2 ! Bessel argument
!! Compute moments and mixing terms !! Compute moments and mixing terms
! term propto N_e^{p,j} ! term propto N_e^{p,j}
TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj
! term propto N_e^{p+2,j} ! term propto N_e^{p+2,j}
IF (ip+2 .LE. ipe_e) THEN IF (ip+2 .LE. pmaxe+1) THEN
TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
ELSE ELSE
TNapp2j = 0. TNapp2j = 0.
ENDIF ENDIF
! term propto N_e^{p-2,j} ! term propto N_e^{p-2,j}
IF (ip-2 .GE. ips_e) THEN IF (ip-2 .GE. 1) THEN
TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
ELSE ELSE
TNapm2j = 0. TNapm2j = 0.
ENDIF ENDIF
! xterm propto N_e^{p,j+1} ! xterm propto N_e^{p,j+1}
IF (ij+1 .LE. ije_e) THEN IF (ij+1 .LE. jmaxe+1) THEN
TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1 TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE ELSE
TNapjp1 = 0. TNapjp1 = 0.
ENDIF ENDIF
! term propto N_e^{p,j-1} ! term propto N_e^{p,j-1}
IF (ij-1 .GE. ijs_e) THEN IF (ij-1 .GE. 1) THEN
TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1 TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE ELSE
TNapjm1 = 0. TNapjm1 = 0.
ENDIF ENDIF
! Collision term completed (DK Dougherty) ! Dougherty Collision terms
IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) .AND. (pmaxe .GE. 2) ) THEN ! kronecker p0 * j1
TColl01 = 2.0/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
TColl20 = 0.0; TColl10 = 0.0;
ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) .AND. (jmaxe .GE. 1)) THEN ! kronecker p2 * j0
TColl20 = -SQRT2/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
TColl10 = 0.0; TColl01 = 0.0;
ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
TColl10 = moments_e(2,1,ikr,ikz,updatetlevel)
TColl20 = 0.0; TColl01 = 0.0;
ELSE
TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
ENDIF
! Total collisional term
TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) & TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
+ TColl01 + TColl10 + TColl20) + TColl01 + TColl10 + TColl20)
!! Electrical potential term !! Electrical potential term
Tphi = 0 IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker delta_p^0, delta_p^2
IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0)
kernelj = b_e2**(ij-1) * exp(-b_e2)/factj kernelj = b_e2**(ij-1) * exp(-b_e2)/factj
kerneljp1 = kernelj * b_e2 /(ij_dp + 1.) kerneljp1 = kernelj * b_e2 /(ij_dp + 1.)
kerneljm1 = kernelj * ij_dp / b_e2 kerneljm1 = kernelj * ij_dp / b_e2
Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz) Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0
ENDIF ENDIF
! Sum of all precomputed terms ! Sum of all precomputed terms
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = & moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
END DO kzloope END DO kzloope
END DO krloope END DO krloope
...@@ -160,37 +161,41 @@ SUBROUTINE moments_eq_rhs ...@@ -160,37 +161,41 @@ SUBROUTINE moments_eq_rhs
END DO ploope END DO ploope
!!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!!
Tphi = 0 ! electrostatic potential term
ploopi : DO ip = ips_i, ipe_i ploopi : DO ip = ips_i, ipe_i
ip_dp = REAL(ip-1.,dp) ip_dp = REAL(ip-1,dp)
! x N_i^{p+2,j} ! x N_i^{p+2,j}
IF (ip+2 .LE. ipe_i) THEN IF (ip+2 .LE. pmaxi+1) THEN
xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.)) xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE ELSE
xNapp2j = 0. xNapp2j = 0.
ENDIF ENDIF
! x N_i^{p-2,j} ! x N_i^{p-2,j}
IF (ip-2 .GE. ips_i) THEN IF (ip-2 .GE. 1) THEN
xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.)) xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.))
ELSE ELSE
xNapm2j = 0. xNapm2j = 0.
ENDIF ENDIF
factj = 1.0 ! Start of the recursive factorial
jloopi : DO ij = ijs_i, ije_i jloopi : DO ij = ijs_i, ije_i
ij_dp = REAL(ij-1.,dp) ij_dp = REAL(ij-1,dp)
IF (ij_dp .GT. 0) THEN
factj = factj * ij_dp; ! Recursive factorial
ENDIF
! x N_i^{p,j+1} ! x N_i^{p,j+1}
IF (ij+1 .LE. ije_i) THEN IF (ij+1 .LE. jmaxi+1) THEN
xNapjp1 = taui_qi_etaB * (ij_dp + 1.) xNapjp1 = taui_qi_etaB * (ij_dp + 1.)
ELSE ELSE
xNapjp1 = 0. xNapjp1 = 0.
ENDIF ENDIF
! x N_i^{p,j-1} ! x N_i^{p,j-1}
IF (ij-1 .GE. ijs_i) THEN IF (ij-1 .GE. 1) THEN
xNapjm1 = taui_qi_etaB * ij_dp xNapjm1 = taui_qi_etaB * ij_dp
ELSE ELSE
xNapjm1 = 0. xNapjm1 = 0.
...@@ -199,38 +204,19 @@ SUBROUTINE moments_eq_rhs ...@@ -199,38 +204,19 @@ SUBROUTINE moments_eq_rhs
! x N_i^{pj} ! x N_i^{pj}
xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.) xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)
! Collision term completed (DK Dougherty) ! Collision operator (DK Lenard-Bernstein basis)
TColl = -nu * (xColl * moments_i(ip,ij,ikr,ikz,updatetlevel) & xColl = ip_dp + 2.*ij_dp
+ TColl01 + TColl10 + TColl20)
! ... adding Dougherty terms
IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) ) THEN ! kronecker p0 * j1
TColl01 = 2.0/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
TColl20 = 0.0; TColl10 = 0.0;
ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) ) THEN ! kronecker p2 * j0
TColl20 = -SQRT2/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
TColl10 = 0.0; TColl01 = 0.0;
ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
TColl10 = moments_i(2,1,ikr,ikz,updatetlevel)
TColl20 = 0.0; TColl01 = 0.0;
ELSE
TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
ENDIF
! x phi ! x phi
IF (ip .EQ. 1) THEN !(krokecker delta_p^0) IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) ) xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
xphijp1 = -(eta_T - eta_B)*(ij_dp+1.) xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
xphijm1 = -(eta_T - eta_B)* ij_dp xphijm1 = -(eta_T - eta_B)* ij_dp
factj = REAL(Factorial(ij-1),dp)
ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2) ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
xphij = (eta_T/SQRT2 - SQRT2*eta_B) xphij = (eta_T/SQRT2 - SQRT2*eta_B)
factj = REAL(Factorial(ij-1),dp) xphijp1 = 0.; xphijm1 = 0.
ELSE ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0. xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
factj = 1.
ENDIF ENDIF
! Loop on kspace ! Loop on kspace
...@@ -239,52 +225,68 @@ SUBROUTINE moments_eq_rhs ...@@ -239,52 +225,68 @@ SUBROUTINE moments_eq_rhs
kr = krarray(ikr) ! Poloidal wavevector kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector kz = kzarray(ikz) ! Toroidal wavevector
kperp2 = kr**2 + kz**2 ! perpendicular wavevector kperp2 = kr**2 + kz**2 ! perpendicular wavevector
b_i2 = kperp2 * xb_i2 ! Bessel argument b_i2 = kperp2 * sigmai2_taui_o2 ! Bessel argument
!! Compute moments and mixing terms !! Compute moments and mixing terms
! term propto N_i^{p,j} ! term propto N_i^{p,j}
TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj
! term propto N_i^{p+2,j} ! term propto N_i^{p+2,j}
IF (ip+2 .LE. ipe_i) THEN IF (ip+2 .LE. pmaxi+1) THEN
TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
ELSE ELSE
TNapp2j = 0. TNapp2j = 0.
ENDIF ENDIF
! term propto N_i^{p-2,j} ! term propto N_i^{p-2,j}
IF (ip-2 .GE. ips_i) THEN IF (ip-2 .GE. 1) THEN
TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
ELSE ELSE
TNapm2j = 0. TNapm2j = 0.
ENDIF ENDIF
! xterm propto N_i^{p,j+1} ! xterm propto N_i^{p,j+1}
IF (ij+1 .LE. ije_i) THEN IF (ij+1 .LE. jmaxi+1) THEN
TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1 TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE ELSE
TNapjp1 = 0. TNapjp1 = 0.
ENDIF ENDIF
! term propto N_i^{p,j-1} ! term propto N_i^{p,j-1}
IF (ij-1 .GE. ijs_i) THEN IF (ij-1 .GE. 1) THEN
TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1 TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE ELSE
TNapjm1 = 0. TNapjm1 = 0.
ENDIF ENDIF
! Collision term completed (Dougherty) ! Dougherty Collision terms
TColl = -nu* (xColl * moments_i(ip,ij,ikr,ikz,updatetlevel)) IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) .AND. (pmaxi .GE. 2)) THEN ! kronecker p0 * j1
TColl01 = 2.0/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
TColl20 = 0.0; TColl10 = 0.0;
ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) .AND. (jmaxi .GE. 1)) THEN ! kronecker p2 * j0
TColl20 = -SQRT2/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
- 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
TColl10 = 0.0; TColl01 = 0.0;
ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
TColl10 = moments_i(2,1,ikr,ikz,updatetlevel)
TColl20 = 0.0; TColl01 = 0.0;
ELSE
TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
ENDIF
! Total collisional term
TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
+ TColl01 + TColl10 + TColl20)
!! Electrical potential term !! Electrical potential term
Tphi = 0 IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! kronecker delta_p^0, delta_p^2
IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0, delta_p^2)
kernelj = b_i2**(ij-1) * exp(-b_i2)/factj kernelj = b_i2**(ij-1) * exp(-b_i2)/factj
kerneljp1 = kernelj * b_i2 /(ij_dp + 1.) kerneljp1 = kernelj * b_i2 /(ij_dp + 1.)
kerneljm1 = kernelj * ij_dp / b_i2 kerneljm1 = kernelj * ij_dp / b_i2
Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz) Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0
ENDIF ENDIF
! Sum of all precomputed terms ! Sum of all precomputed terms
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = & moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
END DO kzloopi END DO kzloopi
END DO krloopi END DO krloopi
......
src/poisson.F90
+ 28
− 36
View file @ f94fe388
...@@ -13,82 +13,74 @@ SUBROUTINE poisson ...@@ -13,82 +13,74 @@ SUBROUTINE poisson
INTEGER :: ikr,ikz, ini,ine, i0j INTEGER :: ikr,ikz, ini,ine, i0j
REAL(dp) :: ini_dp, ine_dp REAL(dp) :: ini_dp, ine_dp
REAL(dp) :: kr, kz, kperp2 REAL(dp) :: kr, kz, kperp2
REAL(dp) :: Kn, Kn2 ! sub kernel factor for recursive build REAL(dp) :: Kne, Kni ! sub kernel factor for recursive build
REAL(dp) :: b_e2, b_i2, alphaD REAL(dp) :: b_e2, b_i2, alphaD
REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation
COMPLEX(dp) :: gammaD, gammaphi REAL(dp) :: sum_kernel2_e, sum_kernel2_i ! Store sum Kn^2
COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i ! Store to compute sum Kn^2 COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store sum Kn*Napn
COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store to compute sum Kn*Napn REAL(dp) :: gammaD
COMPLEX(dp) :: gammaD_phi
!Precompute species dependant factors !Precompute species dependant factors
sigmae2_taue_o2 = sigma_e**2 * tau_e/2. sigmae2_taue_o2 = sigma_e**2 * tau_e/2. ! factor of the Kernel argument
sigmai2_taui_o2 = sigma_i**2 * tau_i/2. sigmai2_taui_o2 = sigma_i**2 * tau_i/2.
qe2_taue = (q_e**2)/tau_e qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum
qi2_taui = (q_i**2)/tau_i qi2_taui = (q_i**2)/tau_i
DO ikr=ikrs,ikre DO ikr=ikrs,ikre
DO ikz=ikzs,ikze DO ikz=ikzs,ikze
kr = krarray(ikr) kr = krarray(ikr)
kz = kzarray(ikz) kz = kzarray(ikz)
kperp2 = kr**2 + kz**2 kperp2 = kr**2 + kz**2
!! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2) !! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
b_e2 = kperp2 * sigmae2_taue_o2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a)) b_e2 = kperp2 * sigmae2_taue_o2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a))
! Initialization for n = 0 (ine = 1) ! Initialization for n = 0 (ine = 1)
Kn = 1. Kne = exp(-b_e2)
Kn2 = 1. sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
sum_kernel_mom_e = moments_e(1,1,ikr,ikz,updatetlevel) sum_kernel2_e = Kne**2
sum_kernel2_e = Kn2
! loop over n only if the max polynomial degree is 1 or more ! loop over n only if the max polynomial degree is 1 or more
if (jmaxe .GT. 0) then if (jmaxe .GT. 0) then
DO ine=1,jmaxe+1 ! ine = n+1 DO ine=2,jmaxe+1 ! ine = n+1
ine_dp = REAL(ine-1,dp) ine_dp = REAL(ine-1,dp)
Kn = Kn * b_e2/(ine_dp+1) ! We update iteratively the kernel functions (to spare factorial computations) Kne = Kne * b_e2/ine_dp ! We update iteratively the kernel functions (to spare factorial computations)
Kn2 = Kn2 *(b_e2/(ine_dp+1))**2 ! the exp term is still missing but does not depend on ine ...
sum_kernel_mom_e = sum_kernel_mom_e + Kn * moments_e(1,ine,ikr,ikz,updatetlevel) sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
sum_kernel2_e = sum_kernel2_e + Kn2 ! ... sum recursively ... sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
END DO END DO
endif endif
sum_kernel_mom_e = sum_kernel_mom_e * exp(-b_e2) ! ... multiply the final sum with the missing exponential term
sum_kernel2_e = sum_kernel2_e * exp(-2.*b_e2) ! and its squared using exp(x)^2 = exp(2x).
!! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2) !! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
b_i2 = kperp2 * sigmai2_taui_o2 b_i2 = kperp2 * sigmai2_taui_o2
! Initialization for n = 0 (ini = 1) ! Initialization for n = 0 (ini = 1)
Kn = 1. Kni = exp(-b_i2)
Kn2 = 1. sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel)
sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel) sum_kernel2_i = Kni**2
sum_kernel2_i = Kn2
! loop over n only if the max polynomial degree is 1 or more ! loop over n only if the max polynomial degree is 1 or more
if (jmaxi .GT. 0) then if (jmaxi .GT. 0) then
DO ini=1,jmaxi+1 DO ini=2,jmaxi+1
ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax) ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
Kn = Kn * b_i2/(ini_dp+1) ! We update iteratively the kernel functions this time for ions ... Kni = Kni * b_i2/ini_dp ! We update iteratively the kernel functions this time for ions ...
Kn2 = Kn2 *(b_i2/(ini_dp+1.))**2
sum_kernel_mom_i = sum_kernel_mom_i + Kn * moments_i(1,ini,ikr,ikz,updatetlevel) sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
sum_kernel2_i = sum_kernel2_i + Kn2 ! ... sum recursively ... sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ...
END DO END DO
endif endif
sum_kernel_mom_i = sum_kernel_mom_i * exp(-b_i2) ! ... multiply the final sum with the missing exponential term.
sum_kernel2_i = sum_kernel2_i * exp(-2.*b_i2)
!!! Assembling the poisson equation !!! Assembling the poisson equation
alphaD = kperp2 * lambdaD**2. alphaD = kperp2 * lambdaD**2
gammaD = alphaD + qe2_taue * (1. - sum_kernel2_e) & gammaD = alphaD + qe2_taue * (1. - sum_kernel2_e) &
+ qi2_taui * (1. - sum_kernel2_i) + qi2_taui * (1. - sum_kernel2_i)
gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i ! gamma_D * phi term
gammaphi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i ! gamma_D * phi term
phi(ikr, ikz) = gammaphi/gammaD phi(ikr, ikz) = gammaD_phi/gammaD
END DO END DO
END DO END DO
......
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