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

linear coeff computation is generalized to any species

parent f4f7fae1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 133 additions and 166 deletions
src/numerics_mod.F90
+ 133
− 166
View file @ 2ec5c035
...@@ -217,178 +217,145 @@ END SUBROUTINE evaluate_ampere_op ...@@ -217,178 +217,145 @@ END SUBROUTINE evaluate_ampere_op
!******************************************************************************! !******************************************************************************!
SUBROUTINE compute_lin_coeff SUBROUTINE compute_lin_coeff
USE array
USE model, ONLY: taue_qe, taui_qi, & USE array, ONLY: xnepj, &
k_Te, k_Ti, k_Ne, k_Ni, CurvB, GradB, KIN_E,& ynepp1j, ynepm1j, ynepp1jm1, ynepm1jm1,&
tau_e, tau_i, sigma_e, sigma_i zNepm1j, zNepm1jp1, zNepm1jm1,&
xnepp1j, xnepm1j, xnepp2j, xnepm2j,&
xnepjp1, xnepjm1,&
xphij_e, xphijp1_e, xphijm1_e,&
xpsij_e, xpsijp1_e, xpsijm1_e,&
xnipj, &
ynipp1j, ynipm1j, ynipp1jm1, ynipm1jm1,&
zNipm1j, zNipm1jp1, zNipm1jm1,&
xnipp1j, xnipm1j, xnipp2j, xnipm2j,&
xnipjp1, xnipjm1,&
xphij_i, xphijp1_i, xphijm1_i,&
xpsij_i, xpsijp1_i, xpsijm1_i
USE model, ONLY: k_Te, k_Ti, k_Ne, k_Ni, k_cB, k_gB, KIN_E,&
tau_e, tau_i, sigma_e, sigma_i, q_e, q_i
USE prec_const USE prec_const
USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, & USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i
IMPLICIT NONE
INTEGER :: p_int, j_int ! polynom. degrees IF(KIN_E) THEN
REAL(dp) :: p_dp, j_dp CALL lin_coeff(k_Te,k_Ne,k_cB,k_gB,tau_e,q_e,sigma_e,&
!! Electrons linear coefficients for moment RHS !!!!!!!!!! parray_e(ips_e:ipe_e),jarray_e(ijs_e:ije_e),ips_e,ipe_e,ijs_e,ije_e,&
IF(KIN_E)THEN xnepj,xnepp1j,xnepm1j,xnepp2j,xnepm2j,xnepjp1,xnepjm1,&
DO ip = ips_e, ipe_e ynepp1j,ynepm1j,ynepp1jm1,ynepm1jm1,zNepm1j,zNepm1jp1,zNepm1jm1,&
p_int= parray_e(ip) ! Hermite degree xphij_e,xphijp1_e,xphijm1_e,xpsij_e,xpsijp1_e,xpsijm1_e)
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! All Napj terms
xnepj(ip,ij) = taue_qe*(CurvB*(2._dp*p_dp + 1._dp) &
+GradB*(2._dp*j_dp + 1._dp))
! Mirror force terms
ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp)
ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp)
ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp)
ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
zNepm1j (ip,ij) = +SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp)
zNepm1jp1(ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(p_dp)
zNepm1jm1(ip,ij) = -SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
ENDDO
ENDDO
DO ip = ips_e, ipe_e
p_int= parray_e(ip) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
! Landau damping coefficients (ddz napj term)
xnepp1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp + 1_dp)
xnepm1j(ip) = SQRT(tau_e)/sigma_e * SQRT(p_dp)
! Magnetic curvature term
xnepp2j(ip) = taue_qe * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
xnepm2j(ip) = taue_qe * CurvB * SQRT(p_dp * (p_dp - 1._dp))
ENDDO
DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! Magnetic gradient term
xnepjp1(ij) = -taue_qe * GradB * (j_dp + 1._dp)
xnepjm1(ij) = -taue_qe * GradB * j_dp
ENDDO
ENDIF ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ions linear coefficients for moment RHS !!!!!!!!!! CALL lin_coeff(k_Ti,k_Ni,k_cB,k_gB,tau_i,q_i,sigma_i,&
DO ip = ips_i, ipe_i parray_i(ips_i:ipe_i),jarray_i(ijs_i:ije_i),ips_i,ipe_i,ijs_i,ije_i,&
p_int= parray_i(ip) ! Hermite degree xnipj,xnipp1j,xnipm1j,xnipp2j,xnipm2j,xnipjp1,xnipjm1,&
p_dp = REAL(p_int,dp) ! REAL of Hermite degree ynipp1j,ynipm1j,ynipp1jm1,ynipm1jm1,zNipm1j,zNipm1jp1,zNipm1jm1,&
DO ij = ijs_i, ije_i xphij_i,xphijp1_i,xphijm1_i,xpsij_i,xpsijp1_i,xpsijm1_i)
j_int= jarray_i(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree CONTAINS
! All Napj terms SUBROUTINE lin_coeff(k_Ta,k_Na,k_cB,k_gB,tau_a,q_a,sigma_a,&
xnipj(ip,ij) = taui_qi*(CurvB*(2._dp*p_dp + 1._dp) & parray_a,jarray_a,ips_a,ipe_a,ijs_a,ije_a,&
+GradB*(2._dp*j_dp + 1._dp)) xnapj,xnapp1j,xnapm1j,xnapp2j,xnapm2j,xnapjp1,xnapjm1,&
! Mirror force terms ynapp1j,ynapm1j,ynapp1jm1,ynapm1jm1,zNapm1j,zNapm1jp1,zNapm1jm1,&
ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp+1._dp) xphij_a,xphijp1_a,xphijm1_a,xpsij_a,xpsijp1_a,xpsijm1_a)
ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1)*SQRT(p_dp) IMPLICIT NONE
ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp+1._dp) ! INPUTS
ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp) REAL(dp), INTENT(IN) :: k_Ta,k_Na,k_cB,k_gB,tau_a,q_a,sigma_a
! Trapping terms INTEGER, DIMENSION(ips_a:ipe_a), INTENT(IN) :: parray_a
zNipm1j (ip,ij) = +SQRT(tau_i)/sigma_i* (2._dp*j_dp+1_dp)*SQRT(p_dp) INTEGER, DIMENSION(ijs_a:ije_a), INTENT(IN) :: jarray_a
zNipm1jp1(ip,ij) = -SQRT(tau_i)/sigma_i* (j_dp+1_dp)*SQRT(p_dp) INTEGER, INTENT(IN) :: ips_a,ipe_a,ijs_a,ije_a
zNipm1jm1(ip,ij) = -SQRT(tau_i)/sigma_i* j_dp*SQRT(p_dp) ! OUTPUTS (linear coefficients used in moment_eq_rhs_mod.F90)
ENDDO REAL(dp), DIMENSION(ips_a:ipe_a,ijs_a:ije_a), INTENT(OUT) :: xnapj
ENDDO REAL(dp), DIMENSION(ips_a:ipe_a), INTENT(OUT) :: xnapp1j, xnapm1j, xnapp2j, xnapm2j
DO ip = ips_i, ipe_i REAL(dp), DIMENSION(ijs_a:ije_a), INTENT(OUT) :: xnapjp1, xnapjm1
p_int= parray_i(ip) ! Hermite degree REAL(dp), DIMENSION(ips_a:ipe_a,ijs_a:ije_a), INTENT(OUT) :: ynapp1j, ynapm1j, ynapp1jm1, ynapm1jm1
p_dp = REAL(p_int,dp) ! REAL of Hermite degree REAL(dp), DIMENSION(ips_a:ipe_a,ijs_a:ije_a), INTENT(OUT) :: zNapm1j, zNapm1jp1, zNapm1jm1
! Landau damping coefficients (ddz napj term) REAL(dp), DIMENSION(ips_a:ipe_a,ijs_a:ije_a), INTENT(OUT) :: xphij_a, xphijp1_a, xphijm1_a
xnipp1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp + 1._dp) REAL(dp), DIMENSION(ips_a:ipe_a,ijs_a:ije_a), INTENT(OUT) :: xpsij_a, xpsijp1_a, xpsijm1_a
xnipm1j(ip) = SQRT(tau_i)/sigma_i * SQRT(p_dp) INTEGER :: p_int, j_int ! polynom. dagrees
! Magnetic curvature term REAL(dp) :: p_dp, j_dp
xnipp2j(ip) = taui_qi * CurvB * SQRT((p_dp + 1._dp) * (p_dp + 2._dp)) !! linear coefficients for moment RHS !!!!!!!!!!
xnipm2j(ip) = taui_qi * CurvB * SQRT(p_dp * (p_dp - 1._dp)) DO ip = ips_a, ipe_a
ENDDO p_int= parray_a(ip) ! Hermite degree
DO ij = ijs_i, ije_i p_dp = REAL(p_int,dp) ! REAL of Hermite degree
j_int= jarray_i(ij) ! Laguerre degree DO ij = ijs_a, ije_a
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree j_int= jarray_a(ij) ! Laguerre degree
! Magnetic gradient term j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
xnipjp1(ij) = -taui_qi * GradB * (j_dp + 1._dp) ! All Napj terms
xnipjm1(ij) = -taui_qi * GradB * j_dp xnapj(ip,ij) = tau_a/q_a*(k_cB*(2._dp*p_dp + 1._dp) &
ENDDO +k_gB*(2._dp*j_dp + 1._dp))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Mirror force terms
!! ES linear coefficients for moment RHS !!!!!!!!!! ynapp1j (ip,ij) = -SQRT(tau_a)/sigma_a * (j_dp+1._dp)*SQRT(p_dp+1._dp)
IF (KIN_E) THEN ynapm1j (ip,ij) = -SQRT(tau_a)/sigma_a * (j_dp+1._dp)*SQRT(p_dp)
DO ip = ips_e, ipe_e ! ynapp1jm1(ip,ij) = +SQRT(tau_a)/sigma_a * j_dp*SQRT(p_dp+1._dp) ! Version of BJF
p_int= parray_e(ip) ! Hermite degree ! ynapm1jm1(ip,ij) = +SQRT(tau_a)/sigma_a * j_dp*SQRT(p_dp)
DO ij = ijs_e, ije_e ynapp1jm1(ip,ij) = +SQRT(tau_a)/sigma_a *(2._dp*j_dp-1._dp)*SQRT(p_dp+1._dp)
j_int= jarray_e(ij) ! REALof Laguerre degree ynapm1jm1(ip,ij) = +SQRT(tau_a)/sigma_a *(2._dp*j_dp-1._dp)*SQRT(p_dp)
j_dp = REAL(j_int,dp) ! REALof Laguerre degree ! Trapping terms
!! Electrostatic potential pj terms zNapm1j (ip,ij) = +SQRT(tau_a)/sigma_a *(2._dp*j_dp+1._dp)*SQRT(p_dp)
IF (p_int .EQ. 0) THEN ! kronecker p0 zNapm1jp1(ip,ij) = -SQRT(tau_a)/sigma_a * (j_dp+1._dp)*SQRT(p_dp)
xphij_e(ip,ij) = +k_Ne+ 2.*j_dp*k_Te zNapm1jm1(ip,ij) = -SQRT(tau_a)/sigma_a * j_dp*SQRT(p_dp)
xphijp1_e(ip,ij) = -k_Te*(j_dp+1._dp) ENDDO
xphijm1_e(ip,ij) = -k_Te* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
xphij_e(ip,ij) = +k_Te/SQRT2
xphijp1_e(ip,ij) = 0._dp; xphijm1_e(ip,ij) = 0._dp;
ELSE
xphij_e(ip,ij) = 0._dp; xphijp1_e(ip,ij) = 0._dp
xphijm1_e(ip,ij) = 0._dp;
ENDIF
ENDDO ENDDO
ENDDO DO ip = ips_a, ipe_a
ENDIF p_int= parray_a(ip) ! Hermite degree
DO ip = ips_i, ipe_i p_dp = REAL(p_int,dp) ! REAL of Hermite degree
p_int= parray_i(ip) ! Hermite degree ! Landau damping coefficients (ddz napj term)
DO ij = ijs_i, ije_i xnapp1j(ip) = SQRT(tau_a)/sigma_a * SQRT(p_dp+1._dp)
j_int= jarray_i(ij) ! REALof Laguerre degree xnapm1j(ip) = SQRT(tau_a)/sigma_a * SQRT(p_dp)
j_dp = REAL(j_int,dp) ! REALof Laguerre degree ! Magnetic curvature term
!! Electrostatic potential pj terms xnapp2j(ip) = tau_a/q_a * k_cB * SQRT((p_dp+1._dp)*(p_dp + 2._dp))
IF (p_int .EQ. 0) THEN ! kronecker p0 xnapm2j(ip) = tau_a/q_a * k_cB * SQRT( p_dp *(p_dp - 1._dp))
xphij_i(ip,ij) = +k_Ni + 2._dp*j_dp*k_Ti
xphijp1_i(ip,ij) = -k_Ti*(j_dp+1._dp)
xphijm1_i(ip,ij) = -k_Ti* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
xphij_i(ip,ij) = +k_Ti/SQRT2
xphijp1_i(ip,ij) = 0._dp; xphijm1_i(ip,ij) = 0._dp;
ELSE
xphij_i(ip,ij) = 0._dp; xphijp1_i(ip,ij) = 0._dp
xphijm1_i(ip,ij) = 0._dp;
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! EM linear coefficients for moment RHS !!!!!!!!!!
IF (KIN_E) THEN
DO ip = ips_e, ipe_e
p_int= parray_e(ip) ! Hermite degree
DO ij = ijs_e, ije_e
j_int= jarray_e(ij) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
xpsij_e (ip,ij) = +(k_Ne + (2._dp*j_dp+1._dp)*k_Te)* SQRT(tau_e)/sigma_e
xpsijp1_e(ip,ij) = - k_Te*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
xpsijm1_e(ip,ij) = - k_Te* j_dp * SQRT(tau_e)/sigma_e
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
xpsij_e (ip,ij) = + k_Te*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
xpsijp1_e(ip,ij) = 0._dp; xpsijm1_e(ip,ij) = 0._dp;
ELSE
xpsij_e (ip,ij) = 0._dp; xpsijp1_e(ip,ij) = 0._dp
xpsijm1_e(ip,ij) = 0._dp;
ENDIF
ENDDO ENDDO
ENDDO DO ij = ijs_a, ije_a
ENDIF j_int= jarray_a(ij) ! Laguerre degree
DO ip = ips_i, ipe_i j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
p_int= parray_i(ip) ! Hermite degree ! Magnetic gradient term
DO ij = ijs_i, ije_i xnapjp1(ij) = -tau_a/q_a * k_gB * (j_dp + 1._dp)
j_int= jarray_i(ij) ! REALof Laguerre degree xnapjm1(ij) = -tau_a/q_a * k_gB * j_dp
j_dp = REAL(j_int,dp) ! REALof Laguerre degree ENDDO
!! Electrostatic potential pj terms !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF (p_int .EQ. 1) THEN ! kronecker p1 !! ES linear coefficients for moment RHS !!!!!!!!!!
xpsij_i (ip,ij) = +(k_Ni + (2._dp*j_dp+1._dp)*k_Ti)* SQRT(tau_i)/sigma_i DO ip = ips_a, ipe_a
xpsijp1_i(ip,ij) = - k_Ti*(j_dp+1._dp) * SQRT(tau_i)/sigma_i p_int= parray_a(ip) ! Hermite degree
xpsijm1_i(ip,ij) = - k_Ti* j_dp * SQRT(tau_i)/sigma_i DO ij = ijs_a, ije_a
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3 j_int= jarray_a(ij) ! REALof Laguerre degree
xpsij_i (ip,ij) = + k_Ti*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i j_dp = REAL(j_int,dp) ! REALof Laguerre degree
xpsijp1_i(ip,ij) = 0._dp; xpsijm1_i(ip,ij) = 0._dp; !! Electrostatic potential pj terms
ELSE IF (p_int .EQ. 0) THEN ! kronecker p0
xpsij_i (ip,ij) = 0._dp; xpsijp1_i(ip,ij) = 0._dp xphij_a(ip,ij) = +k_Na + 2._dp*j_dp*k_Ta
xpsijm1_i(ip,ij) = 0._dp; xphijp1_a(ip,ij) = -k_Ta*(j_dp+1._dp)
ENDIF xphijm1_a(ip,ij) = -k_Ta* j_dp
ENDDO ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
ENDDO xphij_a(ip,ij) = +k_Ta/SQRT2
xphijp1_a(ip,ij) = 0._dp; xphijm1_a(ip,ij) = 0._dp;
ELSE
xphij_a(ip,ij) = 0._dp; xphijp1_a(ip,ij) = 0._dp
xphijm1_a(ip,ij) = 0._dp;
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Electromagnatic linear coefficients for moment RHS !!!!!!!!!!
DO ip = ips_a, ipe_a
p_int= parray_a(ip) ! Hermite degree
DO ij = ijs_a, ije_a
j_int= jarray_a(ij) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
IF (p_int .EQ. 1) THEN ! kronecker p1
xpsij_a (ip,ij) = +(k_Na + (2._dp*j_dp+1._dp)*k_Ta)* SQRT(tau_a)/sigma_a
xpsijp1_a(ip,ij) = - k_Ta*(j_dp+1._dp) * SQRT(tau_a)/sigma_a
xpsijm1_a(ip,ij) = - k_Ta* j_dp * SQRT(tau_a)/sigma_a
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
xpsij_a (ip,ij) = + k_Ta*SQRT3/SQRT2 * SQRT(tau_a)/sigma_a
xpsijp1_a(ip,ij) = 0._dp; xpsijm1_a(ip,ij) = 0._dp;
ELSE
xpsij_a (ip,ij) = 0._dp; xpsijp1_a(ip,ij) = 0._dp
xpsijm1_a(ip,ij) = 0._dp;
ENDIF
ENDDO
ENDDO
END SUBROUTINE lin_coeff
END SUBROUTINE compute_lin_coeff END SUBROUTINE compute_lin_coeff
END MODULE numerics END MODULE numerics
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