Homec4science
Diffusion COSOlver (improved_sugama) Browse

cosolver/src/model_mod.f90
cd1286e2d410improved_sugama

Last Change
Raw File

model_mod.f90
View Options

module model
! phyical parameters
use prec_const
use FMZM
use basic
IMPLICIT NONE
PRIVATE
! physical parameters
real(dp), public, protected :: nu = 1._dp ! electron-ion collision freq.
real(dp), public, protected :: mu =1._dp ! sqr(m_e/m_i)
real(dp), public, protected :: tau = 1._dp ! T_i0/T_e0
real(dp), public, protected :: kperp = 1._dp ! perpendicular wavenumber normalized to cs
! electron-ion mass/temperature ratio
real(dp), public, protected :: sigmaei, sigmaaa,sigmaie
real(dp), public, protected :: tauei, tauie, tauaa
real(dp), public, protected :: sigma_ = 1._dp !
real(dp), public, protected :: tau_ = 1._dp !
real(dp), public, protected :: chi_ = 1._dp !
real(dp), public, protected :: be_ = 0._dp ! perpendicular wavenumber normalized to electron thermal Larmor radius
real(dp), public, protected :: bi_ = 0._dp ! perpendicular wavenumber normalized to electron thermal Larmor radius
real(dp), public, protected :: ba_ = 0._dp
real(dp), public, protected :: bb_ = 0._dp
! Multiple-Precision FM type
TYPE(FM), public, protected :: sigmafm_ ! mass ratio
TYPE(FM), public, protected :: taufm_ ! temperature ratio
TYPE(FM), public, protected :: chifm_ ! SQRT(tau/sigma) = vTa/vTb
TYPE(FM), public, protected :: thetaabfm_ ! Sugama theta_{ab} coefficient [see Eq. (29) in Sugama et al. (2009)]
TYPE(FM), public, protected :: thetabafm_ ! Sugma theta_{ba} coefficient
TYPE(FM), public, protected :: bafm_
TYPE(FM), public, protected :: bbfm_ !
! Poisson
TYPE(FM), public :: poisson_coeff ! ion polarization and adiabatic electron coeff.
! label *.h5 file
character(len=15), public :: strnu
character(len=15), public:: strmu
character(len=15), public :: strkperp
character(len=15), public :: strtau
character(len=1), PUBLIC :: strETEST, strITEST,strEBACK,strIBACK,strESELF,strISELF,strGKE,strGKI
! Test/Back-Reaction Model Collision Operator: 0 => Coulomb
INTEGER,PUBLIC :: ETEST=0 ! Electron Test Model
INTEGER,PUBLIC :: ITEST=0 ! Ion Test Model
INTEGER,PUBLIC :: EBACK=0 ! Electron Back-Reaction Model
INTEGER,PUBLIC :: IBACK=0 ! Ion Back-Reaction Model
INTEGER, PUBLIC :: ESELF=0 ! Electron self Model
INTEGER, PUBLIC :: ISELF=0 ! Ion Self Model
INTEGER, PUBLIC :: GKE =0 ! Gyrokinetic Electrons
INTEGER, PUBLIC :: GKI =0 ! Gyrokinetic Ions
! Colliding species
CHARACTER(len=2), PUBLIC, PROTECTED :: coll_ab
! Sugama Flag: used in AUXVAL for error functions necessary for back-reaction parts
LOGICAL, PUBLIC, PROTECTED :: IFSUGAMA = .FALSE.
! Improved Sugama Flag: used in AUXVAL and MEMORY for precomputation of CurlI(para/perp)(a/b)pjk in GK mode
LOGICAL, PUBLIC, PROTECTED :: IFIMPROVEDSUGAMA = .FALSE.
! Coulomb flag
LOGICAL, PUBLIC, PROTECTED :: IFCOULOMB = .FALSE.
! Polarization flag
LOGICAL, PUBLIC, PROTECTED :: IFPOL = .FALSE.
! Zeroth order Hirshman Sigmar Clarke flag
logical, PUBLIC, protected :: IF0HSC = .false.
! GK Abel flag: used in AUXVAL for error functions necessary for back-reaction parts
logical, PUBLIC, protected :: IFABEL = .false.
! DK/GK test or/and field: used only if GKI and/or GKE
! Note: implented only for Abel collision operator
logical, PUBLIC, protected :: DKBACK = .false. ! use DK back reaction
logical, PUBLIC, protected :: DKTEST = .false. ! use DK test part
! Flag turn on/off Test/Field parts
logical, public, protected :: ADDBACK =.true. ! IF 1 add field part, 0 otherwise
logical, public, protected :: ADDTEST=.true. ! IF 1 add test part, 0 otherwise
! Flag turn on/off spatial diff in test part (implemetend of GK Abel)
logical, public, protected :: only_gk_part = .false.
! Gyrokinetic flag
LOGICAL, PUBLIC, PROTECTED :: IFGK = .FALSE.
public :: model_readinputs, set_model,set_model_fm,operatormodel_readinput,invert_model
contains
!================================================================================
subroutine model_readinputs
! Read inputs physical input parameters
use basic, only: lu_in,ierr
use prec_const
implicit none
character(len=8) :: fmt = '(f6.4)'
NAMELIST /MODEL_PAR/ nu,mu,tau,kperp
IF( me .ne. 0) GOTO 111
read(lu_in,MODEL_PAR)
write(*,MODEL_PAR)
111 CONTINUE
! We compute the normalized collisional matrix
nu = 1._dp
!
! Broadcast from PE 0 to
!
CALL MPI_bcast(nu,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(mu,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(tau,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(kperp,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
! convert to string
write(strnu,fmt) nu
write(strmu,fmt) mu
write(strtau,fmt) tau
write(strkperp,fmt) kperp
! set electron-ion mass/temperature ratio
sigmaei = mu ! defined as me/mi !!!!!!
sigmaie = 1/sigmaei
sigmaaa = 1._dp
tauei = 1/tau
tauie = tau
tauaa = 1._dp
! Set electron/ion perpendicular wavenumber
be_ = kperp*SQRT(2._dp)*SQRT(sigmaei)
bi_ = kperp*SQRT(2._dp)*SQRT(tauei)
end subroutine model_readinputs
!================================================================================
subroutine model_outputinputs
! output physical input parameters to *.h5 file
use basic, only: lu_in,me
use prec_const
implicit none
! code here
end subroutine model_outputinputs
!================================================================================
subroutine set_model(ab)
! set the species mass/temperature ratios
implicit none
character(len=2), intent(in) :: ab
! Ion-Ion collisions
if(ab .eq. 'ii') then
chi_ = 1._dp
tau_ = 1._dp
sigma_ = 1._dp
ba_ = bi_
bb_ = bi_
coll_ab = 'ii'
end if
! Electron- Electron collisions
if(ab .eq. 'ee') then
chi_ = 1._dp
tau_ = 1._dp
sigma_ = 1._dp
ba_ = be_
bb_ = be_
coll_ab = 'ee'
end if
! Electron-ion collisions
if(ab .eq. 'ei') then
chi_ = SQRT(tauei/sigmaei)
tau_ = tauei
sigma_ = sigmaei
ba_ = be_
bb_ = bi_
coll_ab = 'ei'
end if
! Ion-electron collisions
if(ab .eq. 'ie') then
chi_ = SQRT(tauie/sigmaie)
tau_ = tauie
sigma_ = sigmaie
ba_ = bi_
bb_ = be_
coll_ab = 'ie'
end if
! set model parameter to FM type
CALL set_model_fm
!
end subroutine set_model
!================================================================================
subroutine set_model_fm
! set mass/temperature ratio to FM type
implicit none
!
CALL FM_ENTER_USER_ROUTINE
! Set temperature and mass ratios to FM type
sigmafm_ = TO_FM(sigma_)
taufm_ = TO_FM(tau_)
chifm_ = SQRT(taufm_/sigmafm_)
! Set perpendicular waves numbers to FM type
bafm_ = TO_FM(ba_)
bbfm_ = TO_FM(bb_)
!
!
! Set Sugama parameter
thetaabfm_ = SQRT(TO_FM( taufm_ + chifm_**2)/(TO_FM(1) + chifm_**2))
thetabafm_ = SQRT(TO_FM( TO_FM(1)/taufm_ + TO_FM(1)/chifm_**2)/(TO_FM(1) + TO_FM(1)/chifm_**2))
!
CALL FM_EXIT_USER_ROUTINE
end subroutine set_model_fm
!================================================================================
subroutine invert_model
! inverse ab --> ba
! Used in the Sugama collision operator
IMPLICIT NONE
!
CALL FM_ENTER_USER_ROUTINE
!
sigmafm_ = 1/sigmafm_
taufm_ = 1/taufm_
chifm_ = 1/chifm_
!
! Set Sugama parameter
thetaabfm_ = SQRT( (TO_FM(1) + sigmafm_)/(TO_FM(1) + sigmafm_/taufm_))
thetabafm_ = SQRT( (TO_FM(1) + TO_FM(1)/sigmafm_)/(TO_FM(1) + taufm_/sigmafm_))
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE invert_model
!================================================================================
SUBROUTINE operatormodel_readinput
! Read Model Operator from input file
USE BASIC
IMPLICIT NONE
! Namelist
NAMELIST /OPERATOR_MODEL/ ETEST,ITEST,EBACK,IBACK,ESELF,ISELF,GKE,GKI,DKTEST,DKBACK,ADDBACK,ADDTEST,only_gk_part
!
IF(me .ne. 0) GOTO 111
READ(lu_in,OPERATOR_MODEL)
WRITE(*,OPERATOR_MODEL)
!
111 CONTINUE
!
! MPI COMM
CALL MPI_bcast(ETEST,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(EBACK,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ITEST,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(IBACK,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ESELF,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ISELF,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(GKI,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(GKE,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(DKTEST,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(DKBACK,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ADDBACK,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast(ADDTEST,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
CALL MPI_bcast( only_gk_part,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
!
!
write(strETEST,'(I1)') ETEST
write(strITEST,'(I1)') ITEST
write(strEBACK,'(I1)') EBACK
write(strIBACK,'(I1)') IBACK
write(strESELF,'(I1)') ESELF
write(strISELF,'(I1)') ISELF
write(strGKE,'(I1)') GKE
write(strGKI,'(I1)') GKI
! Set different flags
! If Sugama COOP
IF(ETEST .eq. 3 .or. EBACK .eq. 3 .or. ITEST .eq. 3 .or. IBACK .eq. 3 .or. ESELF .eq. 3 .or. ISELF .eq. 3) THEN
IFSUGAMA = .true.
ENDIF
! If Improved Sugama COOP
IF(ETEST .eq. 5 .or. EBACK .eq. 5 .or. ITEST .eq. 5 .or. IBACK .eq. 5 .or. ESELF .eq. 5 .or. ISELF .eq. 5) THEN
IFSUGAMA = .true.
IFIMPROVEDSUGAMA = .true.
imaxx = JEmaxx
ENDIF
! IF zeroth order Hirshman-Sigmar-Clarke (0HSC)
IF(ETEST .eq. 6 .or. EBACK .eq. 6 .or. ITEST .eq. 6 .or. IBACK .eq. 6 .or. ESELF .eq. 6 .or. ISELF .eq. 6) THEN
IF0HSC = .true.
imaxx = JEmaxx
ENDIF
! IF GK abel
IF(ISELF .eq. 7 ) THEN
IFABEL = .true.
ENDIF
! IF Coulomb operator
IF(ETEST .eq. 1 .or. EBACK .eq. 1 .or. ITEST .eq. 1 .or. IBACK .eq. 1 .or. ESELF .eq. 1 .or. ISELF .eq. 1) THEN
IFCOULOMB = .true.
ENDIF
! IF gyrokinetic
IF(GKE .eq. 1 .or. GKI .eq. 1) THEN
IFGK = .true.
ENDIF
!
! IF GK Polarization
IF(ISELF .eq. 9) THEN
IFPOL = .true.
ENDIF
!
END SUBROUTINE operatormodel_readinput
END MODULE model
c4science · Help