Homec4science
Diffusion COSOlver (new_coulomb) Browse

cosolver/src/coulomb_parts_mod.f90
41e7a23b16aenew_coulomb

Last Change
Raw File

coulomb_parts_mod.f90
View Options

module coulomb_parts
! contains non isothermal and FLR corrections of the test component of the Coulomb collision operator
use prec_const
use FMZM
use basic
use model
use basis_transformation
use array
use coeff
use speed_functions
implicit none
contains
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)
IF(only_symmetric) icolend = numba(l,k)
!
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)
!______________________________________________
! INVERSE BASIS TRANSFORMATION
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
! NEW VERSION
! NX3 = TO_FM(curlyAshrg(s,h,r,g))
! rCabT(icol) = rCabT(icol) + NX*NX0*NX1*T4inv1*T4inv2*NX3
ENDIF
ENDDO sloop
ENDIF
ENDDO rloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDDO colloop
!
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Rablk
!
SUBROUTINE compute_curlyAshrg
! prec compute
!
! \sum_{s1=0}^s \sum_{h1=0}^h L_{hh1}^g L_{ss1}^r \alpha_\perp^{grs1h1}
implicit none
!
! locs. var
integer :: s1,h1,s,h,r,g
integer :: smaxx,hmaxx,rmaxx,gmaxx
TYPE(FM), save :: buff_,Lrss1,Lghh1,alpha0_
!
CALL FM_ENTER_USER_ROUTINE
!
smaxx = Jmaxa + floor(real(Pmaxa,dp)/2._dp)
hmaxx = Jmaxa + floor(real(Pmaxa,dp)/2._dp)
!
rmaxx = Pmaxa + 2*Jmaxa
gmaxx = Pmaxa + 2*Jmaxa
! allocate array for prec-computed coefficients for spatial diffusion in test part
call allocate_array(curlyAshrg, 0, smaxx, &
0, hmaxx, &
0, rmaxx, &
0, gmaxx )
DO s=0,smaxx
DO h=0,hmaxx
DO r=0,rmaxx
DO g=0,gmaxx
! compute coeff
buff_ = ZEROFM
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))
!
alpha0_ = ZEROFM
IF( r .eq. g) THEN
alpha0_ = (TO_FM(1)/TO_FM(2*g+1) )*( nuabparak(g + 2 + s1 +h1) + nuabDk( g+2 +s1 +h1 ))
ENDIF
IF( r .eq. g .or. r .eq. g-2 .or. r .eq. g+2) THEN
alpha0_ = alpha0_ + X2prodLegendre(r,g)/TO_FM(2)*(nuabDk( (r +g)/2 +2 +s1 +h1 ) - nuabparak((g+r)/2 + 2 + s1 +h1))
ENDIF
! coefficient
buff_ = buff_ + alpha0_*Lrss1*Lghh1
!
ENDDO h1loop
ENDDO s1loop
! store it
curlyAshrg(s,h,r,g) = to_dp(buff_)
ENDDO
ENDDO
ENDDO
ENDDO
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE compute_curlyAshrg
!
END MODULE coulomb_parts
c4science · Help