Homec4science
Diffusion FENNECS lite (master) Browse

fennecslite/src/neutcol_mod.f90
b8eece80d5a5master

Last Change
Raw File

neutcol_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: neutcol
!
!> @author
!> Guillaume Le Bars EPFL/SPC
!> Pierrick Giroud-Garampon EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for handling the electron-neutral collisions and creating electrons
!> by ionisation Based on the paper by Birdsall 1991 and Sengupta et al.
!------------------------------------------------------------------------------
module neutcol
USE constants
IMPLICIT NONE
private
LOGICAL, SAVE :: nlcol=.false. !< Flag to activate or not electron neutral collisions
LOGICAL, SAVE :: adjustv=.false. !< Flag to activate or not ajustment of the velocity
LOGICAL :: nlmaxwellio=.false. !< Flag to define how ionised electrons are created (physically or according to maxwellian)
INTEGER :: itio = 1 !< number of dt between injcetion by the ionisation source
INTEGER :: itfreq = 1 !< number of dt between each evaluation of cross-sections
Real(kind=db) :: neutdens=2.4e16 !< Neutral particle density in m-3
Real(kind=db) :: neuttemp=300 !< Neutral particle temperature in K
Real(kind=db) :: neutpressure !< Neutral particle pressure in mbar
Real(kind=db) :: scatter_fac = 24.2 !< Energy scattering factor for the considered gas (here for Ne) [eV] see Opal 1971 https://doi.org/10.1063/1.1676707
real(kind=db) :: Eion = 21.56 !< Ionisation energy (eV) (here for Ne)
Real(kind=db) :: E0 = 27.21 !< Atomic unit of energy used for calculation of deviation angles [eV]
Real(kind=db) :: dE = 1 !< resolution for the computation of the cross-sections
Real(kind=db) :: Emax = 1 !< resolution for the computation of the cross-sections
real(kind=db) :: collfactor !< Normalised collision factor (n_n \delta t)
Real(kind=db) :: Escale !< Energy normalisation factor used to reduce computation costs
INTEGER :: nb_io_cross=0
Real(kind=db), ALLOCATABLE :: io_cross_sec(:,:) !< Ionisation cross-section table
Real(kind=db), ALLOCATABLE :: io_cross_sec_lin(:,:) !< Ionisation cross-section table for linear interpolation
Real(kind=db), ALLOCATABLE :: io_growth_cross_sec(:) !< Ionisation exponential fitting factor
INTEGER :: nb_ela_cross=0
Real(kind=db), ALLOCATABLE :: ela_cross_sec(:,:) !< Elastic collision cross section table
Real(kind=db), ALLOCATABLE :: ela_cross_sec_lin(:,:) !< Elastic collision cross section table for linear interpolation
Real(kind=db), ALLOCATABLE :: ela_growth_cross_sec(:) !< Elastic collision exponential fitting factor
Real(kind=db), ALLOCATABLE :: vdistrib(:,:) !< Velocity distribution
CHARACTER(len=128) :: io_cross_sec_file=''
CHARACTER(len=128) :: ela_cross_sec_file=''
CHARACTER(len=128) :: vdistrib_file=''
Real(kind=db) :: etemp=22000 !< In case of nlmaxwelio, defines the temperature of created electrons [K]
Real(kind=db) :: vth !< In case of nlmaxwelio, defines the normalised thermal velocity of created electrons
LOGICAL :: nldragio=.false. !< Set if inpinging electrons are affected by ionising collisions
INTEGER :: species(2) !< species(1) contains the specie index in plist which stores the colliding particles, species(2) stores the specie index for the released ion.
LOGICAL :: isotropic = .true. !< is the scattering angle isotropic
type coll_process
Real(kind=db):: Emin !< minimum energy in EV for collision to happen
Real(kind=db):: Emax !< maximum energy in EV for linear interpolation
Real(kind=db):: deltaE !< energy step in ev for interpolation
Real(kind=db), allocatable:: cross_sec(:) !< cross-section as read from file
Real(kind=db), allocatable:: E_levels(:) !< energy values as read from file
Real(kind=db), allocatable:: growth_parameter(:) !< collision exponential fitting factor
Real(kind=db), allocatable:: lin_E(:) !< linear energy levels for linear interpolation
Real(kind=db), allocatable:: lin_cross(:) !< linear cross-sections for linear interpolation
end type
NAMELIST /neutcolparams/ neutdens, Eion, scatter_fac, nlcol, adjustv, io_cross_sec_file, ela_cross_sec_file, vdistrib_file, nlmaxwellio, etemp, nldragio, itio, itfreq, species, isotropic
PUBLIC:: neutcol_init, neutcol_diag, neutdens, comp_neutcol, ionization_inject, itio, itfreq
CONTAINS
subroutine neutcol_init(lu_in, p)
use mpi
Use basic, only: mpirank, dt, nlclassical,rnorm, vnorm, rgrid, zgrid, nr, nz, nnr, nnz, invdr, invdz, nsubr, nsubz
Use beam, only: particles
Use constants
implicit none
INTEGER, INTENT(IN) :: lu_in
TYPE(particles) :: p
INTEGER:: ierr, istat, i, j, k
character(len=1000) :: line
real(kind=db):: xsi
! Default
species(1)=1
species(2)=-1
Rewind(lu_in)
READ(lu_in, neutcolparams, iostat=istat)
if (istat.gt.0) then
backspace(lu_in)
read(lu_in,fmt='(A)') line
write(*,'(A)') &
'Invalid line in neutcolparams: '//trim(line)
call MPI_Abort(MPI_COMM_WORLD, -1, ierr)
stop
end if
IF(mpirank .eq. 0) THEN
WRITE(*, neutcolparams)
END IF
if(.not. nlcol) return
if(nlclassical)THEN
Escale=0.5*p%m/elchar*vlight**2
else
Escale=p%m*vlight**2/elchar
end if
if (nlmaxwellio) vth=sqrt(kb*etemp/p%m)/vnorm
if(io_cross_sec_file .ne.'') then
call read_cross_sec(io_cross_sec_file,io_cross_sec, nb_io_cross)
if(nb_io_cross .gt. 0) then
allocate(io_growth_cross_sec(nb_io_cross-1))
! Normalisations
io_cross_sec(:,2)=io_cross_sec(:,2)/rnorm**2
! Precomputing of exponential fitting factor for faster execution
io_growth_cross_sec=log(io_cross_sec(2:nb_io_cross,2)/io_cross_sec(1:nb_io_cross-1,2))/ &
& log(io_cross_sec(2:nb_io_cross,1)/io_cross_sec(1:nb_io_cross-1,1))
end if
end if
if(ela_cross_sec_file .ne.'') then
call read_cross_sec(ela_cross_sec_file,ela_cross_sec, nb_ela_cross)
if(nb_ela_cross .gt. 0) then
allocate(ela_growth_cross_sec(nb_ela_cross-1))
! Normalisations
ela_cross_sec(:,2)=ela_cross_sec(:,2)/rnorm**2
if(.not. isotropic) then
do i=1,nb_ela_cross
xsi=ela_cross_sec(i,1)/(0.25*E0+ela_cross_sec(i,1))
ela_cross_sec(i,2)=ela_cross_sec(i,2)*(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi))
end do
end if
! Precomputing of exponential fitting factor for faster execution
ela_growth_cross_sec=log(ela_cross_sec(2:nb_ela_cross,2)/ela_cross_sec(1:nb_ela_cross-1,2))/ &
& log(ela_cross_sec(2:nb_ela_cross,1)/ela_cross_sec(1:nb_ela_cross-1,1))
end if
END IF
IF(vdistrib_file .ne. '') THEN
CALL read_vdistrib(vdistrib_file,vdistrib)
END IF
nlcol=nlcol .and. (allocated(io_cross_sec) .or. allocated(ela_cross_sec))
adjustv=adjustv .and. allocated(vdistrib)
! Collision factor depending on neutral gas parameters
collfactor=neutdens*dt*rnorm**3*itio
neutpressure=neutdens*kb*300/100
end subroutine neutcol_init
Subroutine neutcol_diag(File_handle, str, vnorm)
use mpi
Use futils
Integer:: File_handle
Real(kind=db):: vnorm
Character(len=*):: str
CHARACTER(len=256):: grpname
Integer:: ierr, mpirank,i
Real(kind=db)::xsi
Real(kind=db),allocatable:: nonisotropic_rescale(:,:)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
IF(mpirank .eq. 0 .and. nlcol) THEN
Write(grpname,'(a,a)') trim(str),"/neutcol"
If(.not. isgroup(File_handle, trim(grpname))) THEN
CALL creatg(File_handle, trim(grpname))
END IF
Call attach(File_handle, trim(grpname), "neutdens", neutdens)
Call attach(File_handle, trim(grpname), "neuttemp", neuttemp)
Call attach(File_handle, trim(grpname), "neutpressure", neutpressure)
Call attach(File_handle, trim(grpname), "scatter_fac", scatter_fac)
Call attach(File_handle, trim(grpname), "Eion", Eion)
Call attach(File_handle, trim(grpname), "E0", E0)
Call attach(File_handle, trim(grpname), "Escale", Escale)
Call putarr(File_handle,trim(grpname)//"species", species)
if (allocated(io_cross_sec)) Call putarr(File_handle, trim(grpname)//"/io_cross_sec", io_cross_sec)
allocate(nonisotropic_rescale(nb_ela_cross,2))
nonisotropic_rescale=1
if(.not. isotropic) then
do i=1,nb_ela_cross
xsi=ela_cross_sec(i,1)/(0.25*E0+ela_cross_sec(i,1))
nonisotropic_rescale(i,2)=(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi))
end do
end if
if (allocated(ela_cross_sec)) Call putarr(File_handle, trim(grpname)//"/ela_cross_sec", ela_cross_sec/nonisotropic_rescale)
END IF
End subroutine neutcol_diag
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Simulates the ionisation process
!---------------------------------------------------------------------------
SUBROUTINE comp_neutcol(p)
USE random_mod
use omp_lib
USE basic
USE beam, ONLY: particles
INTEGER:: i, j, iend, k, omp_thread, num_threads
real(kind=db):: v2, v, Ek, sig_io, sig_ela, Uperp
REAL(kind=db) :: rand(0:24)
TYPE(particles) :: p
IF(.not. nlcol .or. p%Nploc .le. 0) RETURN
num_threads=omp_get_max_threads()
omp_thread=omp_get_thread_num()+1
!$OMP DO
DO j=1,p%Nploc,25
! for each particle
iend=min(p%Nploc,j+24)
CALL random_array(rand,iend-j+1,ran_index(omp_thread),ran_array(:,omp_thread))
! we calculate the kinetic energy and norm of the velocity
DO i=j,iend
v2=(p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2)
IF(adjustv) THEN
Uperp=sqrt(v2-(p%Upar(i))**2)
DO k=1,size(vdistrib,1)-1
IF(rand(i-j) .ge. vdistrib(k,2) .and. rand(i-j) .lt. vdistrib(k+1,2)) THEN
Uperp=Uperp*vdistrib(k,1)
EXIT
END IF
END DO
v2=(Uperp**2+p%Upar(i)**2)
IF(nlclassical) THEN
Ek=v2*Escale
v=sqrt(v2)
ELSE
Ek=(p%gamma(i)-1)*Escale
v=sqrt(v2)/p%gamma(i)
END IF
sig_io=0
sig_ela=0
! computes the ionisation and elastic collision frequencies at this kinetic energy
IF(Ek .gt. Eion .and. nb_io_cross .gt. 1) THEN
sig_io=sig_fit(io_cross_sec,io_growth_cross_sec,ek, nb_io_cross)
END IF
IF(nb_ela_cross .gt. 1) THEN
sig_ela=sig_fit(ela_cross_sec,ela_growth_cross_sec,ek, nb_ela_cross)
END IF
p%ionfreq(i)=sig_io*v*vnorm*rnorm**2*neutdens
p%collfreq(i)=(sig_ela+sig_io)*v*vnorm*rnorm**2*neutdens
ELSE
IF(nlclassical) THEN
Ek=v2*Escale
v=sqrt(v2)
ELSE
Ek=(p%gamma(i)-1)*Escale
v=sqrt(v2)/p%gamma(i)
END IF
p%ionfreq(i)=(-8.996644e-15+7.225088e-21*v*vnorm-3.693399e-28*(v*vnorm)**2)*neutdens
p%collfreq(i)=(7.942905e-14-8.886946e-21*v*vnorm+4.948035e-28*(v*vnorm)**2)*neutdens
END IF
END DO
END DO
!$OMP END DO NOWAIT
END SUBROUTINE comp_neutcol
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Simulates the ionization process
!---------------------------------------------------------------------------
SUBROUTINE ionization_inject(plist)
! Inject the particles created by ionization of the neutrals
USE random_mod
USE omp_lib
USE basic
USE beam
IMPLICIT NONE
TYPE(particles), TARGET :: plist(:)
TYPE(particles), POINTER :: p
TYPE(linked_part_row) :: ins_p
TYPE(linked_part), POINTER :: created
INTEGER :: npartsadd, nbionisation, i, j, k, iend, nbelec, omp_thread, num_threads
REAL(kind=db) :: randnb(1:1)
REAL(kind=db) :: lambda
REAL(kind=db), DIMENSION(15) :: proba, CDF
REAL(kind=db), ALLOCATABLE:: R(:),Z(:)
p=>plist(species(1))
IF(.not. nlcol .or. p%Nploc .le. 0) RETURN
nbionisation=0
num_threads=omp_get_max_threads()
omp_thread=omp_get_thread_num()+1
ALLOCATE(ins_p%start)
ins_p%n=0
created=>ins_p%start
!$OMP DO
DO j=1,p%Nploc,25
! for each particle
iend=min(p%Nploc,j+24)
DO i=j,iend
IF(p%ionfreq(i) .gt. 0) THEN
lambda=p%ionfreq(i)*dt*tnorm*itio ! number of electrons to add
! Compute probabilities and CDF values for the Poisson distribution of parameter lambda
DO k = 1,15
proba(k) = exp(-lambda)*lambda**(k-1)/factorial_fun(k-1);
CDF(k) = sum(proba(1:k));
END DO
! Generate the number of electrons according to Poisson distribution
CALL random_array(randnb,1,ran_index(omp_thread),ran_array(:,omp_thread)) ! uniform distribution in [0,1]
DO k = 1,size(CDF)-1
IF(randnb(1) .lt. CDF(1)) THEN
nbelec = 0
ELSE IF((CDF(k) .le. randnb(1)) .and. (randnb(1) .lt. CDF(k+1))) THEN
nbelec = k
END IF
END DO
IF(nbelec .le. 0) THEN
CYCLE
ELSE
! nbelec electrons are added in the simulation by the source
nbionisation=nbionisation+nbelec
DO k=1,nbelec
! Prepare the memory for the released electron
ins_p%n=ins_p%n+1
ALLOCATE(created%next)
created%next%prev=>created
! Fill the created particle
IF(larmorrad) THEN
CALL random_array(randnb,1,ran_index(omp_thread),ran_array(:,omp_thread)) ! uniform distribution in [0,1]
created%p%pos=(/p%pos(1,i)-p%rho(i)*cos(2*pi*randnb),p%pos(2,i),p%pos(3,i)/)
ELSE
created%p%pos=(/p%pos(1,i),p%pos(2,i),p%pos(3,i)/)
END IF
created%p%rho=0
created%p%U=0
created%p%Upar=0
created%p%gamma=1.0
! Prepare the next created particle
ins_p%end=>created
created=>created%next
END DO
END IF
END IF
END DO
END DO
!$OMP END DO NOWAIT
! Clean up the memory after the loop
IF(associated(created%prev)) THEN
created=>created%prev
ins_p%end=>created
DEALLOCATE(created%next)
ELSE
DEALLOCATE(ins_p%start)
END IF
!$OMP CRITICAL
p%nbionisation=p%nbionisation+nbionisation
!$OMP END CRITICAL
!$OMP BARRIER
! We collect all created particules into one linked list for easier insertion in plist
! Insert ions
!$OMP CRITICAL
IF(species(2) .gt. 0 .and. ins_p%n .gt. 0) THEN
CALL add_created_part(plist(species(2)),ins_p,.false.,.true.)
END IF
!$OMP END CRITICAL
! Insert electrons
!$OMP CRITICAL
IF(ins_p%n .gt.0) THEN
CALL add_created_part(plist(species(1)),ins_p,.true.,.false.)
END IF
!$OMP END CRITICAL
END SUBROUTINE ionization_inject
FUNCTION sig_fit(sig_vec,growth_vec,ek,nb_cross)
use distrib, ONLY: closest
real(kind=db)::sig_fit, ek
real(kind=db):: sig_vec(:,:), growth_vec(:)
Integer:: k, nb_cross
sig_fit=0
k=closest(sig_vec(:,1),ek, nb_cross-1)
if(k.lt.1) return
! Exponential fitting relevant at high energies
sig_fit=sig_vec(k,2)*exp(log(ek/sig_vec(k,1))*growth_vec(k))
END FUNCTION sig_fit
FUNCTION sig_fit_lin(sig_vec,growth_vec,ek,nb_cross)
use distrib, ONLY: closest
real(kind=db)::sig_fit_lin, ek
real(kind=db):: sig_vec(:,:), growth_vec(:)
Integer:: k, nb_cross
sig_fit_lin=0
k=(k-sig_vec(k,1))/dE
if(k.lt.1) return
! Exponential fitting relevant at high energies
sig_fit_lin=sig_vec(k,2)+(ek-sig_vec(k,1))/(sig_vec(k+1,1)-sig_vec(k,1))*(sig_vec(k+1,2)-sig_vec(k,2))
END FUNCTION sig_fit_lin
SUBROUTINE read_cross_sec(filename,cross_sec, nb_cross)
CHARACTER(len=*) ::filename
Real(kind=db), ALLOCATABLE :: cross_sec(:,:)
INTEGER:: nb_cross
INTEGER :: lu_cross_sec=9999
INTEGER:: i, openerr, reason
CHARACTER(len=256) :: header
real(kind=db):: t1,t2
nb_cross=0
OPEN(UNIT=lu_cross_sec,FILE=trim(filename),ACTION='READ',IOSTAT=openerr)
header=' '
IF(openerr .ne. 0) THEN
CLOSE(unit=lu_cross_sec)
RETURN
END IF
! The cross section table is defined as a two column energy and cross_section
DO WHILE(.true.)
READ(lu_cross_sec,'(a)',IOSTAT=reason) header
header=adjustl(header)
if(reason .lt. 0 ) exit ! We reached end of file
if( header(1:1) .ne. '!') then
READ(header,*) t1, t2
if(t1 .ne. 0 .and. t2.ne. 0) nb_cross=nb_cross+1
end if
END DO
if (allocated(cross_sec)) deallocate(cross_sec)
allocate(cross_sec(nb_cross,2))
REWIND(lu_cross_sec)
! The cross section table is defined as a two column energy and cross_section
i=1
DO WHILE(i .le. nb_cross)
READ(lu_cross_sec,'(a)',IOSTAT=reason) header
header=adjustl(header)
if(reason .lt. 0 ) exit ! We reached end of file
if( header(1:1) .ne. '!') then
READ(header,*) cross_sec(i,1), cross_sec(i,2)
if(cross_sec(i,1) .ne. 0 .and. cross_sec(i,2).ne. 0) i=i+1
end if
END DO
CLOSE(unit=lu_cross_sec)
END subroutine read_cross_sec
subroutine init_cross_sec_type(cross_sec_data, nb, cross_sec_array, exp_growth_array, Emin, Emax, deltaE)
type(coll_process):: cross_sec_data
INTEGER:: nb
real(kind=db):: cross_sec_array(:,:), exp_growth_array(:)
real(kind=db):: Emin, Emax, deltaE
allocate(cross_sec_data%cross_sec(nb),cross_sec_data%E_levels(nb))
cross_sec_data%cross_sec=cross_sec_array(:,2)
cross_sec_data%E_levels=cross_sec_array(:,1)
allocate(cross_sec_data%growth_parameter(nb-1))
cross_sec_data%growth_parameter=exp_growth_array
cross_sec_data%Emin=Emin
cross_sec_data%Emax=Emax
end subroutine
SUBROUTINE read_vdistrib(filename,vdistrib)
CHARACTER(len=*) ::filename
REAL(kind=db), ALLOCATABLE :: vdistrib(:,:)
INTEGER :: lu_vdistrib=1111
INTEGER:: i, openerr, reason, nbv
CHARACTER(len=256) :: header
nbv=0
OPEN(UNIT=lu_vdistrib,FILE=trim(filename),ACTION='READ',IOSTAT=openerr)
header=' '
IF(openerr .ne. 0) THEN
CLOSE(unit=lu_vdistrib)
RETURN
END IF
DO WHILE(.true.)
READ(lu_vdistrib,'(a)',IOSTAT=reason) header
header=adjustl(header)
IF(reason .lt. 0 ) EXIT ! We reached end of file
IF(header(1:1) .ne. '!') THEN
nbv=nbv+1
END IF
END DO
IF(allocated(vdistrib)) DEALLOCATE(vdistrib)
ALLOCATE(vdistrib(nbv,2))
REWIND(lu_vdistrib)
i=1
DO WHILE(i .le. nbv)
READ(lu_vdistrib,'(a)',IOSTAT=reason) header
header=adjustl(header)
IF(reason .lt. 0 ) EXIT ! We reached end of file
IF(header(1:1) .ne. '!') THEN
READ(header,*) vdistrib(i,1), vdistrib(i,2)
i=i+1
END IF
END DO
CLOSE(unit=lu_vdistrib)
END subroutine read_vdistrib
!---------------------------------------------------------------------------
! FUNCTION factorial_fun(n)
!
! DESCRIPTION
!> Gives the factorial of an integer
!--------------------------------------------------------------------------
INTEGER FUNCTION factorial_fun(n)
INTEGER :: ii
INTEGER :: n
factorial_fun = 1
IF (n .lt. 0) THEN
RETURN
ELSE IF (n .eq. 0) THEN
factorial_fun = 1
ELSE
DO ii=1,n
factorial_fun = factorial_fun*ii
END DO
END IF
END FUNCTION factorial_fun
end module neutcol
c4science · Help