!------------------------------------------------------------------------------ ! 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) 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. 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(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) ! fill the added particles buffer 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=newparts(i)%UR*newparts(i)%gamma newparts(i)%UZ=newparts(i)%UZ*newparts(i)%gamma newparts(i)%UTHET=newparts(i)%UTHET*newparts(i)%gamma END IF 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