Homec4science
Diffusion COSOlver (llechot) Browse

cosolver/src/coulomb_mod.f90
df55999e752fllechot

Last Change
Raw File

coulomb_mod.f90
View Options

MODULE coulomb
!
! Linearized Full DK Coulomb Collision operator between species a and b
!
USE prec_const
USE basic
USE basic_mpi
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
IMPLICIT NONE
PUBLIC
CONTAINS
!===============================================================================
SUBROUTINE CoulombTDK(l,k,rCabT)
!
! Construct the test drift-kinetic Coulomb operator
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
! Declare local variables
INTEGER :: p,j,n,g,h,s1,np,t,f,d
TYPE(FM),SAVE :: NX0,NX1,NX2,T4,SPJ,bdmnkf,nusT,T4inv,Lptd
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_, T4invstr_
!
!
CALL FM_ENTER_USER_ROUTINE
!
!
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
irow = numba(l,k)
!
jloop: DO j=0,JEmaxx
ploop: DO p=0,l+2*k ! ... min(pamaxx,l+2*k)
SPJ = sigmapj(p,j)
NX2 = (TO_FM(2)**p)*(FACTORIAL(TO_FM(p)))*(FACTORIAL(TO_FM(p)))/FACTORIAL(TO_FM(2*p))/TO_FM(2*p+1)
tloop: DO t =0, k+ floor(real(l/2))
T4invstr_ = GET_T4(p,t,l,k) ! ... (T^-1)^pt_lk = ... T^lk_pt
T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*TO_FM(2)**l*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM(2*p + 1)/2/FACTORIAL(TO_FM(2*p+2*t +1)/2)
IF( T4inv .ne. ZEROFM) THEN
hloop: DO h=0, min(j+floor(real(p)/2),Jmaxa)
gloop: DO g=0,min(p+2*j,Pmaxa)
NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF( T4 .ne. TO_FM('0.0')) THEN
dloop: DO d=0,t
nusT = nusTabpjgd(p,j,p,d)
Lptd=Lpjl(real(p,dp),real(t,dp),real(d,dp))
icol = numba(g,h)
rCabT(icol) = rCabT(icol ) + nusT*Lptd*T4*T4inv*NX0*NX1*NX2/SPJ
ENDDO dloop
ENDIF
ENDDO gloop
ENDDO hloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
END SUBROUTINE CoulombTDK
! ===============================================================================
SUBROUTINE CoulombFDK(l,k,rCabF)
!
! Construct the field (back-reaction) DRIFT-KINETIC Coulomb operator
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabF
! Declare local variables
INTEGER :: p,j,n,g,h,s1,np,t,f,d
TYPE(FM),SAVE :: NX0,NX1,NX2,T4,SPJ,nusF,T4inv,Lptd
INTEGER :: icol
CHARACTER(len=T4len) :: T4str_, T4invstr_
!
!
CALL FM_ENTER_USER_ROUTINE
!
!
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
jloop: DO j=0,JEmaxx
ploop: DO p=0,l+2*k ! ... min(pamaxx,l+2*k)
SPJ = sigmapj(p,j)
NX2 = TO_FM(2)**p*(FACTORIAL(TO_FM(p)))*(FACTORIAL(TO_FM(p)))/FACTORIAL(TO_FM(2*p))/TO_FM(2*p+1)
tloop: DO t =0, k+ floor(real(l)/2)
T4invstr_ = GET_T4(p,t,l,k) ! ... (T^-1)^pt_lk = ... T^lk_pt
T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*TO_FM(2)**l*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM(2*p + 1)/2/FACTORIAL(TO_FM(2*p+2*t +1)/2)
IF( T4inv .ne. TO_FM('0.0')) THEN
gloop: DO g=0,min(p+2*j,Pmaxb)
NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0, min(j+floor(real(p)/2),Jmaxb)
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF( T4 .ne. TO_FM('0.0')) THEN
dloop: DO d=0,t
nusF = nusFabpjgd(p,j,p,d)
Lptd=Lpjl(real(p,dp),real(t,dp),real(d,dp))
icol = numbb(g,h)
rCabF(icol) = rCabF(icol ) + nusF*Lptd*T4*T4inv*NX0*NX1*NX2/SPJ
ENDDO dloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
END SUBROUTINE CoulombFDK
!
!_______________________________________________________________
!
SUBROUTINE CoulombSelfDK(l,k,rCaa)
!
! Construct the drift-kinetic like-species Coulomb collisional matrix
!
! NOTE : Tested and benchmarked against Mathematica code.
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCaa
! Declare local variables
INTEGER :: p,j,n,g,h,s1,np,t,f,d
TYPE(FM),SAVE :: NX0,NX1,NX2,T4,SPJ,nusT,nusF,T4inv,Lptd
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4invstr_,T4str_
!
!
CALL FM_ENTER_USER_ROUTINE
!
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
irow = numba(l,k)
!
jloop: DO j=0,JEmaxx
ploop: DO p=0,l+2*k ! ... min(pamaxx,l+2*k)
SPJ = sigmapj(p,j)
NX2 = TO_FM(2)**p*(FACTORIAL(TO_FM(p)))*(FACTORIAL(TO_FM(p)))/FACTORIAL(TO_FM(2*p))/TO_FM(2*p+1)
tloop: DO t =0, k+ floor(real(l)/2)
T4invstr_ = GET_T4(p,t,l,k) ! ... (T^-1)^pt_lk = ... T^lk_pt
T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*TO_FM(2)**l*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM(2*p + 1)/2/FACTORIAL(TO_FM(2*p+2*t +1)/2)
IF( T4inv .ne. ZEROFM) THEN
gloop: DO g=0,min(p+2*j,Pmaxa)
NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0, min(j+floor(real(p)/2),Jmaxa)
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF( T4 .ne. ZEROFM) THEN
dloop: DO d=0,t
nusT = ZEROFM
IF(ADDTEST) THEN
nusT= nusTabpjgd(p,j,p,d)
ENDIF
nusF = ZEROFM
IF(ADDBACK) THEN
nusF = nusFabpjgd(p,j,p,d)
ENDIF
Lptd=Lpjl(real(p,dp),real(t,dp),real(d,dp))
icol = numba(g,h)
rCaa(icol) = rCaa(icol ) + (nusT + nusF)*Lptd*T4*T4inv*NX0*NX1*NX2/SPJ
ENDDO dloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE CoulombSelfDK
!
!------------------------------------------------
!
SUBROUTINE Coulomb_isothermalTDK(l,k,rCabT)
!
! Construct the isothermal part of the test part of the DK linearized Coulomb. Used in sugama_mod and abel_mod
!
! Ref: Sugama et al., Physics of Plasmas 16, 112503 (2009)
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
! Declare local variables
INTEGER :: p,h,j,d,t,g
TYPE(FM),SAVE :: NX0,NX1,NX2,SPJ,Lphd,T4,T4inv,nusS
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_, T4invstr_
!
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
NX0 = 1/SQRT((TO_FM(2)**l)*FACTORIAL(TO_FM(l)))
!
jloop: DO j=0,JEmaxx
ploop: DO p=0,l + 2*k
SPJ = sigmapj(p,j)
NX1 = (TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM(2*p +1)
tloop: DO t=0, k + floor(real(l)/2)
T4invstr_ = GET_T4(p,t,l,k) ! ... (T^-1)^pt_lk = ... T^lk_pt
T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM(2*p + 1)/2/FACTORIAL(TO_FM(2*p+2*t +1)/2)
IF( T4inv .ne. TO_FM('0.0')) THEN
gloop: DO g=0, min(p + 2*j,Pmaxa)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0, min(j + floor(real(p)/2),Jmaxa)
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
icol = numba(g,h)
IF( T4 .ne. TO_FM('0.0')) THEN
dloop: DO d=0, t
Lphd=Lpjl(real(p,dp),real(t,dp),real(d,dp))
nusS = nusSabpjd(p,j,d)
rCabT(icol) = rCabT(icol) + NX0*NX1*T4inv*NX2*T4*Lphd*nusS/SPJ
ENDDO dloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
CALL FM_EXIT_USER_ROUTINE
!
!
END SUBROUTINE Coulomb_isothermalTDK
!
!-------------------------------------------
!
SUBROUTINE Rablk(l,k,rCabT)
!
! Computes the FLR corrections to isothermal part of test DK Coulomb
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! Local vars.
!
INTEGER :: p, j,g,h,r,s,s1,h1
INTEGER :: pj(2)
INTEGER :: icolend
TYPE(FM),SAVE :: Lrss1,Lghh1,T4inv1,T4inv2,NX0,NX1,NX3
TYPE(FM), SAVE :: NX
CHARACTER(len=T4len) :: T4str_
INTEGER :: icol
!
!
CALL FM_ENTER_USER_ROUTINE
!
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
NX = -TO_FM(1)/SQRT(PIFM)*bafm_*bafm_
!
!
icolend = numba(Pmaxa,Jmaxa)
!
colloop: do icol = 1,icolend
pj = invnumba(icol)
p = pj(1)
j = pj(2)
NX1 = TO_FM(1)/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
gloop: DO g=0, l + 2*k
hloop: DO h=0, k + floor(real(l)/2._dp)
T4str_ = GET_T4(g,h,l,k) ! ... T^{lk}_{gh}
T4inv1 = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM(2*g + 1)/2/FACTORIAL(TO_FM(2*g+2*h +1)/2)*TO_FM(T4str_)
IF(T4inv1 .ne. ZEROFM) THEN
rloop: DO r=0, p + 2*j
IF( (r .eq. g) .or. (r .eq. g-2 ) .or. (r .eq. g+2)) THEN
sloop: DO s=0, j + floor(real(p)/2._dp)
T4str_ = GET_T4(r,s,p,j) ! ... T^{pj}_{rs}
T4inv2 = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(s))*TO_FM(2*r + 1)/2/FACTORIAL(TO_FM(2*r+2*s +1)/2)*TO_FM(T4str_)
IF(T4inv2 .ne. ZEROFM) THEN
! FIRST VERSION
s1loop: DO s1=0, s
Lrss1 = Lpjl(real(r,dp),real(s,dp),real(s1,dp))
h1loop: DO h1=0,h
Lghh1 = Lpjl(real(g,dp),real(h,dp),real(h1,dp))
! compute speed functions
NX3 = TO_FM('0')
IF( r .eq. g) THEN
NX3 = (TO_FM(1)/TO_FM(2*g+1) )*( nuabparak(g + 2 + s1 +h1) + nuabDk( g+2 +s1 +h1 ))
ENDIF
NX3 = NX3 + X2prodLegendre(r,g)/TO_FM(2)*(nuabDk( (r +g)/2 +2 +s1 +h1 ) - nuabparak((g+r)/2 + 2 + s1 +h1))
!
rCabT(icol) = rCabT(icol) + NX*NX0*NX1*T4inv1*T4inv2*Lrss1*Lghh1*NX3
ENDDO h1loop
ENDDO s1loop
ENDIF
ENDDO sloop
ENDIF
ENDDO rloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDDO colloop
!
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Rablk
!
END MODULE coulomb
c4science · Help