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

multiple typos correction

parent 2246519c
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
+ 91
− 86
View file @ c7b37e0c
...@@ -12,74 +12,79 @@ SUBROUTINE moments_eq_rhs ...@@ -12,74 +12,79 @@ SUBROUTINE moments_eq_rhs
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
COMPLEX(dp) :: kr, kz, kperp2 REAL(dp) :: kr, kz, kperp2
COMPLEX(dp) :: taue_qe_etaBi, taui_qi_etaBi REAL(dp) :: taue_qe_etaB, taui_qi_etaB
COMPLEX(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
COMPLEX(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables REAL(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables
COMPLEX(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
COMPLEX(dp) :: xphij, xphijp1, xphijm1, xColl REAL(dp) :: xphij, xphijp1, xphijm1, xColl
COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi, TColl ! terms of the rhs COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi, TColl ! terms of the rhs
!!!!!!!!! Electrons moments RHS !!!!!!!!! !Precompute species dependant factors
Tphi = 0 ! electrostatic potential term taue_qe_etaB = tau_e/q_e * eta_B
taue_qe_etaBi = tau_e/q_e * eta_B * imagu ! one-time computable factor xb_e2 = sigma_e**2 * tau_e!/2.0 ! species dependant factor of the Kernel, squared
xb_e2 = sigma_e**2 * tau_e/2. ! species dependant factor of the Kernel, squared taui_qi_etaB = tau_i/q_i * eta_B
xb_i2 = sigma_i**2 * tau_i!/2.0 ! species dependant factor of the Kernel, squared
ploope : DO ip = ips_e, ipe_e !!!!!!!!! Electrons moments RHS !!!!!!!!!
ip_dp = real(ip-1,dp) ! real index is one minus the loop index (0 to pmax) Tphi = 0 ! electrostatic potential term
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)
! x N_e^{p+2,j} ! N_e^{p+2,j} multiplicator
IF (ip+2 .LE. pmaxe + 1) THEN IF (ip+2 .LE. ipe_e) THEN
xNapp2j = taue_qe_etaBi * SQRT(ip_dp + 1._dp) * SQRT(ip_dp + 2._dp) xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE ELSE
xNapp2j = 0. xNapp2j = 0.
ENDIF ENDIF
! x N_e^{p-2,j} ! N_e^{p-2,j} multiplicator
IF (ip-2 .GE. 1) THEN IF (ip-2 .GE. ips_e) THEN
xNapm2j = taue_qe_etaBi * SQRT(ip_dp) * SQRT(ip_dp - 1.) xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.))
ELSE ELSE
xNapm2j = 0. xNapm2j = 0.
ENDIF ENDIF
jloope : DO ij = ijs_e, ije_e 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)
! x N_e^{p,j+1} ! N_e^{p,j+1} multiplicator
IF (ij+1 .LE. jmaxe+1) THEN IF (ij+1 .LE. ije_e) THEN
xNapjp1 = -taue_qe_etaBi * (ij_dp + 1.) xNapjp1 = taue_qe_etaB * (ij_dp + 1.)
ELSE ELSE
xNapjp1 = 0. xNapjp1 = 0.
ENDIF ENDIF
! x N_e^{p,j-1} ! N_e^{p,j-1} multiplicator
IF (ij-1 .GE. 1) THEN IF (ij-1 .GE. ijs_e) THEN
xNapjm1 = -taue_qe_etaBi * ij_dp xNapjm1 = taue_qe_etaB * ij_dp
ELSE ELSE
xNapjm1 = 0. xNapjm1 = 0.
ENDIF ENDIF
! x N_e^{pj} ! N_e^{pj} multiplicator
xNapj = taue_qe_etaBi * 2.*(ip_dp + ij_dp + 1.) xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)
! first part of collision operator (Lenard-Bernstein) ! first part of collision operator (Lenard-Bernstein)
xColl = -nu*(ip_dp + 2.*ij_dp) xColl = -(ip_dp + 2.*ij_dp)
! x phi ! phi multiplicator for different kernel numbers
IF (ip .eq. 1) THEN !(krokecker delta_p^0) IF (ip .EQ. 1) THEN !(kronecker delta_p^0)
xphij = imagu * (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 = -imagu * (eta_B + eta_T)*(ij_dp+1.) xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
xphijm1 = -imagu * (eta_B + eta_T)* ij_dp xphijm1 = -(eta_T - eta_B)* ij_dp
factj = real(Factorial(ij-1),dp) factj = REAL(Factorial(ij-1),dp)
ELSE IF (ip .eq. 3) THEN !(krokecker delta_p^2) ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
xphij = imagu * (SQRT2*eta_B + eta_T/SQRT2) xphij = (eta_T/SQRT2 - SQRT2*eta_B)
factj = real(Factorial(ij-1),dp) factj = REAL(Factorial(ij-1),dp)
ELSE ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0. xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
factj = 1 factj = 1
ENDIF ENDIF
! Loop on kspace !write(*,*) '(ip,ij) = (', ip,',', ij,')'
! 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
...@@ -91,33 +96,33 @@ SUBROUTINE moments_eq_rhs ...@@ -91,33 +96,33 @@ SUBROUTINE moments_eq_rhs
! 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 (xNapp2j .NE. (0.,0.)) THEN IF (ip+2 .LE. ipe_e) 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 (xNapm2j .NE. (0.,0.)) THEN IF (ip-2 .GE. ips_e) 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 (xNapjp1 .NE. (0.,0.)) THEN IF (ij+1 .LE. ije_e) THEN
TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapp2j 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 (xNapjm1 .NE. (0.,0.)) THEN IF (ij-1 .GE. ijs_e) THEN
TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapp2j TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE ELSE
TNapjm1 = 0. TNapjm1 = 0.
ENDIF ENDIF
! Collision term completed (Lenard-Bernstein) ! Collision term completed (Lenard-Bernstein)
IF (xNapp2j .NE. (0.,0.)) THEN IF (ip+2 .LE. ipe_e) THEN
TColl = (xColl - nu * b_e2/4.) * moments_e(ip+2,ij,ikr,ikz,updatetlevel) TColl = nu*(xColl - b_e2) * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
ELSE ELSE
TColl = 0. TColl = 0.
ENDIF ENDIF
...@@ -128,13 +133,12 @@ SUBROUTINE moments_eq_rhs ...@@ -128,13 +133,12 @@ SUBROUTINE moments_eq_rhs
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)* kz * phi(ikr,ikz) Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
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) = &
kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1) & imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
+ Tphi + TColl
END DO kzloope END DO kzloope
END DO krloope END DO krloope
...@@ -142,65 +146,64 @@ SUBROUTINE moments_eq_rhs ...@@ -142,65 +146,64 @@ SUBROUTINE moments_eq_rhs
END DO jloope END DO jloope
END DO ploope END DO ploope
!!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!!
Tphi = 0 ! electrostatic potential term Tphi = 0 ! electrostatic potential term
taui_qi_etaBi = tau_i/real(q_i)*eta_B * imagu ! one-time computable factor
xb_i2 = sigma_i**2 * tau_i/2.0 ! species dependant factor of the Kernel, squared
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. pmaxi + 1) THEN IF (ip+2 .LE. ipe_i) THEN
xNapp2j = taui_qi_etaBi * SQRT(ip_dp + 1.) * SQRT(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. 1) THEN IF (ip-2 .GE. ips_i) THEN
xNapm2j = taui_qi_etaBi * SQRT(ip_dp) * SQRT(ip_dp - 1.) xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.))
ELSE ELSE
xNapm2j = 0. xNapm2j = 0.
ENDIF ENDIF
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)
! x N_i^{p,j+1} ! x N_i^{p,j+1}
IF (ij+1 .LE. jmaxi+1) THEN IF (ij+1 .LE. ije_i) THEN
xNapjp1 = -taui_qi_etaBi * (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. 1) THEN IF (ij-1 .GE. ijs_i) THEN
xNapjm1 = -taui_qi_etaBi * ij_dp xNapjm1 = taui_qi_etaB * ij_dp
ELSE ELSE
xNapjm1 = 0. xNapjm1 = 0.
ENDIF ENDIF
! x N_i^{pj} ! x N_i^{pj}
xNapj = taui_qi_etaBi * 2.*(ip_dp + ij_dp + 1.) xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)
! first part of collision operator (Lenard-Bernstein) ! first part of collision operator (Lenard-Bernstein)
xColl = -nu*(ip_dp + 2.0*ij_dp) xColl = -(ip_dp + 2.0*ij_dp)
! x phi ! x phi
IF (ip .EQ. 1) THEN !(krokecker delta_p^0) IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
xphij = imagu * (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 = -imagu * (eta_B + eta_T)*(ij_dp+1.) xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
xphijm1 = -imagu * (eta_B + eta_T)* ij_dp xphijm1 = -(eta_T - eta_B)* ij_dp
factj = REAL(Factorial(ij-1),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 = imagu * (SQRT2*eta_B + eta_T/SQRT2) xphij = (eta_T/SQRT2 - SQRT2*eta_B)
factj = REAL(Factorial(ij-1),dp) factj = REAL(Factorial(ij-1),dp)
ELSE ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0. xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
factj = 1.
ENDIF ENDIF
! Loop on kspace
! Loop on kspace
krloopi : DO ikr = ikrs,ikre krloopi : DO ikr = ikrs,ikre
kzloopi : DO ikz = ikzs,ikze kzloopi : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector kr = krarray(ikr) ! Poloidal wavevector
...@@ -211,31 +214,34 @@ SUBROUTINE moments_eq_rhs ...@@ -211,31 +214,34 @@ SUBROUTINE moments_eq_rhs
!! 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}
IF (xNapp2j .NE. (0.,0.)) THEN ! term propto N_i^{p+2,j} IF (ip+2 .LE. ipe_i) 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
IF (xNapm2j .NE. (0.,0.)) THEN ! term propto N_i^{p-2,j} ! term propto N_i^{p-2,j}
IF (ip-2 .GE. ips_i) 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
IF (xNapjp1 .NE. (0.,0.)) THEN ! xterm propto N_a^{p,j+1} ! xterm propto N_i^{p,j+1}
TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapp2j IF (ij+1 .LE. ije_i) THEN
TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE ELSE
TNapjp1 = 0. TNapjp1 = 0.
ENDIF ENDIF
IF (xNapjm1 .NE. (0.,0.)) THEN ! term propto N_a^{p,j-1} ! term propto N_i^{p,j-1}
TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapp2j IF (ij-1 .GE. ijs_i) THEN
TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE ELSE
TNapjm1 = 0. TNapjm1 = 0.
ENDIF ENDIF
! Collision term completed (Lenard-Bernstein) ! Collision term completed (Lenard-Bernstein)
IF (xNapp2j .NE. (0.,0.)) THEN IF (ip+2 .LE. ipe_i) THEN
TColl = (xColl - nu * b_i2/4.) * moments_i(ip+2,ij,ikr,ikz,updatetlevel) TColl = nu*(xColl - b_i2) * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
ELSE ELSE
TColl = 0. TColl = 0.
ENDIF ENDIF
...@@ -244,15 +250,14 @@ SUBROUTINE moments_eq_rhs ...@@ -244,15 +250,14 @@ SUBROUTINE moments_eq_rhs
Tphi = 0 Tphi = 0
IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker 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._dp) 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)* kz * phi(ikr,ikz) Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
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) = &
kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1) & imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
+ Tphi + TColl
END DO kzloopi END DO kzloopi
END DO krloopi END DO krloopi
......
src/poisson.F90
+ 45
− 39
View file @ c7b37e0c
...@@ -12,16 +12,19 @@ SUBROUTINE poisson ...@@ -12,16 +12,19 @@ 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
COMPLEX(dp) :: kr, kz, kperp2 REAL(dp) :: kr, kz, kperp2
COMPLEX(dp) :: xkm_, xk2_ ! sub kernel factor for recursive build REAL(dp) :: Kn, Kn2 ! sub kernel factor for recursive build
COMPLEX(dp) :: b_e2, b_i2, alphaD REAL(dp) :: b_e2, b_i2, alphaD
COMPLEX(dp) :: sigmaetaue2, sigmaitaui2 ! To avoid redondant computation REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation
COMPLEX(dp) :: gammaD, gammaphi COMPLEX(dp) :: gammaD, gammaphi
COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i ! Store to compute sum Kn^2
COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store to compute sum Kn*Napn
sigmaetaue2 = sigma_e**2 * tau_e/2. !Precompute species dependant factors
sigmaitaui2 = sigma_i**2 * tau_i/2. sigmae2_taue_o2 = sigma_e**2 * tau_e/2.
sigmai2_taui_o2 = sigma_i**2 * tau_i/2.
qe2_taue = (q_e**2)/tau_e
qi2_taui = (q_i**2)/tau_i
DO ikr=ikrs,ikre DO ikr=ikrs,ikre
DO ikz=ikzs,ikze DO ikz=ikzs,ikze
...@@ -29,57 +32,60 @@ SUBROUTINE poisson ...@@ -29,57 +32,60 @@ SUBROUTINE poisson
kz = kzarray(ikz) kz = kzarray(ikz)
kperp2 = kr**2 + kz**2 kperp2 = kr**2 + kz**2
! Compute electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2) !! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
b_e2 = kperp2 * sigmaetaue2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a)/2) b_e2 = kperp2 * sigmae2_taue_o2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a))
xkm_ = 1. ! Initialization for n = 0 ! Initialization for n = 0 (ine = 1)
xk2_ = 1. Kn = 1.
Kn2 = 1.
sum_kernel_mom_e = moments_e(1,1,ikr,ikz,updatetlevel) sum_kernel_mom_e = moments_e(1,1,ikr,ikz,updatetlevel)
sum_kernel2_e = xk2_ sum_kernel2_e = Kn2
! 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=2,jmaxe+1 DO ine=1,jmaxe+1 ! ine = n+1
ine_dp = REAL(ine-1,dp) ine_dp = REAL(ine-1,dp)
xkm_ = xkm_ * b_e2/ine_dp Kn = Kn * b_e2/(ine_dp+1) ! We update iteratively the kernel functions (to spare factorial computations)
xk2_ = xk2_ *(b_e2/ine_dp)**2 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 + xkm_ * moments_e(1,ine,ikr,ikz,updatetlevel) sum_kernel_mom_e = sum_kernel_mom_e + Kn * moments_e(1,ine,ikr,ikz,updatetlevel)
sum_kernel2_e = sum_kernel2_e + xk2_ sum_kernel2_e = sum_kernel2_e + Kn2 ! ... sum recursively ...
END DO END DO
endif endif
sum_kernel2_e = sum_kernel2_e * exp(-b_e2) sum_kernel_mom_e = sum_kernel_mom_e * exp(-b_e2) ! ... multiply the final sum with the missing exponential term
sum_kernel_mom_e = sum_kernel_mom_e * exp(-2.*b_e2) sum_kernel2_e = sum_kernel2_e * exp(-2.*b_e2) ! and its squared using exp(x)^2 = exp(2x).
! Compute ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2) !! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
b_i2 = kperp2 * sigmaitaui2 b_i2 = kperp2 * sigmai2_taui_o2
xkm_ = 1.
xk2_ = 1.
! Initialization for n = 0 (ini = 1)
Kn = 1.
Kn2 = 1.
sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel) sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel)
sum_kernel2_i = xk2_ sum_kernel2_i = Kn2
! 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=2,jmaxi + 1 DO ini=1,jmaxi+1
ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax) ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
xkm_ = xkm_ * b_i2/ini_dp Kn = Kn * b_i2/(ini_dp+1) ! We update iteratively the kernel functions this time for ions ...
xk2_ = xk2_ *(b_i2/ini_dp)**2 Kn2 = Kn2 *(b_i2/(ini_dp+1.))**2
sum_kernel_mom_i = sum_kernel_mom_i + xkm_ * moments_i(1,ini,ikr,ikz,updatetlevel) sum_kernel_mom_i = sum_kernel_mom_i + Kn * moments_i(1,ini,ikr,ikz,updatetlevel)
sum_kernel2_i = sum_kernel2_i + xk2_ sum_kernel2_i = sum_kernel2_i + Kn2 ! ... sum recursively ...
END DO END DO
endif endif
sum_kernel2_i = sum_kernel2_i * exp(-b_i2) sum_kernel_mom_i = sum_kernel_mom_i * exp(-b_i2) ! ... multiply the final sum with the missing exponential term.
sum_kernel_mom_i = sum_kernel_mom_i * exp(-2.*b_i2) 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 + (q_e**2)/tau_e * sum_kernel2_e & gammaD = alphaD + qe2_taue * (1. - sum_kernel2_e) &
+ (q_i**2)/tau_i * sum_kernel2_i + qi2_taui * (1. - sum_kernel2_i)
gammaphi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i 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) = gammaphi/gammaD
......
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