diff --git a/src/neutcol_mod.f90 b/src/neutcol_mod.f90 index 15a5bf8..362421a 100644 --- a/src/neutcol_mod.f90 +++ b/src/neutcol_mod.f90 @@ -1,440 +1,440 @@ !------------------------------------------------------------------------------ ! 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) 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 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 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. !< NAMELIST /neutcolparams/ neutdens, neuttemp, neutpressure, Eion, & & scatter_fac, nlcol, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, & & nldragio PUBLIC:: neutcol_init, neutcol_step, neutcol_diag 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 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 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) 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(p) ! USE random USE beam USE omp_lib USE basic, ONLY: nlclassical, rnorm USE distrib, ONLY: lodgaus type(particles)::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 - real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3), B + real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3) 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,B), reduction(+:nbcolls_ela,nbcolls_io, nucol) + !$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) 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 ! The ionisation event can only occur if the incoming electron energy is above the binding energy if (ek .gt. Eion) then sig_io=sig_fit(io_cross_sec,io_growth_cross_sec,ek, nb_io_cross) else sig_io=0 end if 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)) 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 IF( nlmaxwellio ) THEN ! the new electron is created according to a Maxwellian CALL lodgaus(0, Rand(1:3)) ! Fill created particle physical quantities created%p%R=p%R(i) created%p%THET=p%THET(i) created%p%Z=p%Z(i) ! 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 B=sqrt(scatter_fac/(Ek-Eion)) - Es=B*tan(Rand(1)*B*atan((Ek-Eion)/(2*scatter_fac))) + Es=B*tan(Rand(1)*B*atan((Ek-Eion)/(2*scatter_fac)))*Ek ! Compute scattering angles for created electron - cosChi=1-2*Rand(2)/(1+8*Ek/E0*(1-Rand(2))) + 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%R=p%R(i) created%p%THET=p%THET(i) created%p%Z=p%Z(i) 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 deallocate(created%next) else deallocate(ins_p(omp_thread)%start) end if !$OMP END PARALLEL j=0 Do i=1,num_threads if(associated(ins_p(i)%start)) then j=i exit end if end do if(j.ge.1) then created=>ins_p(j)%end Do i=j+1,num_threads created%next=>ins_p(i)%start ins_p(j)%n=ins_p(j)%n+ins_p(i)%n IF(ASSOCIATED(created%next)) then created=>ins_p(i)%end END IF End Do CALL add_created_part(p,ins_p(j)) end if 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) real(kind=db)::sig_fit, ek real(kind=db):: sig_vec(:,:), growth_vec(:) Integer:: k, nb_cross sig_fit=0 if(nb_cross .lt. 1 .or. ek .eq. 0) return k=2 DO while(sig_vec(k,1).le. ek .and. k .lt. nb_cross) k=k+1 end do !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-1,2)*(ek/sig_vec(k-1,1))**growth_vec(k-1) 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 U0=(/Ur,Uthet,Uz/) norm=sqrt(sum(U0**2)) U=U0/norm ! Find a vector perpendicular to U for chi rotation perp=(/1,0,0/) IF( U(2) .ne. 0) THEN perp(2)=-U0(1)/U0(2) perp=perp/sqrt(1+perp(2)**2) else if (U(3) .eq. 0)then ! the velocity is purely radial (U(2) and U(3) is 0) perp=(/0,1,0/) END IF ! 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)*U0(2) + perp(2)*sinchi*U0(3) Uthet = perp(1)*perp(2)*(1-coschi)*U0(1) + (coschi + perp(2)**2*(1-coschi))*U0(2) -perp(1)*sinchi*U0(3) Uz = -perp(2)*sinchi*U0(1) +perp(1)*sinchi*U0(2) + coschi*U0(3) U0 =(/Ur,Uthet,Uz/) ! 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