Homec4science
Diffusion COSOlver (master) Browse

cosolver/src/gkcoulomb_coeffs_mod.f90
5b66a09e0f20master

Last Change
Raw File

gkcoulomb_coeffs_mod.f90
View Options

module gkcoulomb_coeffs
! loads the pre-computed coefficients for gkcoulomb operator: Tljpmf, curlyNpjmn, Pipjm
! Note: - used for GK Coulomb collision operator
use prec_const
use basic
use model
use array
use coeff
implicit none
PUBLIC
INTEGER, SAVE :: fid_
INTEGER, SAVE :: IDXMAX ! ... maximum indices for Tljpmf
! Tljpmf maximum indices
INTEGER, SAVE :: Tidx_l_max, Tidx_j_max, Tidx_p_max, Tidx_f_max
! CurlyNpjmn maximum indices
INTEGER, SAVE :: Nidx_p_max, Nidx_j_max, Nidx_n_max
CHARACTER(len=*),PARAMETER :: Tljpmf_in = 'log.in'
CHARACTER(len=*),PARAMETER :: CurlyNpjmn_in = 'log.in'
CHARACTER(len=*),PARAMETER :: Pipjm_in = 'Pipjm.in'
CHARACTER(len=*), PARAMETER :: Tljpmf_dir = 'gkcoulomb_suite'
CHARACTER(len=*), PARAMETER :: CurlyNpjmn_dir = 'gkcoulomb_suite'
contains
!
!___________________________________________________________
!
subroutine load_gk_coulomb_coeffs(dirCab)
! load the Tljpmf and curlyNpjmn coeffs precomputed for gkcoulomb
implicit none
! directory where the coefficients are stored
CHARACTER(len=*), intent(in) :: dirCab ! species depend directory
IF( me .eq. 0) WRITE(*,*) 'Load GK Coulomb coefficients: ', dirCab, ' ...'
CALL load_Tljpmf(dirCab)
call load_CurlyNpjmn
CALL eval_ALL2L
IF(me .eq. 0) WRITE(*,*) '... done'
CALL mpi_barrier(MPI_COMM_WORLD,ierr)
end subroutine load_gk_coulomb_coeffs
!
!___________________________________________________________
!
subroutine load_Tljpmf(dirCab)
!
! read NN Tljpmf coefficients from Tljpmf directory dir
!
implicit none
!
! directory where the coefficients are stored
CHARACTER(len=*), intent(in) :: dirCab ! species depend directory
! locs. vars.
integer :: NN,ENDOF,iNN
CHARACTER(len=llen ) :: buff_,filename_
INTEGER :: IT1,IT2,IT3,IT4,IT5
INTEGER :: FMAXX
CHARACTER(len=100),DIMENSION(:),ALLOCATABLE :: TableNames
INTEGER :: Ntab,iTab
CHARACTER(len=FMlen) :: buffc_
!
!
! 1. Get the number of NN coefficients
!
FMAXX = max(Jmaxe,Jmaxi) + max(neFLR,niFLR)
IDXMAX = max(FMAXX,JEmaxx,PMMAXX,PMAXI,PMAXE)+1
!
WRITE(filename_, *) trim(adjustl(cosolver_dir))//'/'// Tljpmf_dir // '/' // dirCab // '/'// Tljpmf_in
!
IF( me .eq. 0) WRITE(*, *) 'Read log.in file ', trim(adjustl(filename_))
!
OPEN(fid_,FILE=trim(adjustl(filename_)))
! Count the number of files
NTab = 0
ENDOF = 0
DO WHILE(ENDOF .eq. 0)
READ(fid_,*,IOSTAT =ENDOF) buff_
NTab = NTab +1
ENDDO
!
NTab = NTab -1
!
CLOSE(fid_)
!
! Allocate array of Tljpmf table
IF( allocated(TableNames)) deallocate(TableNames)
ALLOCATE(TableNames(1:NTab))
!
!
! Get the name of the table containing T5 coeff.
!
OPEN(fid_ ,FILE=trim(adjustl(filename_)))
ENDOF = 0
iTab = 1
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) buff_
IF (iTab .le. NTab) THEN
TableNames(iTab) = adjustl(trim(buff_))
iTab = iTab +1
ENDIF
ENDDO
CLOSE(fid_)
!
!
Tidx_l_max = 0
Tidx_j_max = 0
Tidx_p_max = 0
Tidx_f_max = 0
!
NN = 0
DO iTab=1,NTab
! name of table to read
write(filename_, *) trim(adjustl(cosolver_dir))//'/'// Tljpmf_dir // '/' // dirCab // '/'// TableNames(iTab)
OPEN(fid_,FILE = trim(adjustl(filename_)))
ENDOF = 0
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) IT1,IT2,IT3,IT4,IT5
IF (IT1 .le. IDXMAX .and. IT2 .le. IDXMAX .and. IT3 .le. IDXMAX .and. IT4 .le. IDXMAX .and. IT5 .le. IDXMAX) THEN
IF( IT1 .gt. Tidx_l_max ) Tidx_l_max = IT1
IF( IT2 .gt. Tidx_j_max ) Tidx_j_max = IT2
IF( IT3 .gt. Tidx_p_max ) Tidx_p_max = IT3
IF( IT5 .gt. Tidx_f_max ) Tidx_f_max = IT5
NN = NN+1
ENDIF
ENDDO
CLOSE(fid_)
ENDDO
!
! 2. Allocate arrays
!
! revert the indices (column-major): f, m, p, j, l
if( allocated(Tljpmf_STR)) deallocate(Tljpmf_STR)
! allocate( Tljpmf_STR(0: max(Pmaxi, Pmaxe) ,0:JEmaxx, 0: PMMAXX, 0:PMMAXX, 0 : FMAXX )) ! l, j, p , m, f
allocate( Tljpmf_STR( 0 : FMAXX, 0:PMMAXX, 0: PMMAXX , 0:JEmaxx, 0:max(Pmaxi, Pmaxe))) ! f, m, p, j, l
!
! 3. Inquire NN coefficients
!
iNN = 1
!
DO iTab =1, NTab
write(filename_, *) trim(adjustl(cosolver_dir))//'/'// Tljpmf_dir // '/' // dirCab // '/'// TableNames(iTab)
OPEN(FID_,FILE = filename_)
ENDOF = 0
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) IT1,IT2,IT3,IT4,IT5, buffc_
IF (IT1 .le. max(Pmaxi, Pmaxe) .and. &
IT2 .le. JEmaxx .and. &
IT3 .le. PMMAXX .and. &
IT4 .le. PMMAXX .and. &
IT5 .le. FMAXX ) THEN
Tljpmf_STR(IT5, IT4, IT3, IT2, IT1) = buffc_
IF(ENDOF .eq. 0) iNN = iNN+1
ENDIF
ENDDO
!
CLOSE(fid_)
ENDDO
! Free memory
DEALLOCATE(TableNames)
!
end subroutine load_Tljpmf
!
!______________________________________________________________________________________________
!
SUBROUTINE GET_NNLJPMF(l,j,p,m,f,NNLJPMF_)
!
! construct character string l//$//j//$//p//$//m//$//f
!
IMPLICIT NONE
!
INTEGER,intent(in) :: l,j,p,m,f
CHARACTER(len=NNLJPMFlen),intent(inout) :: NNLJPMF_
!
CHARACTER*2 :: SIT1,SIT2,SIT3,SIT4,SIT5
!
! Convert int to string
WRITE(SIT1,'(i2)') l
WRITE(SIT2,'(i2)') j
WRITE(SIT3,'(i2)') p
WRITE(SIT4,'(i2)') m
WRITE(SIT5,'(i2)') f
!
NNLJPMF_ = trim(adjustl(SIT1))//'$'//trim(adjustl(SIT2))//'$'//trim(adjustl(SIT3))//'$'//trim(adjustl(SIT4))//'$'//trim(adjustl(SIT5))
!
END SUBROUTINE GET_NNLJPMF
!______________________________________________________________________________________________
!
SUBROUTINE GET_NNPJM(p,j,m,NNPJM_)
!
! construct character string p//$//j//$//m
!
IMPLICIT NONE
!
INTEGER,intent(in) :: p,j,m
CHARACTER(len=NNPJMlen),intent(inout) :: NNPJM_
!
CHARACTER*2 :: SIT1,SIT2,SIT3
!
! Convert int to string
WRITE(SIT1,'(i2)') p
WRITE(SIT2,'(i2)') j
WRITE(SIT3,'(i2)') m
!
NNPJM_ = trim(adjustl(SIT1))//'$'//trim(adjustl(SIT2))//'$'//trim(adjustl(SIT3))
!
END SUBROUTINE GET_NNPJM
!______________________________________________________________________________________________
!
FUNCTION GET_Tljpmf(l,j,p,m,f) RESULT(Tljpmfstr_)
!
! Returns the Tljpmf string from character string array Tljpmfstr
!
! This routines can be obtained by converting the present implementation into an array.
IMPLICIT NONE
!
INTEGER,intent(in) :: l,j,p,m,f
CHARACTER(len=FMlen) :: Tljpmfstr_
! Local variables
CHARACTER(len=NNLJPMFlen) :: NNLJPMF_
INTEGER :: POS,i,inquireMaxInd
!
!
IF( l .gt. Tidx_l_max .or. j .gt. Tidx_j_max .or. p .gt. Tidx_p_max .or. f .gt. Tidx_f_max) THEN
WRITE(*,* ) 'ERROR Tljpmf: inquire unexisting indices', l,'/' , Tidx_l_max, ',',j,',', Tidx_j_max, ', ', p , '/', Tidx_p_max, ',',f ,'/' ,Tidx_f_max
STop
ENDIF
POS =0
! Get string id
CALL GET_NNLJPMF(l,j,p,m,f,NNLJPMF_)
! Get Tljpmf position
DO i=1,size(NNLJPMF)
IF( NNLJPMF(i) == NNLJPMF_) THEN
POS = i
GOTO 55
ENDIF
ENDDO
55 CONTINUE
! Retrieve coeff string from array
IF( POS .gt. 0) THEN
Tljpmfstr_ = TljpmfSTR(POS)
ELSE
Tljpmfstr_ = '0.0'
ENDIF
END FUNCTION GET_Tljpmf
!
!____________________________________________________________________________________________
!
subroutine read_Pipjm
!
! read NN Pipjm coefficients
!
implicit none
!
integer :: NN,ENDOF,iNN
CHARACTER(len=FMlen) :: buff_,filename_
INTEGER :: IT1,IT2,IT3
CHARACTER(len=FMlen),DIMENSION(:),ALLOCATABLE :: PipjmTableNames
INTEGER :: Ntab,iTab
CHARACTER(len=FMlen) :: coeffstr
!
!
! 1. Get the number of NN coefficients
!
!
OPEN(fid_,FILE=Pipjm_in)
! Count the number of files
NTab = 0
ENDOF = 0
DO WHILE(ENDOF .eq. 0)
READ(fid_,*,IOSTAT =ENDOF) buff_
NTab = NTab +1
ENDDO
!
NTab = NTab -1
!
CLOSE(fid_)
!
! Allocate array of Tljpmf table)
ALLOCATE(PipjmTableNames(1:NTab))
!
!
!
OPEN(fid_,FILE=Pipjm_in)
ENDOF = 0
iTab = 1
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) buff_
IF (iTab .le. NTab) THEN
PipjmTableNames(iTab) = adjustl(trim(buff_))
iTab = iTab +1
ENDIF
ENDDO
CLOSE(fid_)
!
! 3. Inquire NN coefficients
!
iNN = 1
!
DO iTab =1, NTab
filename_ = trim(adjustl('./Pipjm/'//PipjmTableNames(iTab)))
OPEN(FID_,FILE = filename_)
ENDOF = 0
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) IT1,IT2,IT3,coeffstr
IF (IT1 .le. PMMAXX .and. IT2 .le. JEmaxx .and. IT3 .le. PMMAXX) THEN
PipjmSTR(IT1,IT2,IT3) = coeffstr
IF(ENDOF .eq. 0) iNN = iNN+1
ENDIF
ENDDO
!
CLOSE(fid_)
ENDDO
!
!
!
! Free memory
DEALLOCATE(PipjmTableNames)
!
end subroutine read_Pipjm
!______________________________________________________________________________________________
!
subroutine load_curlyNpjmn
! load GK row Vector Npj from *.csv file to character array
! Note: only important the non-zero coefficients
implicit none
!
! directory where the coefficients are stored
CHARACTER(len=4) :: dirCab ! species depend directory
! local vars.
INTEGER :: ENDOF
INTEGER :: IT1,IT2,IT3,IT4,IT5
INTEGER :: PMAXX_,MMAXX_,NMAXX_,ROWMAXX_
CHARACTER(len=100),DIMENSION(:),ALLOCATABLE :: TableNames
INTEGER :: Ntab,iTab,NN,iNN
CHARACTER(len=FMlen) :: buff_
CHARACTER(len=llen) :: filename_
!
! Maximum indices to load
MMAXX_ = PMAXX
NMAXX_ = NPFLRMAXX
ROWMAXX_ = numbb(Pmaxb,Jmaxb)
!
if( allocated( curlyNpjmn_STR)) deallocate(curlyNpjmn_STR)
!
! revert the indices (column-major order): icol, np, m, j, p
! ion
if( coll_ab .eq. 'ie' .or. coll_ab .eq. 'ii' ) THEN
CALL allocate_array(curlyNpjmn_STR,numbi(0,0),numbi(Pmaxi,Jmaxi),0,NPFLRMAXX, 0,PMMAXX, 0,JEmaxx, 0,PMMAXX )
! CALL allocate_array(curlyNpjmn_STR,0,PMMAXX,0,JEmaxx,0,PMMAXX,0,NPFLRMAXX,numbi(0,0),numbi(Pmaxi,Jmaxi))
dirCab = 'Nipj'
ENDIF
! electron
IF( coll_ab .eq. 'ei' .or. coll_ab .eq. 'ee' ) THEN
dirCab = 'Nepj'
CALL allocate_array(curlyNpjmn_STR,numbe(0,0),numbe(Pmaxe,Jmaxe),0,NPFLRMAXX, 0,PMMAXX, 0,JEmaxx, 0,PMMAXX )
! CALL allocate_array(curlyNpjmn_STR,0,PMMAXX,0,JEmaxx,0,PMMAXX,0,NPFLRMAXX,numbe(0,0),numbe(Pmaxe,Jmaxe))
ENDIF
!
WRITE(filename_, *) trim(adjustl(cosolver_dir))//'/' //CurlyNpjmn_dir // '/' // dirCab // '/'// curlyNpjmn_in
!
IF( me .eq. 0) WRITE(*, *) 'Load moments ... '
IF( me .eq. 0) WRITE(*, *) 'Read log.in file ', trim(adjustl(filename_))
!
OPEN(fid_,FILE = trim(adjustl(filename_)))
!
! 1. Count the number of files and coeff
!
NTab = 0
ENDOF = 0
DO WHILE(ENDOF .eq. 0)
READ(fid_,*,IOSTAT =ENDOF) buff_
NTab = NTab +1
ENDDO
!
NTab = NTab -1
!
CLOSE(fid_)
!
! Allocate array of
IF( allocated(TableNames)) deallocate(TableNames)
ALLOCATE(TableNames(1:NTab))
!
OPEN(fid_,FILE= trim(adjustl(filename_)))
!
ENDOF = 0
iTab = 1
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) buff_
IF (iTab .le. NTab) THEN
TableNames(iTab) = adjustl(trim(buff_))
iTab = iTab +1
ENDIF
ENDDO
CLOSE(fid_)
!
Nidx_p_max = 0
Nidx_j_max = 0
Nidx_n_max = 0
!
iNN =1
DO iTab=1,NTab
! name of table to read
write(filename_, *) trim(adjustl(cosolver_dir))//'/' // curlyNpjmn_dir // '/' // dirCab // '/'// TableNames(iTab)
OPEN(fid_,FILE = filename_)
ENDOF = 0
DO WHILE (ENDOF .eq. 0 )
READ(fid_,*,IOSTAT=ENDOF) IT1,IT2,IT3,IT4,IT5,buff_
IF( IT1 .gt. Nidx_p_max ) Nidx_p_max = IT1
IF( IT2 .gt. Nidx_j_max ) Nidx_j_max = IT2
IF( IT4 .gt. Nidx_n_max ) Nidx_n_max = IT4
IF (IT1 .le. PMMAXX .and. IT2 .le. JEmaxx .and. IT3 .le. MMAXX_ .and. IT4 .le. NMAXX_ .and. IT5 .le. ROWMAXX_) THEN
curlyNpjmn_STR(IT5,IT4,IT3,IT2,IT1) = buff_
ENDIF
ENDDO
CLOSE(fid_)
ENDDO
!
IF( PMMAXX .gt. Nidx_p_max .or. JEmaxx .gt. Nidx_j_max .or. NMAXX_ .gt. Nidx_n_max ) THEN
WRITE(*,*) 'Error CurlNpjmn: inquire unexisting indices', PMMAXX , '/', Nidx_p_max, ', ', JEmaxx, '/', Nidx_j_max, ', :,', NMAXX_, '/', Nidx_n_max
STOP
ENDIF
!
! Free memory
DEALLOCATE(TableNames)
!
end subroutine load_curlyNpjmn
!______________________________________________________________________________________________
!
subroutine GET_PJMNCOL(p,j,m,n,col,PJMNCOL_)
!
! construct character string p//$//j//$//m//$//n//$//col
!
implicit none
!
INTEGER,intent(in) :: p,j,m,n,col
CHARACTER(len=PJMNCOLlen),intent(inout) :: PJMNCOL_
!
CHARACTER*2 :: SIT1,SIT2,SIT3,SIT4,SIT5
! Convert int to string
WRITE(SIT1,'(i2)') p
WRITE(SIT2,'(i2)') j
WRITE(SIT3,'(i2)') m
WRITE(SIT4,'(i2)') n
WRITE(SIT5,'(i2)') col
!
PJMNCOL_ = trim(adjustl(SIT1))//'$'//trim(adjustl(SIT2))//'$'//trim(adjustl(SIT3))//'$'//trim(adjustl(SIT4))//'$'//trim(adjustl(SIT5))
!
end subroutine get_PJMNCOL
!______________________________________________________________________________________________
!
subroutine eval_ALL2L
! precomputes ALL2L coefficients before entering the prec_GK_coulomb subroutine
implicit none
INTEGER :: ff,mm,kk,nn
character(len= FMlen) :: val_str
! allocate array
call allocate_array(array_ALL2L,0, max(Jmaxi,Jmaxe) + max(niFLR,neFLR),0,PMMAXX,0,max(Jmaxe,Jmaxi),0,max(niFLR,neFLR))
!array_ALL2L = ZEROFM
CALL FM_ENTER_USER_ROUTINE
DO nn=0,max(neFLR,niFLR)
DO kk=0,max(Jmaxe,Jmaxi)
DO mm=0,PMMAXX
DO ff=0,(max(Jmaxi,Jmaxe) + max(niFLR,neFLR))
CALL FMFORM('F60.40',ALL2L(nn,kk,ff,mm), val_str)
array_ALL2L(ff,mm,kk,nn) = val_str
ENDDO
ENDDO
ENDDO
ENDDO
CALL FM_EXIT_USER_ROUTINE
end subroutine eval_ALL2L
!
end module gkcoulomb_coeffs
c4science · Help