Homec4science
Diffusion COSOlver (new_coulomb) Browse

cosolver/src/improved_sugama_mod.f90.bkp
83b8f38f6c51new_coulomb

Last Change
Raw File

improved_sugama_mod.f90.bkp
View Options

MODULE improved_Sugama
!
! Linearized DK/GK Improved Sugamam Collision operator between species a and b
!
! Note: as defined below, the Sugama operator differs from a 2 factor from the original paper due to the definition of the isothermal part of the linearized Coulomb
!
! Ref: Sugama et. al., Phsyics of Plasmas 26 (2019)
USE prec_const
USE basic
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
USE auxval
use coulomb_parts
use sugama
IMPLICIT NONE
PUBLIC
CONTAINS
!_________________________________________________________________
!
! Drift kientic Imoroved Sugama collision operator
SUBROUTINE SugamaImprovedTDK(l,k,rCabT)
! Add the test correction matrix to DK Sugama. The improved Sugama collision operator reproduces the same friction-flow relations as the lineared DK full Coulomb collision operator [see Eqs. (12) and (13) in Ref. [1]]
!
!
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabT
!
!
! local vars.
INTEGER :: j,i,g,h,icol
TYPE(FM), SAVE :: NX0,NX1,T4,T4inv,NX2,cj_,ci_
TYPE(FM), SAVE :: DeltaMablk_
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! construct DK Test Sugama
CALL SugamaTDK(l,k, rCabT)
!
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
jloop: DO j=0, k + floor(real(l)/2)
NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
cj_ = CJ(j)
T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
IF( T4inv .ne. TO_FM(0)) THEN
iloop: DO i=0,imaxx ! .... energy moment
ci_ = CJ(i)
! Complute test correction matrix
DeltaMablk_= DELTAMlk(j,i)
gloop: DO g =0,min(2*i +1,Pmaxa)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0,min(i,Jmaxa)
T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
T4 = TO_FM(T4str_)
icol = numba(g,h)
rCabT(icol) = rCabT(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*DeltaMablk_
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedTDK
!
!_________________________________________________________________
!
SUBROUTINE SugamaImprovedFDK(l,k,rCabF)
! Add field correction matrix to DK Sugama. The improved Sugama collision operator reproduces the same friction-flow relations as the lineared DK full Coulomb collision operator [see Eqs. (12) and (13) in Ref. [1]]
!
!
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabF
!
!
! local vars.
INTEGER :: j,i,g,h,icol
TYPE(FM), SAVE :: NX0,NX1,T4,T4inv,NX2,cj_,ci_
TYPE(FM), SAVE :: DeltaNablk_
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! construct Sugama Field part
CALL SugamaFDK(l,k, rCabF)
!
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
jloop: DO j=0, k + floor(real(l)/2)
NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
cj_ = CJ(j)
T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
IF( T4inv .ne. TO_FM(0)) THEN
iloop: DO i=0,imaxx ! ... energy moment
ci_ = CJ(i)
! Complute test correction matrix
DeltaNablk_= DELTANlk(j,i)
gloop: DO g =0,min(2*i +1,Pmaxb)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0,min(i,Jmaxb)
T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
T4 = TO_FM(T4str_)
icol = numbb(g,h)
rCabF(icol) = rCabF(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*DeltaNablk_
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedFDK
!_________________________________________________________________
!
!
SUBROUTINE SugamaImprovedSelfDK(l,k,rCaa)
! Add correction matrix to self DK Sugama. The improved Sugama collision operator reproduces the same friction-flow relations as the lineared DK full Coulomb collision operator [see Eqs. (12) and (13) in Ref. [1]]
!
!
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCaa
!
!
! local vars.
INTEGER :: j,i,g,h,icol
TYPE(FM), SAVE :: NX0,NX1,T4,T4inv,NX2,cj_,ci_
TYPE(FM), SAVE :: DeltaMaalk_,DeltaNaalk_
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! Construct self DK Sugama
CALL SugamaSelfDK(l,k,rCaa)
!
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
jloop: DO j=0, k + floor(real(l)/2)
NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
cj_ = CJ(j)
T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
IF( T4inv .ne. TO_FM(0)) THEN
iloop: DO i=0,imaxx
ci_ = CJ(i)
! Complute test and field correction matrix
DeltaNaalk_= DeltaNlk(j,i)
DeltaMaalk_= DeltaMlk(j,i)
gloop: DO g =0,min(2*i +1,Pmaxa)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0,min(i,Jmaxa)
T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
T4 = TO_FM(T4str_)
icol = numba(g,h)
rCaa(icol) = rCaa(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*(DeltaNaalk_ + DeltaMaalk_)
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedSelfDK
!
!
! GYROKINETIC IMPROVED FIELD SUGAMA
!
!
!_________________________________________________________________
SUBROUTINE SugamaImprovedTGK(l,k,rCabT)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCabT
!
CALL FM_ENTER_USER_ROUTINE
!
! CODE HERE ...... !!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedTGK
!
!__________________________________________
!
SUBROUTINE SugamaImprovedFGK(l,k,rCabF)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCabF
!
CALL FM_ENTER_USER_ROUTINE
!
! CODE HERE ...... !!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedFGK
END MODULE improved_Sugama
c4science · Help