!------------------------------------------------------------------------------ ! 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 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 !< 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) :: 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_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 [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 NAMELIST /neutcolparams/ neutdens, Eion, & & scatter_fac, nlcol, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, & & nldragio, itcol, species, isotropic PUBLIC:: neutcol_init, neutcol_step, neutcol_diag, itcol, neutdens PROCEDURE(rotate_vel), POINTER:: change_dir => NULL()!< Function evaluating the weight for Dirichelt boundary conditions ABSTRACT INTERFACE SUBROUTINE rotate_vel(Ur, Uthet, Uz, coschi, thet) use constants real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet END SUBROUTINE end interface 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 if (.not. isotropic)then change_dir=> rotate else change_dir=> scatter end if 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 !------------------------------------------------------------- !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Simulates the elastic and ionising collisions for each particles in plist(species(1)) ! !> @param [inout] plist list of particle species considered in the code !--------------------------------------------------------------------------- SUBROUTINE neutcol_step(plist) ! USE random USE beam USE omp_lib USE basic, ONLY: nlclassical 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):: ins_p type(linked_part), POINTER:: created 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() nbcolls_ela=0 nbcolls_io=0 nucol=0 !!$OMP 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 !omp_thread=1 allocate(ins_p%start) ins_p%n=0 created=>ins_p%start !$OMP DO DO i=1,p%Nploc !for each particle CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread)) ! we calculate the kinetic energy and norm of the velocity v2=(p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2) if(nlclassical) THEN ek=v2*escale v=sqrt(v2) vinit=p%U(:,i) ! (/p%UR(i),p%UTHET(i),p%UZ(i)/) ELSE ek=(p%gamma(i)-1)*escale v=sqrt(v2)/p%gamma(i) vinit=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i) end if sig_io=0 sig_ela=0 ! computes the ionisation and elastic collision cross-sections at this kinetic energy ! 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) end if collisionfact=1-exp(-collfactor*(sig_io+sig_ela)*v) ! If we have a collision event if (Rand(1) .lt.collisionfact) THEN CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread)) ! Check if elastic or ionising event is happening IF(Rand(1).gt. sig_ela/(sig_io+sig_ela)) THEN ! An ionisation collision happened and we create the necessary electron ! prepare the memory for the released electron ins_p%n=ins_p%n+1 allocate(created%next) created%next%prev=>created ! Fill created particle new position created%p%pos=p%pos(:,i)!(/p%R(i), p%THET(i), p%Z(i)/) IF( nlmaxwellio ) THEN ! the new electron velocity is defined according to a Maxwellian CALL lodgaus(0, Rand(1:3)) ! get random velocity created%p%U=vth*Rand(1: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 if (isotropic) then coschi=cos(Rand(2)*pi) else cosChi=1-2*Rand(2)/(1+8*Es/E0*(1-Rand(2))) end if 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 velocity created%p%U=vfact*p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/) ! rotate the velocity vector due to the collision call change_dir(created%p%U(1),created%p%U(2), created%p%U(3), coschi, thet) END IF vend=created%p%U 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%U(1)**2+created%p%U(2)**2+created%p%U(3)**2) vend=vend/created%p%gamma END IF ! We prepare the next created particle ins_p%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 ! its new kinetic energy 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 if (isotropic) then coschi=cos(Rand(1)*pi) else cosChi=1-2*Rand(1)/(1+8*Es/E0*(1-Rand(1))) end if thet=Rand(2)*2*pi ! Change the incident electron velocity direction and amplitude if necessary p%U(:,i)=p%U(:,i)*vfact !p%UTHET(i)=p%UTHET(i)*vfact !p%UZ(i)=p%UZ(i)*vfact call change_dir(p%U(1,i),p%U(2,i), p%U(3,i), coschi, thet) if(nlclassical) THEN vend=p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/) ELSE p%gamma(i)=sqrt(1+p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2) vend=p%U(:,i)/p%gamma(i)!(/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 !$OMP BARRIER ! 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 END PARALLEL ! We collect all created particules into one linked list for easier insertion in plist !Do i=1,num_threads !$OMP CRITICAL (insertions) 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 (insertions) !$OMP CRITICAL (insertelectrons) if(ins_p%n .gt.0) then CALL add_created_part(plist(species(1)),ins_p,.true.,.false.) !exit end if !$OMP END CRITICAL (insertelectrons) !end do !$OMP CRITICAL(addcolls) p%nbcolls=p%nbcolls+(/nbcolls_io, nbcolls_ela/) p%nudcol=p%nudcol+nucol !$OMP END CRITICAL(addcolls) !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) 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 scatter(Ur, Uthet, Uz, coschi, thet) real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet real(kind=db):: norm real(kind=db):: sinchi, sinthet, costhet norm=sqrt(Ur**2+Uz**2+Uthet**2) ! Compute sinus and cosinus for rotation sinchi=sqrt(1-coschi**2) costhet=cos(thet) sinthet=sin(thet) Ur=norm*sinchi*costhet Uthet=norm*sinchi*sinthet Uz=norm*coschi END SUBROUTINE scatter 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 end module neutcol