Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (geom_testmumps) Browse

fennecs/src/neutcol_mod.f90
e4669700ed0dgeom_testmumps

Last Change
Raw File

neutcol_mod.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 :: 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 each evaluation of neutcol_step
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=1e-6 !< Neutral particle pressure in mbar
Real(kind=db) :: scatter_fac = 24.2 !< Energy scattering factor for the considered gas (here for Ne)
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
real(kind=db) :: collfactor !< collision factor (n_n \sigma \delta t)
INTEGER :: nb_io_cross=0
Real(kind=db), ALLOCATABLE :: io_cross_sec(:,:) !< Ionisation cross-section table
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_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
Real(kind=db) :: vth !< In case of nlmaxwelio, defines the thermal velocity of created electrons
LOGICAL :: nldragio=.true. !<
INTEGER :: species(2)
NAMELIST /neutcolparams/ neutdens, neuttemp, neutpressure, Eion, &
& scatter_fac, nlcol, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, &
& nldragio, itcol, species
PUBLIC:: neutcol_init, neutcol_step, neutcol_diag, itcol
CONTAINS
subroutine neutcol_init(lu_in, time, p)
Use mpi
Use basic, only: mpirank, dt, nlclassical,rnorm,tnorm, vnorm
Use beam, only: particles
Use constants
implicit none
INTEGER, INTENT(IN) :: lu_in
REAL(kind=db), INTENT(IN):: time
TYPE(particles) :: p
INTEGER:: ierr, istat
character(len=1000) :: line
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
! 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/tnorm*rnorm**3*itcol
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
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)
if (allocated(ela_cross_sec)) Call putarr(File_handle, trim(grpname)//"/ela_cross_sec", ela_cross_sec)
END IF
End subroutine neutcol_diag
!-------------------------------------------------------------
SUBROUTINE neutcol_step(plist)
!
USE random
USE beam
USE omp_lib
USE basic, ONLY: nlclassical, rnorm
USE distrib, ONLY: lodgaus
type(particles), TARGET::plist(:)
type(particles),pointer::p
INTEGER:: i, omp_thread, num_threads, j, nbcolls_ela, nbcolls_io
real(kind=db):: Rand(5)
real(kind=db):: v2, v, ek, Everif, es, cosChi, thet, sig_io, sig_ela, vfact, xsi
type(linked_part_row), ALLOCATABLE:: ins_p(:)
type(linked_part), POINTER:: created, created2
real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3)
p=>plist(species(1))
if(.not. nlcol .or. p%nploc .le. 0) return
num_threads=omp_get_max_threads()
Allocate(ins_p(num_threads))
nbcolls_ela=0
nbcolls_io=0
nucol=0
!$OMP PARALLEL DEFAULT(SHARED), private(collisionfact,i,omp_thread,Rand,v2,ek,sig_io,sig_ela,es,coschi,thet,vfact, created, v, everif,xsi,vinit,vend), reduction(+:nbcolls_ela,nbcolls_io, nucol)
omp_thread=omp_get_thread_num()+1
allocate(ins_p(omp_thread)%start)
ins_p(omp_thread)%n=0
created=>ins_p(omp_thread)%start
!$OMP DO
DO i=1,p%Nploc
CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread))
v2=(p%UR(i)**2+p%UTHET(i)**2+p%UZ(i)**2)
if(nlclassical) THEN
ek=v2*escale
v=sqrt(v2)
vinit=(/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
ek=(p%gamma(i)-1)*escale
v=sqrt(v2)/p%gamma(i)
vinit=(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
end if
sig_io=0
sig_ela=0
! The ionisation event can only occur if the incoming electron energy is above the binding 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)
xsi=Ek/(0.25*E0+Ek)
sig_ela=sig_ela*(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi))
end if
collisionfact=1-exp(-collfactor*(sig_io+sig_ela)*v)
! If we have a collision
if (Rand(1) .lt.collisionfact) THEN
CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread))
IF(Rand(1).gt. sig_ela/(sig_io+sig_ela)) THEN ! An ionisation collision happened and we create the necessary electron
ins_p(omp_thread)%n=ins_p(omp_thread)%n+1
allocate(created%next)
created%next%prev=>created
! Fill created particle physical quantities
created%p%R=p%R(i)
created%p%THET=p%THET(i)
created%p%Z=p%Z(i)
IF( nlmaxwellio ) THEN ! the new electron is created according to a Maxwellian
CALL lodgaus(0, Rand(1:3))
! get random velocity
created%p%UR=vth*(Rand(1))
created%p%UTHET=vth*(Rand(2))
created%p%UZ=vth*(Rand(3))
ELSE
CALL random_array(Rand,3,ran_index(omp_thread),ran_array(:,omp_thread))
! Compute created electron energy
Es=scatter_fac*tan(Rand(1)*atan((Ek-Eion)/(2*scatter_fac)))
! Compute scattering angles for created electron
cosChi=1-2*Rand(2)/(1+8*Es/E0*(1-Rand(2)))
thet=Rand(3)*2*pi
if(nlclassical)THEN
! new velocity factor for created particle
vfact=sqrt(Es/Ek)
ELSE
! new velocity factor for created particle
vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale)))
END IF
! Fill created particle physical quantities
created%p%UR=vfact*(p%UR(i))
created%p%UTHET=vfact*(p%UTHET(i))
created%p%UZ=vfact*(p%UZ(i))
call rotate(created%p%UR,created%p%UTHET, created%p%UZ, coschi, thet)
END IF
vend=(/created%p%UR,created%p%UTHET,created%p%UZ/)
if(nlclassical)THEN
! Lorentz factor for created particle
created%p%gamma=1.0
ELSE
! Lorentz factor for created particle
created%p%gamma=sqrt(1+created%p%UZ**2+created%p%UR**2+created%p%UTHET**2)
vend=vend/created%p%gamma
END IF
! We prepare the next created particle
ins_p(omp_thread)%end=>created
created=>created%next
! We keep track of what changed
nbcolls_io=nbcolls_io+1
nucol=nucol-vend/vinit
! If we want the incoming electron to be scattered, we need to compute the
if (nldragio) THEN
! We store the lossed energy in pot for keeping track of energy conservation
created%prev%p%pot=Eion+Es
CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread))
Es=Ek-Eion-Es
if(nlclassical)THEN
! new velocity factor for scattered particle
vfact=sqrt(Es/Ek)
ELSE
! new velocity factor for scattered particle
vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale)))
END IF
ELSE
CYCLE
END IF
ELSE ! An elastic collision event happens
CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread))
Es=Ek
vfact=1.0
nbcolls_ela=nbcolls_ela+1
END IF
! We calculate the scattered velocity angle for the scattered electron
cosChi=1-2*Rand(1)/(1+8*Es/E0*(1-Rand(1)))
thet=Rand(2)*2*pi
! Change the incident electron velocity direction and amplitude if necessary
p%UR(i)=p%UR(i)*vfact
p%UTHET(i)=p%UTHET(i)*vfact
p%UZ(i)=p%UZ(i)*vfact
call rotate(p%UR(i),p%UTHET(i), p%UZ(i), coschi, thet)
if(nlclassical) THEN
vend=(/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
p%gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2)
vend=(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
END IF
nucol=nucol+1-vend/vinit
END IF
END DO
!$OMP END DO
if(associated(created%prev)) then
created=>created%prev
ins_p(omp_thread)%end=>created
deallocate(created%next)
else
deallocate(ins_p(omp_thread)%start)
end if
!$OMP END PARALLEL
! We collect all created particules into one linked list
Do i=1,num_threads
if(associated(ins_p(i)%start)) then
created=>ins_p(i)%end
Do j=i+1,num_threads
created%next=>ins_p(j)%start
ins_p(i)%n=ins_p(i)%n+ins_p(j)%n
IF(ASSOCIATED(created%next)) then
created=>ins_p(j)%end
END IF
End Do
if(species(2).gt.0) then
CALL add_created_part(plist(species(2)), ins_p(i),.false.,.true.)
end if
CALL add_created_part(p,ins_p(i),.true.,.false.)
exit
end if
end do
DEALLOCATE(ins_p)
p%nbcolls=p%nbcolls+(/nbcolls_io, nbcolls_ela/)
p%nudcol=nucol
!Write(*,*)"mpirank: ", mpirank, " Nb colls ela, io: ",nbcolls_ela, nbcolls_io
!
END SUBROUTINE neutcol_step
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)
END FUNCTION sig_fit
SUBROUTINE rotate(Ur, Uthet, Uz, coschi, thet)
real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet
real(kind=db):: norm, perp(3), U(3), U0(3), normf
real(kind=db):: sinchi, sinthet, costhet
Integer :: iperp1,iperp2
U0=(/Ur,Uthet,Uz/)
norm=sqrt(sum(U0**2))
U=U0/norm
! Find a vector perpendicular to U for chi rotation
! find the direction with maximum amplitude
perp=(/1,1,1/)
iperp1=maxloc(abs(U),1)
! find second direction with next max amplitude
perp(iperp1)=0
iperp2=maxloc(abs(perp*U),1)
perp=0
perp(iperp2)=U(iperp1)
perp(iperp1)=-U(iperp2)
! Normalise the rotation vector
perp=perp/sqrt(sum(perp**2))
! Compute sinus and cosinus for rotation
sinchi=sqrt(1-coschi**2)
costhet=cos(thet)
sinthet=sin(thet)
! Rotation of angle chi around perp
Ur = (coschi+perp(1)**2*(1-coschi))*U0(1) + (perp(1)*perp(2)*(1-coschi)-perp(3)*sinchi)*U0(2) + (perp(1)*perp(3)*(1-coschi) + perp(2)*sinchi)*U0(3)
Uthet = (perp(1)*perp(2)*(1-coschi)+perp(3)*sinchi)*U0(1) + (coschi + perp(2)**2*(1-coschi))*U0(2) +(perp(2)*perp(3)*(1-coschi)-perp(1)*sinchi)*U0(3)
Uz = (perp(1)*perp(3)*(1-coschi)-perp(2)*sinchi)*U0(1) +(perp(3)*perp(2)*(1-coschi)+perp(1)*sinchi)*U0(2) +( coschi+perp(3)**2*(1-coschi))*U0(3)
U0 =(/Ur,Uthet,Uz/)
! second rotation according to uniform distribution
! Rotation of angle theta around U
Ur = (costhet+U(1)**2*(1-costhet))*U0(1) + (U(1)*U(2)*(1-costhet) - U(3)*sinthet)*U0(2) + (U(1)*U(3)*(1-costhet)+U(2)*sinthet)*U0(3)
Uthet = (U(2)*U(1)*(1-costhet)+U(3)*sinthet)*U0(1) + (costhet + U(2)**2*(1-costhet))*U0(2) + (U(2)*U(3)*(1-costhet)-U(1)*sinthet)*U0(3)
Uz = (U(3)*U(1)*(1-costhet) - U(2)*sinthet)*U0(1) + (U(3)*U(2)*(1-costhet)+U(1)*sinthet)*U0(2) + (costhet +U(3)**2*(1-costhet))*U0(3)
!normf=sqrt(Ur**2+Uthet**2+Uz**2)
!if(abs(norm-normf)/norm .gt. 1e-14) WRITE(*,*) "Error in rotate the norm of v changed"
END SUBROUTINE rotate
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, j
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
end module neutcol
c4science · Help