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

fennecs/src/maxwellsource_mod.f90
a2847b9464fageom_testpardiso

Last Change
Raw File

maxwellsource_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 maxwellsource
!
USE constants
USE mpi
USE mpihelper
USE basic, ONLY: mpirank, mpisize, vnorm, rnorm, Zbounds, zgrid, &
& nlclassical, nlmaxwellsource
USE beam
USE distrib
IMPLICIT NONE
PRIVATE
REAL(kind=db), SAVE :: frequency !< Number of macro particles added per second over the spawning region
REAL(kind=db), SAVE :: temperature !< temperature used for the Maxwellian velocity distribution
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 !< 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:: maxwellsource_init, maxwellsource_inject
contains
subroutine maxwellsource_init(lu_in, time)
implicit none
INTEGER, INTENT(IN)::lu_in
REAL(kind=db), INTENT(IN):: time
REAL(kind=db):: surface
REAL(kind=db):: remsurface
INTEGER:: ierr
READ(lu_in, maxwellsourceparams)
IF(mpirank .eq. 0) THEN
WRITE(*, maxwellsourceparams)
END IF
vth=sqrt(kb*temperature/partslist(1)%m)/vnorm !thermal velocity
last_t=time
IF(time_start .gt. last_t) last_t=time_start
surface=(rlimits(2)-rlimits(1))*(zlimits(2)-zlimits(1))
loc_rlimits=rlimits/rnorm
IF (zlimits(1) .lt. zgrid(Zbounds(mpirank))*rnorm) THEN
remsurface=(rlimits(2)-rlimits(1))*(zgrid(Zbounds(mpirank))*rnorm-zlimits(1))
frequency=frequency*(1-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=(rlimits(2)-rlimits(1))*(zlimits(2)-zgrid(Zbounds(mpirank+1))*rnorm)
frequency=frequency*(1-remsurface/surface)
loc_zlimits(2)=zgrid(Zbounds(mpirank+1))
ELSE
loc_zlimits(2)=zlimits(2)/rnorm
END IF
WRITE(*,*) "init frequency : ", 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 maxwellsource_init
subroutine maxwellsource_inject(time)
implicit none
REAL(kind=db), INTENT(IN) :: time
Type(particle), ALLOCATABLE, Dimension(:):: newparts
INTEGER:: npartsadd
INTEGER:: i
REAL(kind=db), ALLOCATABLE:: VR(:),VZ(:),VTHET(:),R(:),Z(:)
! check if source is on
IF(.not. maxwellsource_on(time)) THEN
RETURN
END IF
! Number of particles to add at this time step
npartsadd=floor((time-last_t)*frequency)
IF (npartsadd.gt. 0) THEN
ALLOCATE(newparts(npartsadd))
ALLOCATE(VR(npartsadd),VZ(npartsadd),VTHET(npartsadd),R(npartsadd),Z(npartsadd))
! Initial velocities distribution
CALL lodgaus(0,VZ)
CALL lodgaus(0,VR)
CALL lodgaus(0,VTHET)
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)
DO i=1,npartsadd
newparts(i)%UR=VR(i)*vth
newparts(i)%UZ=VZ(i)*vth
newparts(i)%UTHET=VTHET(i)*vth
newparts(i)%R=R(i)
newparts(i)%Z=Z(i)*(loc_zlimits(2)-loc_zlimits(1)) + loc_zlimits(1)
newparts(i)%THET=0
IF (nlclassical) THEN
newparts(i)%gamma=1.0
ELSE
newparts(i)%gamma=sqrt(1/(1-newparts(i)%UZ**2-newparts(i)%UR**2-newparts(i)%UTHET**2))
newparts(i)%UR=VR(i)*newparts(i)%gamma
newparts(i)%UZ=VZ(i)*newparts(i)%gamma
newparts(i)%UTHET=VTHET(i)*newparts(i)%gamma
END IF
END DO
CALL add_created_part(partslist(1), newparts)
last_t=last_t+npartsadd/frequency
END IF
end subroutine maxwellsource_inject
logical function maxwellsource_on(time)
REAL(kind=db), intent(in):: time
maxwellsource_on=.false.
IF (time_start .lt. 0 .and. time_end .lt. 0) THEN
maxwellsource_on = .true.
ELSE IF (time .gt. time_start .and. (time .lt. time_end .or. time_end .lt. 0) ) THEN
maxwellsource_on = .true.
END IF
maxwellsource_on=nlmaxwellsource .and. maxwellsource_on
end function
End Module maxwellsource
c4science · Help