Homec4science
Diffusion COSOlver (improved_sugama) Browse

cosolver/src/coll_lorentz.h
cd1286e2d410improved_sugama

Last Change
Raw File

coll_lorentz.h
View Options

SUBROUTINE LorentzDK(l,k,rCabT)
! Electron Drift-kinetic Lorentz Pitch-angle scattering Operator
!
! Note:
! * Remove the element p=0 from the summation since Injlkmp is not defined for p=0
USE prec_const
IMPLICIT NONE
!
!
INTEGER, INTENT(in) :: l,k
! REAL(dp), DIMENSION(:,:),INTENT(inout) :: CabT
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
!
INTEGER :: p,j,g,h,t,t1,j1
TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T4,SPJ,Lptt1,Lpjj1,T4inv
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_,T4invstr_
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
! Pre compute quantities
!
! Hermite normalization coeff
NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*(TO_FM(2)**l))
!
jloop:DO j=0,jamaxx
ploop:DO p=1,l+2*k
NX1 = 2*p*(p+1)*TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p+1))/SQRT(PIFM)
SPJ = sigmapj(p,j)
tloop: DO t=0, k + l/2
!______________________________________________
! INVERSE BASIS TRANSFORMATION
T4str_ = GET_T4(p,t,l,k) ! ...
T4inv = 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)*TO_FM(T4str_)
!______________________________________________
IF (T4INV .ne. TO_FM('0.0')) THEN
gloop:DO g=0, min(Pmaxa,p+2*j)
NX2 = SQRT(FACTORIAL(TO_FM(g))*TO_FM(2)**g)
hloop:DO h=0, min(Jmaxa,j +p/2)
!_______________________________________
! BASIS TRANSFORMATION
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF (T4 .ne. TO_FM('0.0')) THEN
!_______________________________________
t1loop: DO t1=0, t
Lptt1=Lpjl(real(p,dp),real(t,dp),real(t1,dp))
j1loop:DO j1=0, j
Lpjj1=Lpjl(real(p,dp),real(j,dp),real(j1,dp))
NX3 = FACTORIAL(TO_FM( p + t1 + j1 -1))
! CONSTRUCT Lorentz COLLISIONAL MATRIX
icol = numba(g,h)
rCabT(icol ) = rCabT(icol ) - NX1*NX0*NX2*NX3*T4*T4inv*Lptt1*Lpjj1/SPJ
ENDDO j1loop
ENDDO t1loop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE LorentzDK
!===============================================================================
SUBROUTINE LorentzGK(l,k,rCabT)
! Electron GYROKINETIC Lorentz Pitch-angle scattering Operator
!
! Note:
! * Remove the element p=0 from the summation since Injlkmp is not defined for p=0
USE prec_const
IMPLICIT NONE
!
!
INTEGER, INTENT(in) :: l,k
! REAL(dp), DIMENSION(:,:),INTENT(inout) :: CabT
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
!
INTEGER :: p,j,n,g,h,s,np,f,t1,j1,t,s1
TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T4,d0nhs_,SPJ,curlyA_,bdmnkf,T4inv,Lptt1,Lpjj1,dnphs1,ker_
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_
CALL FM_ENTER_USER_ROUTINE
!
!
! Pre compute quantities
!
! Hermite normalization coeff
NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*(TO_FM(2)**l))
!
jloop:DO j=0,jamaxx
nloop: DO n=0,naFLR
curlyA_ = curlyAmn(0,n,bafm_)
floop: DO f=0,n +k
bdmnkf= ALL2L(n,k,f,0) ! ... overline{d}_{nkf}^0
ploop:DO p=1,l + 2*f
NX1 = 2*p*(p+1)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p+1))/SQRT(PIFM)
SPJ = sigmapj(p,j)
tloop: DO t=0, f + l/2
!______________________________________________
! INVERSE BASIS TRANSFORMATION
T4str_ = GET_T4(p,t,l,f) ! ...
T4inv = 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)*TO_FM(T4str_)
!______________________________________________
IF (T4INV .ne. TO_FM('0.0')) THEN
gloop:DO g=0, min(Pmaxa,p+2*j)
NX2 = SQRT(FACTORIAL(TO_FM(g))*TO_FM(2)**g)
hloop:DO h=0, j +p/2
!_______________________________________
! BASIS TRANSFORMATION
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF (T4 .ne. TO_FM('0.0')) THEN
!_______________________________________
nploop: DO np=0, naFLR
ker_ = kerneln(bafm_,np)
s1loop: DO s1=0,min(Jmaxa,np+h)
dnphs1 = ALLX2L(np,h,s1,0)
t1loop: DO t1=0, t
Lptt1=Lpjl(real(p,dp),real(t,dp),real(t1,dp))
j1loop:DO j1=0, j
Lpjj1=Lpjl(real(p,dp),real(j,dp),real(j1,dp))
NX3 = FACTORIAL(TO_FM( p + t1 + j1 -1))
! CONSTRUCT Lorentz COLLISIONAL MATRIX
icol = numba(g,s1)
rCabT(icol ) = rCabT(icol ) - NX1*NX0*NX2*NX3*T4*T4inv*ker_*Lptt1*curlyA_*bdmnkf*dnphs1*Lpjj1/SPJ
ENDDO j1loop
ENDDO t1loop
ENDDO s1loop
ENDDO nploop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO floop
ENDDO nloop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE LorentzGK
!============================================================================
SUBROUTINE Cei1FDK(l,k,rCabF)
! Drift-kinetic Leading order sqrt(m_e/m_i) momentum convservative field part of Lorentz pitch-angle operator
!
!
USE prec_const
IMPLICIT NONE
!
!
INTEGER,INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabF
!
TYPE(FM),SAVE :: NX0
INTEGER ::irow
CALL FM_ENTER_USER_ROUTINE
!
!
NX0 = curlyAlk0pareiDK(l,k)
!
rCabF(numbb(1,0)) = rCabF(numbb(1,0)) +8*NX0/3/SQRT(TO_FM(2))/chifm_/PIFM
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Cei1FDK
!=================================================================================
SUBROUTINE Cei1FGK(l,k,rCabF)
! Gyrokinetic Leading order sqrt(m_e/m_i) momentum convservative field part of Lorentz pitch-angle operator
!
!
USE prec_const
IMPLICIT NONE
!
!
INTEGER,INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabF
!
TYPE(FM),SAVE :: NX0,NX1,ker_
INTEGER ::irow,n
!
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
NX0 = curlyAlk0pareiGK(l,k)
NX1 = curlyAlk0perpeiGK(l,k)
nloop: DO n=0,niFLR
ker_ = kerneln(bbfm_,n)
rCabF(numbb(1,n)) = rCabF(numbb(1,n)) +ker_*8*NX0/3/SQRT(TO_FM(2))/chifm_/PIFM
rCabF(numbb(0,n)) = rCabF(numbb(0,n)) - ker_*4/SQRT(PIFM)/chifm_*bbfm_*NX1
rCabF(numbb(0,n+1)) = rCabF(numbb(0,n+1)) + ker_*4/SQRT(PIFM)/chifm_*bbfm_*NX1
!
!
ENDDO nloop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Cei1FGK
!================================================================================
SUBROUTINE DeflectionTDK(l,k,rCabT)
!
! Construct the Deflection pitch-angle scattering operator with velocity dependent deflection frequency \nu_{ab}^D(v)
!
!
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,nusD
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,jamaxx
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))
nusD = nusabDpjd(p,j,d)
rCabT(icol) = rCabT(icol) - TO_FM(p)*(TO_FM(p) +1)*NX0*NX1*T4inv*NX2*T4*Lphd*nusD/SPJ
ENDDO dloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
CALL FM_EXIT_USER_ROUTINE
!
!
END SUBROUTINE DeflectionTDK
!================================================================================
c4science · Help