Homec4science
Diffusion COSOlver (new_coulomb) Browse

cosolver/gkcoulomb_suite/Tljpmf/src/numerics_mod.f90
138972b69f9dnew_coulomb

Last Change
Raw File

numerics_mod.f90
View Options

module numerics
!
! Coulomb coefficient (COCO) module
!
!
use prec_const
use FMZM
use basic
use model
use array
use basis_transformation
implicit none
contains
! include coefficients and Coulomb speed functions
include 'coulomb_speed_functions.h'
!
FUNCTION ekab(k) RESULT(XOUT)
!
! compute int_0^\infty d x Erf'(s_a) s_b^{2k} e^{-x^2) [see Eq.(60) in Ji & Held (2006)]
!
! Returns 0 if k < 0
!
IMPLICIT NONE
INTEGER, INTENT(in) :: k
TYPE(FM) :: XOUT
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM('0.0')
!
IF( k .ge. 0) THEN
XOUT = FACTORIAL(TO_FM(2*k - 1 )/2)/FACTORIAL(TO_FM(-1)/2)/(TO_FM(1._dp + chifm_**2))**(real(k,dp) + 0.5_dp)
ENDIF
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION ekab
!________________________________________________________________
!
FUNCTION Erfkab(k) RESULT(XOUT)
!
! computes int_0^\infty d x Erf(s_a) s_b^{2k+1} e^{-x^2) [see Eq.(61) in Ji & Held (2006)]
!
!
! Note: Returns 0 if k < -1.
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: k
TYPE(FM) :: XOUT
TYPE(FM), SAVE:: ekab_
INTEGER :: i
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM('0.0')
ekab_ = TO_FM('0.0')
!
! ===== IF k > -2 =====
IF( k .ge. 0) THEN
DO i=0, k
ekab_ = ekab(i)
XOUT = XOUT + FACTORIAL(TO_FM(k))/FACTORIAL(TO_FM(i))*ekab_
END DO
!
XOUT = XOUT*chifm_/2
!
ENDIF
!
IF( k .eq. -1 ) THEN
XOUT = ASINH(chifm_)
ENDIF
!
! ======== NOT DEFINED RETURN 0 =====
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
end FUNCTION Erfkab
!________________________________________________________________
!
subroutine compute_coeff
! computes GK COCO coeff
!
implicit none
!
! local vars.
integer :: m,p,f,j,l
! last computed indices
integer :: m_lst,p_lst,f_lst,j_lst,l_lst
integer, parameter :: CSZ = 10 ! size of the routine cache
TYPE(FM), SAVE :: XOUT(CSZ)
integer :: idx_save(5,CSZ) ! indices array of XOUT
integer :: iNN_start,iNN,NN
CHARACTER(100) :: XOUT_str
INTEGER :: idx,i,istart
INTEGER :: fid,buff_
INTEGER :: ENDOF
LOGICAL :: restartisfound
! output file
character*40 :: filename_
!
CALL FM_ENTER_USER_ROUTINE
!
!
! Define resfile
! filename_ = trim(adjustl('Tljpmf_'//'L'//ldim_str//'J'//jdim_str//'P'//pdim_str//'M'//mdim_str//'F'//fdim_str))//'.csv'
!
filename_ = trim(adjustl(me_str))//'.csv'
!
!
! open file
OPEN(fid,FILE = filename_)
!
NN = 1
IF(restart) THEN
ENDOF = 0
NN = 1
DO WHILE (ENDOF .eq. 0)
READ(fid,*,IOSTAT=ENDOF) l_lst,j_lst,p_lst,m_lst,f_lst
NN = NN +1
ENDDO
!
! get last position in
!
restartisfound = .false.
istart = 1
DO WHILE(.not. restartisfound)
IF(idxNNCoeffs(istart,1) .eq. l_lst .and. idxNNCoeffs(istart,2) .eq. j_lst .and. idxNNCoeffs(istart,3) .eq. p_lst .and. idxNNCoeffs(istart,4) .eq. m_lst .and. idxNNCoeffs(istart,5) .eq. f_lst) THEN
restartisfound = .true.
iNN_start = istart + 1
ENDIF
istart = istart +1
ENDDO
!
IF(iNN_start .le. nbNN_l) THEN
write(*,*) me, ' restarts from ',iNN_start
ELSE
write(*,*) me, ' already finised !!'
ENDIF
ELSE
iNN_start = 1;
ENDIF
!
idx = 1
!
! main loop
DO iNN=iNN_start,nbNN_l
!
! Get indices of coefficients
l = idxNNCoeffs(iNN,1)
j = idxNNCoeffs(iNN,2)
p = idxNNCoeffs(iNN,3)
m = idxNNCoeffs(iNN,4)
f = idxNNCoeffs(iNN,5)
!
! save indices
idx_save(1,idx) = l
idx_save(2,idx) = j
idx_save(3,idx) = p
idx_save(4,idx) = m
idx_save(5,idx) = f
!
XOUT(idx) = Tabljpmf(l,j,p,m,f)
! output only the coefficient greater than 1E-50
!! IF(abs(XOUT(idx)) > TO_FM('1E-60')) THEN
!
idx = idx +1
! If the cache is full, empty it to the output file
if(idx > CSZ) then
! write to *.csv file
DO idx =1, CSZ
write(fid,"(5(I3,','),A)") idx_save(1:5,idx),&
trim(FM_FORMAT('E100.90',XOUT(idx)))
ENDDO
!
! Flush to ensure data is written
flush(fid)
!
! Reset idx to the first cache element
idx = 1
!
endif
!! ENDIF
ENDDO
! If elements remain in the cache, flush it.
if( idx>1 ) then
do i = 1, idx-1
write(fid,"(5(I3,','),A)") idx_save(1:5,i), trim(FM_FORMAT('E100.90',XOUT(i)))
end do
! Again, ensure data is written by flushing
flush(fid)
end if
!
! close file
CLOSE(fid)
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine compute_coeff
!
!________________________________________________________________
!
FUNCTION Tabljpmf(l,j,p,m,f) result(XOUT)
! compute Tabljpmf coefficient [see overleaf DOC]
implicit none
!
INTEGER,intent(in) :: l,j,p,m,f
TYPE(FM):: XOUT
TYPE(FM),SAVE :: T5inv,Lptd,nuF_,nuT_
INTEGER :: t,d
CHARACTER(len=T5len) :: T5str_
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM(0)
!
do t= 0,f + floor(real(l+m)/2)
T5inv = GET_T5inv(l,f,p,t,m)
IF(T5inv .ne. ZEROFM) THEN
do d =0,t
Lptd = Lpjl(real(p,dp),real(t,dp),real(d,dp))
if( coll_type .eq. 'self') THEN
nuT_ = nusTabpjgd(p,j,p,d)
nuF_ = nusFabpjgd(p,j,p,d)
XOUT = XOUT + Lptd*( nuT_ + nuF_)*T5inv
ELSEIF( coll_type .eq. 'test' ) THEN
nuT_ = nusTabpjgd(p,j,p,d)
XOUT = XOUT + Lptd*nuT_*T5inv
ELSEIF ( coll_type .eq. 'field' ) THEN
nuF_ = nusFabpjgd(p,j,p,d)
XOUT = XOUT + Lptd*nuF_*T5inv
ENDIF
enddo
ENDIF
enddo
!
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION Tabljpmf
!
end module numerics
c4science · Help