Homec4science
Diffusion FENNECS lite (3Dmpimode) Browse

fennecslite/src/neutcol_mod copy.f90
6b37315661933Dmpimode

Last Change
Raw File

neutcol_mod copy.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: neutcol
!
!> @author
!> Guillaume Le Bars 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 :: nlmaxwellio=.false. !< Flag to define how ionised electrons are created (physically or according to maxwellian)
INTEGER :: itcol = 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) :: rlimits(2) !< radial limits of the ionisation source
REAL(kind=db) :: zlimits(2) !< axial limits of the ionisation source
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)
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) :: Escale !< Energy normalisation factor used to reduce computation costs
CHARACTER(len=128) :: io_cross_sec_file=''
CHARACTER(len=128) :: ela_cross_sec_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=.true. !< 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 = .false. !< 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, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, &
& nldragio, itcol, itfreq, species, isotropic, rlimits, zlimits
PUBLIC:: neutcol_init, neutcol_diag, itcol, neutdens, comp_neutcol, ionisation_inject, itfreq
CONTAINS
subroutine neutcol_init(lu_in, p)
use mpi
Use basic, only: mpirank, dt, nlclassical,rnorm, vnorm
Use beam, only: particles
Use constants
implicit none
INTEGER, INTENT(IN) :: lu_in
TYPE(particles) :: p
INTEGER:: ierr, istat, i
character(len=1000) :: line
real(kind=db):: xsi
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
nlcol=nlcol .and. (allocated(io_cross_sec) .or. allocated(ela_cross_sec))
! Collision factor depending on neutral gas parameters
collfactor=neutdens*dt*rnorm**3*itcol
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
USE basic
INTEGER:: i,j
real(kind=db):: vEB, Ek
DO i=1,size(Er)
! for each grid cell
! we calculate the norm of the ExB velocity and the associated kinetic energy
vEB=abs(Er(i)/Bz(i)*vnorm) !((Er(i)*Bz(i))**2+(Ez(i)*Br(i))**2)/(Bz(i)**2+Br(i)**2)**2*vnorm**2
Ek=0.5*partmass/elchar*vEB**2
! computes the ionisation and elastic collision cross-sections at this kinetic energy
if(Ek .gt. Eion .and. nb_io_cross .gt. 1) then
ionfreq(i)=sig_fit(io_cross_sec,io_growth_cross_sec,Ek, nb_io_cross)*vEB*rnorm**2*neutdens
end if
if (nb_ela_cross .gt. 1) then
collfreq(i)=sig_fit(ela_cross_sec,ela_growth_cross_sec,Ek, nb_ela_cross)*vEB*rnorm**2*neutdens
end if
END DO
END SUBROUTINE comp_neutcol
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Simulates the ionisation process
!---------------------------------------------------------------------------
subroutine ionisation_inject(plist)
! Inject the particles in the correct species
use basic
USE beam
use distrib
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, is, ie, js ,je, mu, nu
REAL(kind=db), ALLOCATABLE:: R(:),Z(:)
REAL(kind=db) :: loc_rlimits(2), loc_zlimits(2)
p=>plist(species(1))
loc_rlimits=rlimits/rnorm
loc_zlimits=zlimits/rnorm
nbionisation=0
allocate(ins_p%start)
ins_p%n=0
created=>ins_p%start
!js=floor((loc_rlimits(1)-rgrid(1))*invdr(1))
!je=floor((loc_rlimits(2)-rgrid(1))*invdr(1))
!is=floor((loc_zlimits(1)-zgrid(1))*invdz(1))
!ie=floor((loc_rlimits(2)-zgrid(1))*invdz(1))
js=1
je=nr+1
is=1
ie=nz+1
DO i=is,ie
DO j=js,je
! Number of particles to add at this time step
mu=(i-1)*(nr+1)+j
nu=i*(nr+femorder(1))+j+femorder(2)-2
IF(ionfreq(mu) .gt. 0) THEN
npartsadd=floor(npartspercell(nu)*ionfreq(mu)*dt*tnorm*itcol)
ELSE
npartsadd=0
END IF
IF (npartsadd.gt. 0) THEN
!WRITE(*,*) npartsadd
nbionisation=nbionisation+npartsadd
loc_rlimits(1)=rgrid(j-1)-dr(1)/2
loc_rlimits(2)=rgrid(j-1)+dr(1)/2
loc_zlimits(1)=zgrid(i-1)-dz(1)/2
loc_zlimits(2)=zgrid(i-1)+dz(1)/2
!WRITE(*,*) loc_rlimits, loc_zlimits
ALLOCATE(R(npartsadd),Z(npartsadd))
! 1/R distribution in R
CALL lodunir(0,R,loc_rlimits(1),loc_rlimits(2))
! Uniform distribution in Z
CALL loduni(0,Z)
! fill the added particles buffer
DO k=1,npartsadd
! prepare the memory for the released electron
ins_p%n=ins_p%n+1
allocate(created%next)
created%next%prev=>created
! fill created particle
created%p%pos=(/R(k), 0.0_db, Z(k)*(loc_zlimits(2)-loc_zlimits(1)) + loc_zlimits(1)/)
created%p%U=0
created%p%gamma=1.0
! We prepare the next created particle
ins_p%end=>created
created=>created%next
END DO
DEALLOCATE(R,Z)
END IF
END DO
END DO
! 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
p%nbionisation=p%nbionisation+nbionisation
if(ins_p%n .gt.0) then
CALL add_created_part(plist(species(1)),ins_p,.true.,.false.)
end if
end subroutine ionisation_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
!sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1)
! Exponential fitting relevant at high energies
!sig_fit=sig_vec(k,2)*(ek/sig_vec(k,1))**growth_vec(k)
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
!sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1)
! 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
end module neutcol
c4science · Help