Homec4science
Diffusion FENNECS lite (3D) Browse

fennecslite/src/maxwsrce_mod.f90
6b37315661933D

Last Change
Raw File

maxwsrce_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: maxwellsource
!
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> Adds particle in the simulation according to a maxwellian distribution in velocity
!> and a uniform distribution in space
!------------------------------------------------------------------------------
MODULE maxwsrce
!
USE constants
use mpi
USE mpihelper
USE basic, ONLY: mpirank, mpisize, vnorm, rnorm, zgrid, &
& nlclassical, nlmaxwellsource
USE beam
USE distrib
IMPLICIT NONE
PRIVATE
REAL(kind=db), SAVE :: frequency = 0 !< Number of macro particles added per second over the spawning region
REAL(kind=db), SAVE :: loc_frequency = 0 !< local number of macro particles added per second over the spawning region
REAL(kind=db), SAVE :: temperature=11000 !< temperature used for the Maxwellian velocity distribution in Kelvin
REAL(kind=db), SAVE :: rlimits(2) !< radial limits in which particles will be spawned
REAL(kind=db), SAVE :: zlimits(2) !< axial limits in which particles will be spawned
REAL(kind=db), SAVE :: last_t=0 !< last time when a macro-particle was added to the simulation
REAL(kind=db), SAVE :: vth = 0 !< Normalized thermal velocity computed from temperature
REAL(kind=db), SAVE :: loc_zlimits(2)!< Local axial limits in which particles will be spawned used for current mpi process
REAL(kind=db), SAVE :: loc_rlimits(2)!< Local radial limits in which particles will be spawned used for current mpi process
REAL(kind=db), SAVE :: time_start=-1.0 !< time at which the source is turned on
REAL(kind=db), SAVE :: time_end=-1.0 !< time at which the source is turned off
INTEGER, SAVE :: radialtype=2 !< type of radial distribution used for creating particles
NAMELIST /maxwellsourceparams/ frequency, temperature, rlimits, zlimits, time_start, time_end, radialtype
PUBLIC:: maxwsrce_init, maxwsrce_inject, maxwsrce_diag, maxwsrce_calcfreq
contains
subroutine maxwsrce_init(lu_in, time, Zbounds)
implicit none
INTEGER, INTENT(IN)::lu_in
INTEGER, INTENT(IN):: Zbounds(0:)
REAL(kind=db), INTENT(IN):: time
INTEGER:: ierr
Rewind(lu_in)
READ(lu_in, maxwellsourceparams)
IF(mpirank .eq. 0) THEN
WRITE(*, maxwellsourceparams)
END IF
! compute normalized thermal velocity and source reference time
vth=sqrt(kb*temperature/partslist(1)%m)/vnorm
last_t=time
IF(time_start .gt. last_t) last_t=time_start
time_start=last_t
CALL maxwsrce_calcfreq(Zbounds)
if (loc_frequency .lt. 0) loc_frequency=0
WRITE(*,*) "init local frequency : ", loc_frequency
IF(radialtype .gt.4 .or. radialtype .lt. 1) THEN
IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
End subroutine maxwsrce_init
SUBROUTINE maxwsrce_calcfreq(Zbounds)
IMPLICIT NONE
INTEGER:: Zbounds(0:)
REAL(kind=db):: surface, frequencyrem
REAL(kind=db):: remsurface
frequencyrem=0
! compute source surface
surface=(zlimits(2)-zlimits(1))
loc_rlimits=rlimits/rnorm
! if the source is in the mpi process compute region
if(zlimits(2) .gt. zgrid(Zbounds(mpirank))*rnorm .and. zlimits(1) .lt. zgrid(Zbounds(mpirank+1))*rnorm) then
! reduce the frequency and source boundaries to match the volume covered by the mpi process
IF (zlimits(1) .lt. zgrid(Zbounds(mpirank))*rnorm) THEN
remsurface=(zgrid(Zbounds(mpirank))*rnorm-zlimits(1))
frequencyrem=frequency*remsurface/surface
loc_zlimits(1)=zgrid(Zbounds(mpirank))
ELSE
loc_zlimits(1)=zlimits(1)/rnorm
END IF
IF (zlimits(2) .gt. zgrid(Zbounds(mpirank+1))*rnorm) THEN
remsurface=(zlimits(2)-zgrid(Zbounds(mpirank+1))*rnorm)
frequencyrem=frequencyrem+frequency*remsurface/surface
loc_zlimits(2)=zgrid(Zbounds(mpirank+1))
ELSE
loc_zlimits(2)=zlimits(2)/rnorm
END IF
loc_frequency=frequency-frequencyrem
else
! otherwise turn off the source for this mpi process
loc_frequency=0
loc_zlimits=(/zgrid(Zbounds(mpirank)), zgrid(Zbounds(mpirank+1))/)
end if
END SUBROUTINE maxwsrce_calcfreq
Subroutine maxwsrce_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. nlmaxwellsource) THEN
Write(grpname,'(a,a)') trim(str),"/maxwellsource"
If(.not. isgroup(File_handle, trim(grpname))) THEN
CALL creatg(File_handle, trim(grpname))
END IF
Call attach(File_handle, trim(grpname), "frequency", frequency)
Call attach(File_handle, trim(grpname), "temperature", temperature)
Call putarr(File_handle, trim(grpname)//"/zlimits", zlimits)
Call putarr(File_handle, trim(grpname)//"/rlimits", rlimits)
Call attach(File_handle, trim(grpname), "vth", vth*vnorm)
Call attach(File_handle, trim(grpname), "time_start", time_start)
Call attach(File_handle, trim(grpname), "time_end",time_end)
Call attach(File_handle, trim(grpname), "radialtype", radialtype)
END IF
End subroutine maxwsrce_diag
subroutine maxwsrce_inject(time)
! Inject the particles in the correct species according to the source definition
Use omp_lib
use basic,ONLY: mpirank
implicit none
REAL(kind=db), INTENT(IN) :: time
Type(particle), ALLOCATABLE, Dimension(:):: newparts
INTEGER:: npartsadd
INTEGER:: i
REAL(kind=db), ALLOCATABLE:: R(:),Z(:)
! check if source is on
IF(.not. maxwsrce_on(time)) THEN
RETURN
END IF
! Number of particles to add at this time step
npartsadd=floor((time-last_t)*loc_frequency)
IF (npartsadd.gt. 0) THEN
ALLOCATE(newparts(npartsadd))
ALLOCATE(R(npartsadd),Z(npartsadd))
SELECT CASE(radialtype)
CASE(1) ! 1/R distribution in R
CALL lodunir(0,R,loc_rlimits(1),loc_rlimits(2))
CASE(2) ! flat top distribution in R
CALL lodlinr(0,R,loc_rlimits(1),loc_rlimits(2))
CASE(3) ! 1/R^2 distribution in R
CALL lodinvr(0,R,loc_rlimits(1),loc_rlimits(2))
CASE(4) ! gaussian distribution in R
CALL lodgausr(0,R,loc_rlimits(1),loc_rlimits(2))
CASE DEFAULT
IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype
END SELECT
CALL loduni (0,Z)
! fill the added particles buffer
DO i=1,npartsadd
newparts(i)%U=0
newparts(i)%pos=(/R(i), 0.0_db, Z(i)*(loc_zlimits(2)-loc_zlimits(1)) + loc_zlimits(1)/)
newparts(i)%gamma=1.0
END DO
! Move the buffer of created particles to the simulation buffer
CALL add_created_part(partslist(1), newparts)
last_t=last_t+npartsadd/loc_frequency
END IF
end subroutine maxwsrce_inject
logical function maxwsrce_on(time)
REAL(kind=db), intent(in):: time
maxwsrce_on=.false.
IF (time_start .lt. 0 .and. time_end .lt. 0) THEN
maxwsrce_on = .true.
ELSE IF (time .gt. time_start .and. (time .lt. time_end .or. time_end .lt. 0) ) THEN
maxwsrce_on = .true.
END IF
maxwsrce_on=nlmaxwellsource .and. maxwsrce_on
end function
END MODULE maxwsrce
c4science · Help