Homec4science
Diffusion FENNECS lite (master) Browse

fennecslite/src/beam_mod.f90
27cb867868c1master

Last Change
Raw File

beam_mod.f90
View Options

MODULE beam
!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: beam
!
!> @author
!> Guillaume Le Bars EPFL/SPC
!> Patryk Kaminski EPFL/SPC
!> Trach Minh Tran EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for loading, advancing and computing the necessary diagnostics for the simulated particles.
!------------------------------------------------------------------------------
!
USE constants
use mpi
USE mpihelper
USE basic, ONLY: mpirank, mpisize
USE distrib
USE particletypes
USE weighttypes
IMPLICIT NONE
!
!TYPE(particles) :: parts !< Storage for all the particles
!SAVE :: parts
TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist
TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist_towrite
! Diagnostics (scalars)
REAL(kind=db) :: ekin=0 !< Total kinetic energy (J)
REAL(kind=db) :: epot=0 !< Total potential energy (J)
REAL(kind=db) :: etot=0 !< Current total energy (J)
REAL(kind=db) :: etot0=0 !< Initial total energy (J)
REAL(kind=db) :: loc_etot0=0 !< theoretical local total energy (J)
REAL(kind=db) :: Energies(4) !< (1) kinetic energy, (2) potential energy, (3) total energy and (4) gained/lossed energy due to gain or loss of particles (J)
!
INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process
INTERFACE add_created_part
MODULE PROCEDURE add_linked_created_part, add_list_created_part
END INTERFACE add_created_part
!
abstract interface
subroutine rloader(nbase,y,rminus,rplus)
USE constants
REAL(kind=db), INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
REAL(kind=db), INTENT(in) :: rplus, rminus
end subroutine
REAL(kind=db) FUNCTION gamma(UZ, UR, UTHET)
USE constants
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
end FUNCTION
end interface
CONTAINS
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Loads the particles at the beginning of the simulation and create the parts variable if necessary
!---------------------------------------------------------------------------
SUBROUTINE load_parts
USE basic, ONLY: nplasma, mpirank, ierr, distribtype, nlclassical, nbspecies, partfile
use mpi
INTEGER:: i
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma))
! Select case to define the type of distribution
SELECT CASE(distribtype)
CASE(1) ! Gaussian distribution in V, uniform in Z and 1/R in R
CALL loaduniformRZ(partslist(1), VR, VZ, VTHET)
CASE(2) !Stable distribution from Davidson 4.95 p.119
CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir)
CASE(3) !Stable distribution from Davidson 4.95 p.119 but with constant distribution in R
CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodinvr)
CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R
CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodgausr)
CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp
CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir)
CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz<V_perp to satisfy magnetic mirror trapping
CALL loaduniformRZ(partslist(1), VR, VZ, VTHET)
VZ = VZ/50
CASE(7) ! Distribution defined in separate input file
CALL read_part_file(partslist(1), partfile(1), VR, VZ, VTHET)
nplasma=partslist(1)%Nptot
CASE DEFAULT
IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END SELECT
Do i=1,nplasma
partslist(1)%partindex(i)=i
END DO
partslist(1)%newindex=nplasma
partslist(1)%U=0
partslist(1)%Upar=0
partslist(1)%larmorrad=0
partslist(1)%E=0
partslist(1)%B=1
IF(nlclassical) THEN
partslist(1)%Gamma(1:nplasma)=1.0_db
ELSE
partslist(1)%Gamma(1:nplasma)=sqrt(1/(1-VR**2-VZ**2-VTHET**2))
END IF
! Normalization of the velocities
partslist(1)%U(1,1:nplasma)=partslist(1)%Gamma(1:nplasma)*VR
partslist(1)%U(3,1:nplasma)=partslist(1)%Gamma(1:nplasma)*VZ
partslist(1)%U(2,1:nplasma)=partslist(1)%Gamma(1:nplasma)*VTHET
DEALLOCATE(VZ, VR, VTHET)
CALL boundary_loss(partslist(1))
partslist(1)%Nptot=partslist(1)%nploc
DO i=2,nbspecies
call load_part_file(partslist(i),partfile(i))
END DO
Do i=1,nbspecies
call creat_parts(partslist_towrite(i),size(partslist(i)%pot,1))
end do
partslist(1)%calc_moments=.true.
partslist(1)%is_field=.true.
END SUBROUTINE load_parts
SUBROUTINE load_part_file(p,partfilename)
use mpi
USE basic, ONLY: nlclassical
type(particles) :: p
CHARACTER(len=*):: partfilename
Real(kind=db), ALLOCATABLE:: VR(:), VZ(:), VTHET(:)
INTEGER j
! Read the actual file and load the position in p
CALL read_part_file(p, partfilename, VR, VZ, VTHET)
Do j=1,p%Nploc
p%partindex(j)=j
END DO
IF(nlclassical) THEN
p%Gamma=1.0_db
ELSE
p%Gamma(1:p%Nptot)=sqrt(1/(1-VR**2-VZ**2-VTHET**2))
END IF
! Normalization of the velocities
p%U(1,1:p%Nptot) = p%Gamma(1:p%Nptot)*VR(1:p%Nptot)
p%U(3,1:p%Nptot) = p%Gamma(1:p%Nptot)*VZ(1:p%Nptot)
p%U(2,1:p%Nptot)= p%Gamma(1:p%Nptot)*VTHET(1:p%Nptot)
DEALLOCATE(VZ, VR, VTHET)
call boundary_loss(p)
END SUBROUTINE load_part_file
!---------------------------------------------------------------------------
! DESCRIPTION:
!> @brief Checks for each particle if the z position is outside of the local/global simulation space.
!> Depending on the boundary conditions, the leaving particles are sent to the correct neighbouring MPI process
!> or deleted.
!
!> @param[in] p particles structure
!
!> @author Guillaume Le Bars EPFL/SPC
!---------------------------------------------------------------------------
SUBROUTINE bound(p)
USE basic, ONLY: zgrid, nz, rgrid, nr, Zbounds, mpirank, step, leftproc, rightproc, partperiodic
use omp_lib
IMPLICIT NONE
type(particles), INTENT(INOUT):: p
INTEGER :: i,j, rsendnbparts, lsendnbparts, nblostparts
INTEGER :: receivednbparts, partdiff, commid(2)
LOGICAL:: leftcomm, rightcomm
INTEGER, ALLOCATABLE:: partstoremove(:)
INTEGER,allocatable :: nblost(:)
allocate(nblost(size(p%nblost,1)))
nblost=0
IF (p%Nploc .gt. 0) THEN
! We communicate with the left processus
leftcomm = leftproc .ne. -1
! We communicate with the right processus
rightcomm = rightproc .ne. -1
commid=(/leftproc,rightproc/)
! Boundary condition at z direction
!$OMP DO SIMD
DO i=1,p%Nploc
p%losthole(i)=0
p%sendhole(i)=0
! If the particle is to the right of the local simulation space, it is sent to the right MPI process
IF (p%pos(3,i) .ge. zgrid(Zbounds(mpirank+1))) THEN
IF(partperiodic) THEN
DO WHILE (p%pos(3,i) .GT. zgrid(nz))
p%pos(3,i) = p%pos(3,i) - zgrid(nz) + zgrid(0)
END DO
END IF
IF(rightcomm) THEN
p%sendhole(i)=i
ELSE IF(.not. partperiodic) THEN
p%losthole(i)=2
END IF
! If the particle is to the left of the local simulation space, it is sent to the left MPI process
ELSE IF (p%pos(3,i) .lt. zgrid(Zbounds(mpirank))) THEN
IF(partperiodic) THEN
DO WHILE (p%pos(3,i) .LT. zgrid(0))
p%pos(3,i) = p%pos(3,i) + zgrid(nz) - zgrid(0)
END DO
END IF
IF(leftcomm) THEN
! We send the particle to the left process
p%sendhole(i)=-i
ELSE IF(.not. partperiodic) THEN
! we destroy the particle
p%losthole(i)=1
END IF
END IF
END DO
!$OMP END DO SIMD
END IF
!$OMP MASTER
receivednbparts=0
nblost=0
j=1
rsendnbparts=0
lsendnbparts=0
Do i=1,p%Nploc
if(p%sendhole(i) .eq. 0) cycle
p%sendhole(j)=p%sendhole(i)
if(p%sendhole(i).gt.0)then
rsendnbparts=rsendnbparts+1
else
lsendnbparts=lsendnbparts+1
end if
j=j+1
end do
j=1
nblostparts=0
Do i=1,p%Nploc
if(p%losthole(i) .eq. 0) cycle
nblost(p%losthole(i))=nblost(p%losthole(i))+1
p%losthole(j)=i
j=j+1
nblostparts=nblostparts+1
end do
p%nblost=nblost+p%nblost
IF(mpisize .gt. 1) THEN
! We send the particles leaving the local simulation space to the closest neighbour
CALL particlescommunication(p, lsendnbparts, rsendnbparts, receivednbparts, commid)
END IF
! If the boundary conditions are not periodic, we delete the corresponding particles
IF(nblostparts .gt. 0 .and. step .ne. 0) THEN
DO i=1,nblostparts
CALL delete_part(p, p%losthole(i), .false. )
END DO
!WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank
END IF
! computes if we received less particles than we sent
partdiff=max(lsendnbparts+rsendnbparts-receivednbparts,0)
IF(nblostparts + partdiff .gt. 0) THEN
ALLOCATE(partstoremove(nblostparts+partdiff))
partstoremove(1:partdiff)=abs(p%sendhole(receivednbparts+1:receivednbparts+partdiff))
partstoremove(partdiff+1:partdiff+nblostparts)=abs(p%losthole(1:nblostparts))
call LSDRADIXSORT(partstoremove,nblostparts + partdiff)
! If we received less particles than we sent, or lost particles we fill the remaining holes with the particles from the end of the parts arrays
DO i=nblostparts+partdiff,1,-1
CALL move_part(p, p%Nploc, partstoremove(i))
p%partindex(p%Nploc)=-1
p%Nploc = p%Nploc-1
END DO
deallocate(partstoremove)
END IF
!$OMP END MASTER
deallocate(nblost)
END subroutine bound
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Check if a particle is outside the simulation domain and remove it if needed
!> @param[in] p particles structure
!---------------------------------------------------------------------------
SUBROUTINE boundary_loss(p)
USE basic, ONLY: rgrid, nr
Use geometry, ONLY: geom_weight, dom_weight,is_insidegeom
Use omp_lib
USE see
USE iiee
IMPLICIT NONE
type(particles), INTENT(INOUT):: p
INTEGER :: i,j,isup, nblostparts, iend,nbunch,rindex,zindex
INTEGER, DIMENSION(16)::idwall
INTEGER :: nblost(size(p%nblost,1)), ii, Nploc_init, Nploc_new
INTEGER :: index(3)
REAL(kind=db) :: pos(3), B
logical::inside
nblost=0
nbunch=16
IF (p%Nploc .le. 0) return
!$OMP DO
DO i=1,p%Nploc,nbunch
! Avoid segmentation fault caused by accessing non relevant data
iend=min(i+nbunch-1,p%Nploc)
p%losthole(i:iend)=0
! calculate the weight do determine if a particle is inside the simulation domain.
do j=i,iend
! If the particle is above or below the simulation domain
IF(p%pos(1,j) .gt. rgrid(nr)) THEN
p%losthole(j)=4
cycle
else if(p%pos(1,j) .lt. rgrid(0))then
p%losthole(j)=3
cycle
end if
! If the particle is to the right of the local simulation space, it is sent to the right MPI process
IF (p%pos(3,j) .ge. zgrid(nz)) THEN
IF(partperiodic) THEN
DO WHILE (p%pos(3,j) .GT. zgrid(nz))
p%pos(3,j) = p%pos(3,j) - zgrid(nz) + zgrid(0)
END DO
ELSE
p%losthole(j)=2
cycle
END IF
! If the particle is to the left of the local simulation space, it is sent to the left MPI process
ELSE IF (p%pos(3,j) .lt. zgrid(0)) THEN
IF(partperiodic) THEN
DO WHILE (p%pos(3,j) .LT. zgrid(0))
p%pos(3,j) = p%pos(3,j) + zgrid(nz) - zgrid(0)
END DO
ELSE
! we destroy the particle
p%losthole(j)=1
cycle
END IF
END IF
call p_calc_cellindex(p,j)
! call is_insidegeom(p%pos(1,j), p%pos(3,j),p%cellindex(1,j),p%cellindex(3,j),idwall(j-i+1),inside)
! if(.not. inside) then
! ! If the particle is outside of the vacuum region it is deleted.
! p%losthole(j)=max(idwall(j-i+1),1)
! end if
B=sqrt(p%B(1,j)**2+p%B(2,j)**2)
pos=(/ p%pos(1,j)+p%larmorrad(j)*p%B(2,j)/B,p%pos(2,j),p%pos(3,j)-p%larmorrad(j)*p%B(1,j)/B /)
call calc_cellindex(index,pos)
call is_insidegeom(pos(1),pos(3),index(1),index(3),idwall(j-i+1),inside)
if(.not. inside) then
! If the particle is outside of the vacuum region it is deleted.
p%losthole(j)=max(idwall(j-i+1),1)
end if
pos=(/ p%pos(1,j)-p%larmorrad(j)*p%B(2,j)/B,p%pos(2,j),p%pos(3,j)+p%larmorrad(j)*p%B(1,j)/B /)
call calc_cellindex(index,pos)
call is_insidegeom(pos(1),pos(3),index(1),index(3),idwall(j-i+1),inside)
if(.not. inside) then
! If the particle is outside of the vacuum region it is deleted.
p%losthole(j)=max(idwall(j-i+1),1)
end if
end do
call geom_weight(p%pos(1,i:iend), p%pos(3,i:iend), p%geomweight(:,i:iend))
END DO
!$OMP END DO
!$OMP MASTER
nblostparts=0
nblost=0
j=1
do i=1,p%Nploc
if(p%losthole(i) .eq. 0) cycle
nblost(p%losthole(i)+4)=nblost(p%losthole(i)+4)+1
p%losthole(j)=i
j=j+1
nblostparts=nblostparts+1
end do
p%nblost=nblost+p%nblost
IF(nblostparts.gt.0) THEN
IF(p%see_id.gt.0) THEN ! Secondary electron emission
Nploc_init = p%Nploc
CALL secondary_electrons(p, p%losthole, nblostparts)
Nploc_new = p%Nploc
if (Nploc_new-Nploc_init .ge. 1) then
DO ii =Nploc_init+1,Nploc_new
Call p_calc_cellindex(p,ii)
call geom_weight(p%pos(1,ii), p%pos(3,ii), p%geomweight(:,ii))
END DO
end if
ELSE IF(p%iiee_id.gt.0) THEN ! Ion induced electron emission
Nploc_init = partslist(p%iiee_id)%Nploc
CALL ion_induced(p, p%losthole, partslist(p%iiee_id), nblostparts)
Nploc_new = partslist(p%iiee_id)%Nploc
if (Nploc_new-Nploc_init .ge. 1) then
DO ii =Nploc_init+1,Nploc_new
Call p_calc_cellindex(partslist(p%iiee_id),ii)
call geom_weight(partslist(p%iiee_id)%pos(1,ii), partslist(p%iiee_id)%pos(3,ii), partslist(p%iiee_id)%geomweight(:,ii))
END DO
end if
END IF
DO i=nblostparts,1,-1
CALL delete_part(p,p%losthole(i))
END DO
END IF
!$OMP END MASTER
!$OMP BARRIER
END SUBROUTINE boundary_loss
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Computes the radial and axial cell index of the particle i
!> @param[in] p particles structure
!> @param[in] i index in p of the particle
!---------------------------------------------------------------------------
subroutine p_calc_cellindex(p,i)
use basic, only: rgrid,zgrid,invdz,invdr, nnr, nsubr,nsubz, nnz, invdthet, nr, nz
integer::i,j,k
type(particles)::p
k=0
p%cellindex(1,i)=0
do j=1,nsubr
IF (p%pos(1,i) .GE. rgrid(k) .AND. p%pos(1,i) .LT. rgrid(k+nnr(j))) THEN
p%cellindex(1,i)=floor((p%pos(1,i)-rgrid(k))*invdr(j))+k
exit
end if
k=k+nnr(j)
end do
if (p%pos(1,i) .GE. rgrid(nr)) p%cellindex(1,i)=nr-1
p%cellindex(2,i)=floor(p%pos(2,i)*invdthet)
k=0
p%cellindex(3,i)=0
do j=1,nsubz
IF (p%pos(3,i) .GE. zgrid(k) .AND. p%pos(3,i) .LT. zgrid(k+nnz(j))) THEN
p%cellindex(3,i)=floor((p%pos(3,i)-zgrid(k))*invdz(j))+k
exit
end if
k=k+nnz(j)
end do
if (p%pos(3,i) .GE. zgrid(nz)) p%cellindex(3,i)=nz-1
end subroutine p_calc_cellindex
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Computes the radial and axial cell index at position pos
!> @param[in] pos position
!> @param[in] index index of the cell
!---------------------------------------------------------------------------
subroutine calc_cellindex(index,pos)
use basic, only: rgrid,zgrid,invdz,invdr, nnr, nsubr,nsubz, nnz, invdthet, nr, nz
integer :: i,j,k
integer :: index(3)
real(kind=db) :: pos(3)
k=0
index=0
do j=1,nsubr
IF(pos(1) .GE. rgrid(k) .AND. pos(1) .LT. rgrid(k+nnr(j))) THEN
index(1)=floor((pos(1)-rgrid(k))*invdr(j))+k
exit
end if
k=k+nnr(j)
end do
if(pos(1) .GE. rgrid(nr)) index(1)=nr-1
index(2)=floor(pos(2)*invdthet)
k=0
do j=1,nsubz
IF(pos(3) .GE. zgrid(k) .AND. pos(3) .LT. zgrid(k+nnz(j))) THEN
index(3)=floor((pos(3)-zgrid(k))*invdz(j))+k
exit
end if
k=k+nnz(j)
end do
if(pos(3) .GE. zgrid(nz)) index(3)=nz-1
end subroutine calc_cellindex
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Computes the magnetic field amplitude at each particle position interpolated from the magnetic field at the closeset grid point
!> @param[in] p particles structure
!---------------------------------------------------------------------------
SUBROUTINE comp_mag_p(p)
USE basic, ONLY: zgrid, rgrid, BZ, BR, nz, nr, invdz
type(particles), INTENT(INOUT):: p
INTEGER :: i
Real(kind=db):: WZ,WR
INTEGER:: j1,j2,j3,j4
!$OMP DO SIMD
DO i=1,p%Nploc
WZ=(p%pos(3,i)-zgrid(p%cellindex(3,i)))/(zgrid(p%cellindex(3,i)+1)-zgrid(p%cellindex(3,i)));
WR=(p%pos(1,i)-rgrid(p%cellindex(1,i)))/(rgrid(p%cellindex(1,i)+1)-rgrid(p%cellindex(1,i)));
J1=(p%cellindex(3,i))*(nr+1) + p%cellindex(1,i)+1
J2=(p%cellindex(3,i))*(nr+1) + p%cellindex(1,i)+2
J3=(p%cellindex(3,i)+1)*(nr+1)+p%cellindex(1,i)+1
J4=(p%cellindex(3,i)+1)*(nr+1)+p%cellindex(1,i)+2
! Interpolation for magnetic field
p%B(2,i)=(1-WZ)*(1-WR)*Bz(J4) + WZ*(1-WR)*Bz(J3) + (1-WZ)*WR*Bz(J2) + WZ*WR*Bz(J1)
p%B(1,i)=(1-WZ)*(1-WR)*Br(J4) + WZ*(1-WR)*Br(J3) + (1-WZ)*WR*Br(J2) + WZ*WR*Br(J1)
END DO
!$OMP END DO SIMD NOWAIT
end subroutine comp_mag_p
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the classical simulations.
!> This routine systematically returns 1.0 to treat the system according to classical dynamic.
!
!> @param[out] gamma the lorentz factor \f$\gamma\f$
!> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity
!> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity
!> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity
!---------------------------------------------------------------------------
REAL(kind=db) FUNCTION gamma_classical(UZ, UR, UTHET)
!!#if __INTEL_COMPILER > 1700
!$OMP declare simd(gamma_classical)
!!#endif
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
gamma_classical=1.0
END FUNCTION gamma_classical
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the relativistic simulations.
!> This routine computes the Lorentz factor \f$\gamma=\sqrt{1+\mathbf{\gamma\beta}^2}\f$
!
!> @param[out] gamma the lorentz factor \f$\gamma\f$
!> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity
!> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity
!> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity
!---------------------------------------------------------------------------
REAL(kind=db) FUNCTION gamma_relativistic(UZ, UR, UTHET)
!!#if __INTEL_COMPILER > 1700
!$OMP declare simd(gamma_relativistic)
!!#endif
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
gamma_relativistic=sqrt(1+UZ**2+UR**2+UTHET**2)
END FUNCTION gamma_relativistic
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief General routine to compute the velocities at time t+1.
!> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint,
!> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle.
!
!> @param[in] p The particles structure being updated
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity(p)
USE basic, ONLY : nlclassical
type(particles), INTENT(INOUT):: p
! Store old Velocities
!$OMP MASTER
CALL swappointer2(p%Uold, p%U)
CALL swappointer(p%Gammaold, p%Gamma)
!$OMP END MASTER
!$OMP BARRIER
IF (nlclassical) THEN
CALL comp_velocity_fun(p, gamma_classical)
ELSE
CALL comp_velocity_fun(p, gamma_relativistic)
END IF
END SUBROUTINE comp_velocity
!---------------------------------------------------------------------------
!> @author
!> Patryk Kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Routine called by comp_velocity to compute the velocities at time t+1.
!> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint,
!> by using the routine computing gamma as an input. This avoid the nlclassical flag check on each particle.
!
!> @param[in] gamma the function used to compute the value of the lorentz factor \f$\gamma\f$
!> @param[in] p The particles structure being updated
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity_fun(p, gammafun)
USE basic, ONLY : vnorm, bnorm, rnorm, tnorm, dt
procedure(gamma)::gammafun
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: B, cosa, sina, Udrift
INTEGER:: i, iend, nbunch, j
nbunch=16
IF (p%Nploc .NE. 0) THEN
!$OMP DO
DO i=1,p%Nploc,nbunch
iend = min(i + nbunch - 1, p%Nploc)
DO j=i,iend
B=sqrt(p%B(1,j)**2+p%B(2,j)**2)
cosa=p%B(2,j)/B
sina=p%B(1,j)/B
p%U(2,j)=-p%E(1,j)/p%B(2,j)
Udrift=p%collfreq(j)/(-p%qmRatio*B*bnorm)*p%U(2,j)
p%larmorrad(j)=sqrt(p%U(2,j)**2+Udrift**2)*vnorm/(-p%qmRatio*B*bnorm*rnorm)
p%Upar(j)=p%Upar(j)+p%qmRatio*dt*bnorm*tnorm*(p%E(3,j)*cosa+p%E(1,j)*sina)
p%U(1,j)=cosa*Udrift+sina*p%Upar(j)
p%U(3,j)=-sina*Udrift+cosa*p%Upar(j)
! Computation of the Lorentz factor
p%Gamma(j)=gammafun(p%U(3,i), p%U(1,i), p%U(2,i))
END DO
END DO
!$OMP END DO NOWAIT
END IF
p%collected=.false.
END SUBROUTINE comp_velocity_fun
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief General routine to compute the velocities at time t+1.
!> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint,
!> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle.
!
!> @param[in] p The particles structure being updated
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity_ions(p)
USE basic, ONLY : vnorm, bnorm, rnorm, tnorm, dt
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: tau
REAL(kind=db):: ZBR(16), ZBZ(16), ZPR(16), ZPZ(16), ZPTHET(16), SQR(16), ZBZ2(16), ZBR2(16)
INTEGER:: i, iend, nbunch, j,k
nbunch=16
! Store old Velocities
!$OMP MASTER
CALL swappointer2(p%Uold, p%U)
CALL swappointer(p%Gammaold, p%Gamma)
!$OMP END MASTER
!$OMP BARRIER
! Normalized time increment
tau=p%qmRatio*bnorm*0.5*dt*tnorm
IF(p%Nploc .NE. 0) THEN
!$OMP DO
DO i=1,p%Nploc,nbunch
iend = min(i + nbunch - 1, p%Nploc)
DO j=i,iend
k=j-i+1
! First half of electric pulse
p%U(3,j)=p%Uold(3,j)+p%E(3,j)*tau
p%U(1,j)=p%Uold(1,j)+p%E(1,j)*tau
p%Gamma(j)=1.0
! Rotation along magnetic field
ZBZ(k)=tau*p%B(2,j)/p%Gamma(j)
ZBR(k)=tau*p%B(1,j)/p%Gamma(j)
ZPZ(k)=p%U(3,j)-ZBR(k)*p%Uold(2,j) ! u'_{z}
ZPR(k)=p%U(1,j)+ZBZ(k)*p%Uold(2,j) ! u'_{r}
ZPTHET(k)=p%Uold(2,j)+(ZBR(k)*p%U(3,j)-ZBZ(k)*p%U(1,j)) ! u'_{theta}
SQR(k)=1+ZBZ(k)*ZBZ(k)+ZBR(k)*ZBR(k)
ZBZ2(k)=2*ZBZ(k)/SQR(k)
ZBR2(k)=2*ZBR(k)/SQR(k)
p%U(3,j)=p%U(3,j)-ZBR2(k)*ZPTHET(k) ! u+_{z}
p%U(1,j)=p%U(1,j)+ZBZ2(k)*ZPTHET(k) ! u+_{r}
p%U(2,j)=p%Uold(2,j)+(ZBR2(k)*ZPZ(k)-ZBZ2(k)*ZPR(k)) ! u+_{theta}
! Second half of acceleration
p%U(3,j)=p%U(3,j)+p%E(3,j)*tau
p%U(1,j)=p%U(1,j)+p%E(1,j)*tau
! Final computation of the Lorentz factor
p%Gamma(j)=1.0
END DO
END DO
!$OMP END DO NOWAIT
END IF
p%collected=.false.
END SUBROUTINE comp_velocity_ions
!---------------------------------------------------------------------------
!> @author
!> Patryk Kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
!DESCRIPTION:
!> @brief Computes the particles position at time t+1
!> This routine computes the particles position at time t+1 according to the Bunemann algorithm.
!
!> @param[in] p The particles structure being updated
! ---------------------------------------------------------------------------
SUBROUTINE push(p)
Use basic, ONLY: dt
type(particles), INTENT(INOUT):: p
REAL(kind=db):: xp, yp, cosa, sina, alpha, vr, vthet
INTEGER :: i
IF(p%Nploc .NE. 0) THEN
!$OMP DO SIMD
DO i=1,p%Nploc
! Local Cartesian coordinates
xp=p%pos(1,i)+dt*p%U(1,i)/p%Gamma(i)
IF(xp.le.0)THEN
xp=abs(xp)
p%U(1,i)=-p%U(1,i)
END IF
yp=dt*p%U(2,i)/p%Gamma(i)
! Conversion to cylindrical coordiantes
p%pos(3,i)=p%pos(3,i)+dt*p%U(3,i)/p%Gamma(i)
p%pos(1,i)=sqrt(xp**2+yp**2)
! Computation of the rotation angle
IF(p%pos(1,i) .EQ. 0) THEN
cosa=1
sina=0
alpha=0
ELSE
cosa=xp/p%pos(1,i)
sina=yp/p%pos(1,i)
alpha=asin(sina)
END IF
! New azimuthal position
p%pos(2,i)=p%pos(2,i)+alpha
IF(p%pos(2,i) .lt. 0.0) p%pos(2,i)=mod(p%pos(2,i),2*pi)+2*pi
IF(p%pos(2,i) .ge. 2*pi) p%pos(2,i)=mod(p%pos(2,i),2*pi)
! Velocity in rotated reference frame
vr=cosa*p%U(1,i)+sina*p%U(2,i)
vthet=-sina*p%U(1,i)+cosa*p%U(2,i)
p%U(1,i)=vr
p%U(2,i)=vthet
END DO
!$OMP END DO SIMD NOWAIT
END IF
!$OMP SINGLE
p%collected=.false.
!$OMP END SINGLE NOWAIT
END SUBROUTINE push
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Computes several diagnostic quantities
!> This routine computes the total kinetic and electric potential energy.
!> It keeps track of the reference energy and the number of particle per mpi node.
!
!---------------------------------------------------------------------------
SUBROUTINE partdiagnostics
!
! Compute energies
!
USE constants, ONLY: vlight
USE basic, ONLY: phinorm, cstep, nlclassical, ierr, nbspecies
INTEGER:: i,j
! Reset the quantities
!$OMP SINGLE
ekin=0
epot=0
etot=0
!$OMP END SINGLE NOWAIT
! Computation of the kinetic and potential energy as well as fluid velocities and density
!!$OMP PARALLEL DO REDUCTION(+:epot, ekin) DEFAULT(SHARED), PRIVATE(i,j)
Do j=1,nbspecies
if(.not. partslist(j)%is_field) CYCLE
!$OMP DO reduction(+:epot,ekin)
DO i=1,partslist(j)%Nploc
! Potential energy
epot=epot+(partslist(j)%pot(i)+partslist(j)%potxt(i))*partslist(j)%q*partslist(j)%weight
! Kinetic energy
IF(.not. nlclassical) THEN
ekin=ekin+(0.5*(partslist(j)%Gammaold(i)+partslist(j)%Gamma(i))-1)*partslist(j)%m*partslist(j)%weight
ELSE
ekin=ekin+0.5*( partslist(j)%U(1,i)*partslist(j)%Uold(1,i) &
& + partslist(j)%U(3,i)*partslist(j)%Uold(3,i) &
& + partslist(j)%U(2,i)*partslist(j)%Uold(2,i) )*partslist(j)%m*partslist(j)%weight
END IF
END DO
!$OMP END DO NOWAIT
END DO
!$OMP BARRIER
!$OMP MASTER
!!$OMP END PARALLEL DO
epot=epot*phinorm*0.5
ekin=ekin*vlight**2
! Shift to Etot at cstep=1 (not valable yet at cstep=0!)
IF(cstep.EQ. 1) THEN
! Compute the local total energy
loc_etot0 = epot+ekin
etot0=0
END IF
! Compute the total energy
etot=epot+ekin
Energies=(/ekin,epot,etot,loc_etot0/)
! The computed energy is sent to the root process
IF(mpisize .gt.1) THEN
IF(mpirank .eq.0 ) THEN
CALL MPI_REDUCE(MPI_IN_PLACE, Energies, 4, db_type, db_sum_op, &
& 0, MPI_COMM_WORLD, ierr)
etot0=etot0+Energies(4)
ekin=Energies(1)
epot=Energies(2)
etot=Energies(3)
ELSE
CALL MPI_REDUCE(Energies, Energies, 4, db_type, db_sum_op, &
& 0, MPI_COMM_WORLD, ierr)
END IF
ELSE
etot0=etot0+loc_etot0
END IF
loc_etot0=0
! Send the local number of particles on each node to the root process
IF(mpisize .gt. 1) THEN
Nplocs_all(mpirank)=partslist(1)%Nploc
IF(mpirank .eq.0 ) THEN
CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
!CALL MPI_REDUCE(MPI_IN_PLACE,partslist(1)%nudcol,3,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr)
partslist(1)%Nptot=sum(Nplocs_all)
!partslist(1)%nudcol=partslist(1)%nudcol/partslist(1)%Nptot
ELSE
CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
!CALL MPI_REDUCE(partslist(1)%nudcol,partslist(1)%nudcol,3,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr)
END IF
ELSE
partslist(1)%Nptot=partslist(1)%Nploc
END IF
!$OMP END MASTER
!$OMP BARRIER
end subroutine partdiagnostics
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Collect the particles positions and velocities on the root process.
!> If the collection has already been performed at this time step, the routine does nothing.
!
!---------------------------------------------------------------------------
SUBROUTINE collectparts(p,p_towrite)
USE basic, ONLY: mpirank, mpisize, ierr
type(particles), INTENT(INOUT):: p
type(particles), INTENT(INOUT), OPTIONAL:: p_towrite
INTEGER, DIMENSION(:), ALLOCATABLE :: displs, Nploc
INTEGER:: i
INTEGER:: particles_type(0:mpisize-1) !< Stores the MPI data type used for particles gathering on node 0 and broadcast from node 0
INTEGER :: part_requests(0:mpisize-1)
INTEGER:: stats(MPI_STATUS_SIZE,mpisize)
INTEGER:: root_send_particles_type
INTEGER:: root_send_part_request
part_requests=MPI_REQUEST_NULL
root_send_part_request=MPI_REQUEST_NULL
particles_type=MPI_DATATYPE_NULL
root_send_particles_type=MPI_DATATYPE_NULL
IF(p%collected) RETURN ! exit subroutine if particles have already been collected during this time step
ALLOCATE(Nploc(0:mpisize-1))
ALLOCATE(displs(0:mpisize-1))
displs=0
Nploc(mpirank)=p%Nploc
CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, Nploc, 1, MPI_INTEGER,&
& MPI_COMM_WORLD, ierr)
p%Nptot=sum(Nploc)
IF(p%Nptot .eq. 0 ) THEN
p%partindex=-1
p%collected=.true.
RETURN
END IF
Do i=1,mpisize-1
displs(i)=displs(i-1)+Nploc(i-1)
END DO
IF (present(p_towrite))then
IF(mpirank.eq.0 .and. p%Nptot .gt. size(p_towrite%pos,2)) THEN
CALL change_parts_allocation(p_towrite,max(p%Nptot-size(p_towrite%pos,2),floor(0.5*size(p_towrite%pos,2))))
END IF
else
IF(mpirank.eq.0 .and. p%Nptot .gt. size(p%pos,2)) THEN
CALL change_parts_allocation(p,max(p%Nptot-size(P%pos,2),floor(0.5*size(P%pos,2))))
END IF
end if
IF(mpirank .ne. 0) THEN
if(Nploc(mpirank) .gt. 0) THEN
Call init_particles_gather_mpi(p,1,Nploc(mpirank),particles_type(mpirank))
! Send Particles informations to root process
CALL MPI_SEND(p, 1, particles_type(mpirank), 0, partsgather_tag, MPI_COMM_WORLD, ierr)
CALL MPI_TYPE_FREE(particles_type(mpirank),ierr)
END IF
ELSE
! Receive particle information from all processes
if(present(p_towrite))then
if(Nploc(mpirank) .gt. 0) THEN
Call init_particles_gather_mpi(p,1,Nploc(mpirank),root_send_particles_type)
! Send Particles informations to root process
CALL MPI_ISEND(p, 1, root_send_particles_type, 0, partsgather_tag, MPI_COMM_WORLD, root_send_part_request, ierr)
end if
DO i=0,mpisize-1
if(Nploc(i) .lt. 1) cycle
Call init_particles_gather_mpi(p_towrite,displs(i)+1,Nploc(i),particles_type(i))
CALL MPI_IRECV(p_towrite,1,particles_type(i),i,partsgather_tag,MPI_COMM_WORLD, part_requests(i), ierr)
END DO
CALL MPI_WAITALL(mpisize,part_requests, stats, ierr)
if(Nploc(mpirank).gt. 0) CALL MPI_TYPE_FREE(root_send_particles_type,ierr)
p_towrite%partindex(sum(Nploc)+1:)=-1
else
DO i=1,mpisize-1
if(Nploc(i) .lt. 1) cycle
Call init_particles_gather_mpi(p,displs(i)+1,Nploc(i),particles_type(i))
CALL MPI_IRECV(p,1,particles_type(i),i,partsgather_tag,MPI_COMM_WORLD, part_requests(i), ierr)
END DO
CALL MPI_WAITALL(mpisize,part_requests, stats, ierr)
p%partindex(sum(Nploc)+1:)=-1
end if
Do i=1,mpisize-1
if(Nploc(i) .lt. 1) cycle
CALL MPI_TYPE_FREE(particles_type(i),ierr)
END DO
END IF
p%collected=.TRUE.
END SUBROUTINE collectparts
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Computes the velocities at time t-1/2 delta t to keep the second order precision in time on the velocity.
!> This should only be used at particle initialisation time, ot in the case of a restart.
!
!---------------------------------------------------------------------------
SUBROUTINE adapt_vinit(p)
!! Computes the velocity at time -dt/2 from velocities computed at time 0
!
USE basic, ONLY : bnorm, dt, tnorm, nlclassical, phinorm, distribtype, vnorm
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, &
& SQR, Vperp, v2
INTEGER :: i
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
! In case Davidson distribution is used the longitudinal and radial velocities are adapted to take into account the
! electric potential.
IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4 .or. p%Davidson) THEN
ALLOCATE(VR(p%Nploc),VZ(p%Nploc),VTHET(p%Nploc))
CALL loduni(7,VZ)
VZ=VZ*2*pi
VTHET=p%U(2,:)/p%Gamma*vnorm
DO i=1,p%Nploc
Vperp=sqrt(MAX(2*p%H0/p%m-2*p%qmRatio*p%pot(i)*phinorm-VTHET(i)**2,0.0_db))
VR(i)=Vperp*sin(VZ(i))
VZ(i)=Vperp*cos(VZ(i))
IF(nlclassical) THEN
p%Gamma(i)=1
ELSE
v2=VR(i)**2+VZ(i)**2+VTHET(i)**2
p%Gamma(i)=sqrt(1/(1-v2/vnorm**2))
END IF
p%U(1,i)=p%Gamma(i)*VR(i)/vnorm
p%U(3,i)=p%Gamma(i)*VZ(i)/vnorm
p%U(2,i)=p%Gamma(i)*VTHET(i)/vnorm
END DO
DEALLOCATE(VR,VZ,VTHET)
END IF
! Normalised time increment
!tau=-omegac/2/omegap*dt/tnorm
tau=-p%qmRatio*bnorm*0.5*dt*tnorm
! Store old Velocities
CALL swappointer2(p%Uold, p%U)
!CALL swappointer(p%URold, p%UR)
!CALL swappointer(p%UTHETold, p%UTHET)
CALL swappointer(p%Gammaold, p%Gamma)
IF (p%Nploc .NE. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR)
DO i=1,p%Nploc
! Half inverse Rotation along magnetic field
ZBZ=tau*p%B(2,i)/p%Gammaold(i)
ZBR=tau*p%B(1,i)/p%Gammaold(i)
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZPZ=(p%Uold(3,i)-ZBR*p%Uold(2,i))/SQR !u-_{z}
ZPR=(p%Uold(1,i)+ZBZ*p%Uold(2,i))/SQR !u-_{r}
ZPTHET=p%Uold(2,i)+(ZBR*p%Uold(3,i)-ZBZ*p%Uold(1,i))/SQR !u-_{theta}
p%U(3,i)=ZPZ
p%U(1,i)=ZPR
p%U(2,i)=ZPTHET
! half of decceleration
p%U(3,i)=p%U(3,i)+p%E(3,i)*tau
p%U(1,i)=p%U(1,i)+p%E(1,i)*tau
IF(.not. nlclassical) THEN
p%Gamma(i)=sqrt(1+p%U(3,i)**2+p%U(1,i)**2+p%U(2,i)**2)
END IF
END DO
!$OMP END PARALLEL DO SIMD
END IF
END SUBROUTINE adapt_vinit
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Calculates the number of particles per column of the spatial grid ( at fixed axial cell position)
!> This facilitate the computation of the axial grid limits for each MPI worker
!
!---------------------------------------------------------------------------
SUBROUTINE calcnbperz(p,nbperz)
USE basic, only: nz
IMPLICIT NONE
type(particles):: p
INTEGER, INTENT(INOUT):: nbperz(0:)
Integer::i, zindex
!$OMP PARALLEL DO DEFAULT(SHARED) reduction(+:nbperz), private(zindex,i)
Do i=1,p%Nploc
! we make sure zindex is in [0, nz-1] to avoid segmentation faults
zindex=min(nz-1,max(p%cellindex(3,i),0))
nbperz(zindex)=nbperz(zindex)+1
END DO
!$OMP END PARALLEL DO
END SUBROUTINE calcnbperz
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief In the case of MPI parallelism, computes the axial limits of the local domain.
!---------------------------------------------------------------------------
SUBROUTINE calc_Zbounds(p, Zbounds, norder)
! Computes the start and end indices for the Z boundaries on local processus
! Computes the particle indices from initial particle loading vector, that stay in current process
USE basic, ONLY: nz, cstep, mpirank, mpisize, step, nbspecies
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER:: Zbounds(0:)
INTEGER:: norder(2)
INTEGER:: old_Zbounds(0:size(Zbounds,1)-1)
INTEGER:: k, i, nbparts
REAL(kind=db):: idealnbpartsperproc
INTEGER, DIMENSION(0:nz-1):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1)
INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction
INTEGER:: Zperproc, ierr, remparts
CHARACTER(12)::fmt
partspercol=0
! calculate the axial distribution integrated along the radial direction
do i=1,nbspecies
call calcnbperz(partslist(i),partspercol)
end do
! gather this data on all nodes
if(step .gt. 0 .and. mpisize .gt. 1) THEN
old_Zbounds=Zbounds
CALL MPI_ALLREDUCE(MPI_IN_PLACE, partspercol, nz, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
END IF
! estimate the ideal number of particle per MPI worker
idealnbpartsperproc = 0
DO i=1,nbspecies
idealnbpartsperproc = idealnbpartsperproc+partslist(i)%Nptot/mpisize
END DO
! find the start and end indices where particles are present
Zmin=0
Zmax=nz-1
Do k=0,nz-1
if(partspercol(k) .gt.0) then
Zmin=k
exit
end if
end do
Do k=nz-1,0,-1
if(partspercol(k) .gt.0) then
Zmax=k
exit
end if
end do
! Find naive axial limits assuming uniform axial distribution
IF(Zmax .le. 0) Zmax=nz-1
IF(Zmin .gt. nz) Zmin=0
Zperproc=(Zmax-Zmin)/mpisize
IF (Zperproc .lt. 1) THEN
!! No particles are present initially
Zperproc=nz/mpisize
Zmin=0
! Define boundaries using naive guess on start or restart (allow to start with 0 parts)
DO k=1,mpisize-1
IF(k .lt. mpisize-1-MODULO(Zmax-Zmin,mpisize)) THEN
Zbounds(k)=Zmin+k*Zperproc-1
ELSE
Zbounds(k)=Zmin+k*Zperproc-1+k-mpisize+2+MODULO(Zmax-Zmin,mpisize)
END IF
END DO
ELSE
i=0
! Define axial boundaries using the axial distribution information.
! the subdomains are not equal
remparts=p%Nptot
DO k=1,mpisize-1
nbparts=0
DO WHILE(nbparts<0.98*idealnbpartsperproc .and. i .lt. Zmax .and. (nbparts+partspercol(i)).lt.1.05*idealnbpartsperproc)
nbparts=nbparts+partspercol(i)
i=i+1
END DO
remparts=remparts-nbparts
Zbounds(k)=i
END DO
END IF
IF(step .gt. 0 .and. mpirank .eq. 0) THEN
Do i=1,mpisize-1
!We check that the new limits will not exceed the old limits of the left and right process
! This avoids particle communications with process >mpirank+2 and < mpirank-2
! However this should converge over time
IF(Zbounds(i) .lt. old_Zbounds(i-1)) Zbounds(i)=old_Zbounds(i-1)
if(Zbounds(i) .gt. old_Zbounds(i+1))Zbounds(i)=old_Zbounds(i+1)
! If a process would have an axial domain shoter than axial norder, we revert to the old boundaries.
IF((Zbounds(i)-Zbounds(i-1)).lt. norder(1) .or. (Zbounds(i+1)-Zbounds(i)).lt. norder(1)) THEN
Zbounds=old_Zbounds
EXIT
END IF
END DO
END IF
! send the new boundaries to all the workers
CALL MPI_Bcast(Zbounds,mpisize+1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr)
DO k=0,mpisize-1
Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1))
END DO
if(mpirank .eq. 0) THEN
WRITE(fmt,'(a,i3,a)')"(a,",mpisize+1, "i5)"
WRITE(*,fmt) "Zbounds: ", Zbounds
WRITE(fmt,'(a,i3,a)')"(a,",mpisize, "i8)"
WRITE(*,fmt) "Nplocs: ", Nplocs_all
END IF
END SUBROUTINE calc_Zbounds
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief After a restart keep only the particles in the local domain of the current MPI worker
!---------------------------------------------------------------------------
SUBROUTINE keep_mpi_self_parts(p,Zbounds)
TYPE(particles),INTENT(INOUT):: p
INTEGER,INTENT(in)::Zbounds(0:)
INTEGER :: i, partstart, old_sum,ierr
partstart=1
p%Nploc=0
Do i=1,p%Nptot
IF(p%cellindex(3,i) .ge.Zbounds(mpirank).and.p%cellindex(3,i).lt.Zbounds(mpirank+1)) THEN
p%Nploc=p%Nploc+1
CALL move_part(p,i,p%Nploc)
END IF
END DO
old_sum=p%Nptot
CALL MPI_REDUCE(p%Nploc, p%Nptot,1,MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
IF(p%Nptot .ne. old_sum) THEN
WRITE(*,*) "Error in particle distribution kept: ", p%Nptot, "/",old_sum
!call MPI_Abort(MPI_COMM_WORLD, -1, ierr)
!stop
END IF
END SUBROUTINE keep_mpi_self_parts
!_______________________________________________________________________________
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Manage the particle communication between neighbours.
!> This routine is responsible to receive the incoming particles from the MPI neighbours and to send its outgoing
!> particles to these neighbours
!
!> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1)
!> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1)
!---------------------------------------------------------------------------
SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, receivednbparts, procs)
USE mpihelper, ONLY: particle_type
#ifdef _DEBUG
USE basic, ONLY: step
#endif
type(particles), INTENT(INOUT):: p
INTEGER, INTENT(inout) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(out) :: receivednbparts
INTEGER, INTENT(in) :: procs(2)
INTEGER :: rrecvnbparts=0, lrecvnbparts=0
INTEGER :: sendrequest(2), recvrequest(2)
INTEGER :: sendstatus(MPI_STATUS_SIZE,2), recvstatus(MPI_STATUS_SIZE,2)
TYPE(particle), ALLOCATABLE :: rrecvpartbuff(:), lrecvpartbuff(:), rsendpartbuff(:), lsendpartbuff(:) ! buffers to send and receive particle from left and right processes
INTEGER :: lsentnbparts, rsentnbparts
INTEGER :: lreceivednbparts, rreceivednbparts, ierr
lsentnbparts=lsendnbparts
rsentnbparts=rsendnbparts
sendrequest=MPI_REQUEST_NULL
recvrequest=MPI_REQUEST_NULL
lrecvnbparts=0
rrecvnbparts=0
! Send and receive the number of particles to exchange
if (procs(1).ge. 0) then
CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, procs(1), nbpartsexchange_tag, MPI_COMM_WORLD, recvrequest(1), ierr)
CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, procs(1), nbpartsexchange_tag, MPI_COMM_WORLD, sendrequest(1), ierr)
end if
if (procs(2).ge. 0) then
CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, procs(2), nbpartsexchange_tag, MPI_COMM_WORLD, recvrequest(2), ierr)
CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, procs(2), nbpartsexchange_tag, MPI_COMM_WORLD, sendrequest(2), ierr)
end if
CALL MPI_Waitall(2,recvrequest(1:2), recvstatus(:,1:2), ierr)
recvrequest=MPI_REQUEST_NULL
lreceivednbparts=lrecvnbparts
rreceivednbparts=rrecvnbparts
! Re/allocate enough memory to store the incoming particles
ALLOCATE(rrecvpartbuff(rreceivednbparts))
ALLOCATE(lrecvpartbuff(lreceivednbparts))
! Receive particles from left and right processes to the corresponding buffers
IF ( lrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, procs(1), partsexchange_tag, MPI_COMM_WORLD, recvrequest(1), ierr)
END IF
IF( rrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, procs(2), partsexchange_tag, MPI_COMM_WORLD, recvrequest(2), ierr)
END IF
ALLOCATE(rsendpartbuff(rsendnbparts))
ALLOCATE(lsendpartbuff(lsendnbparts))
! Copy the leaving particles to the corresponding send buffers
IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN
CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, lsendpartbuff, rsendpartbuff)
END IF
CALL MPI_Waitall(2,sendrequest(1:2), sendstatus(:,1:2), ierr)
! Send the particles to the left and right neighbours
IF( lsendnbparts .gt. 0) THEN
CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, procs(1), partsexchange_tag, MPI_COMM_WORLD, sendrequest(1), ierr)
#ifdef _DEBUG
!WRITE(*,*)"snding ", lsendnbparts , " to left at step: ",step
#endif
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, procs(2), partsexchange_tag, MPI_COMM_WORLD, sendrequest(2), ierr)
#ifdef _DEBUG
!WRITE(*,*)"snding ", rsendnbparts , " to right at step: ",step
#endif
END IF
receivednbparts=rreceivednbparts+lreceivednbparts
IF(p%Nploc+receivednbparts-lsendnbparts-rsendnbparts .gt. size(p%pos,2)) THEN
CALL change_parts_allocation(p,receivednbparts)
END IF
! Receive the incoming parts in the receive buffers
IF ( lreceivednbparts .gt. 0) THEN
CALL MPI_Wait(recvrequest(1), recvstatus(:,1), ierr)
IF(ierr .ne. MPI_SUCCESS) THEN
WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,1)
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
#ifdef _DEBUG
!WRITE(*,*)"rcvd ", lreceivednbparts , " from left at step: ",step
#endif
END IF
IF ( rreceivednbparts .gt. 0) THEN
CALL MPI_Wait(recvrequest(2), recvstatus(:,2), ierr)
IF(ierr .ne. MPI_SUCCESS) THEN
WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,2)
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
#ifdef _DEBUG
!WRITE(*,*)"rcvd ", rreceivednbparts , " from right at step: ",step
#endif
END IF
! Copy the incoming particles from the receive buffers to the simulation parts variable
CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, &
& lrecvpartbuff, rrecvpartbuff)
! Wait for the outgoing particles to be fully received by the neighbours
IF( lsendnbparts .gt. 0) THEN
CALL MPI_Wait(sendrequest(1), sendstatus(:,1), ierr)
#ifdef _DEBUG
!WRITE(*,*)"sent ", lsentnbparts , " to left at step: ",step
#endif
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_Wait(sendrequest(2), sendstatus(:,2), ierr)
#ifdef _DEBUG
!WRITE(*,*)"sent ", rsentnbparts , " to right at step: ",step
#endif
END IF
!
!
END SUBROUTINE particlescommunication
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Copy the particles from the receive buffers to the local simulation variable parts.
!> The incoming particles will first be stored in the holes left by the outgoing particles, then they
!> will be added at the end of the parts variable
!
!> @param [in] rrecvnbparts number of particles received from the right neighbour (mpirank+1)
!> @param [in] lrecvnbparts number of particles received from the left neighbour (mpirank-1)
!> @param [in] sendnbparts total number of particles having left the local domain
!---------------------------------------------------------------------------
SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts,lrecvpartbuff, rrecvpartbuff)
!
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts
TYPE(particle), INTENT(IN) :: rrecvpartbuff(:), lrecvpartbuff(:)
INTEGER k,partpos
! First import the particles coming from the right
IF(rrecvnbparts .gt. 0) THEN
Do k=1,rrecvnbparts
IF(k .le. sendnbparts) THEN
! Fill the holes left by sent parts
partpos=abs(p%sendhole(k))
ELSE
! Add at the end of parts and keep track of number of parts
p%Nploc=p%Nploc+1
partpos=p%Nploc
END IF
CALL Insertincomingpart(p, rrecvpartbuff(k), partpos)
END DO
END IF
! Then import the particles coming from the left
IF(lrecvnbparts .gt. 0) THEN
Do k=1,lrecvnbparts
IF(k+rrecvnbparts .le. sendnbparts) THEN
! Fill the holes left by sent parts
partpos=abs(p%sendhole(k+rrecvnbparts))
ELSE
! Add at the end of parts and keep track of number of parts
p%Nploc=p%Nploc+1
partpos=p%Nploc
END IF
CALL Insertincomingpart(p, lrecvpartbuff(k), partpos)
END DO
END IF
!
END SUBROUTINE Addincomingparts
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Copy the particles from the local parts variable to the left and right send buffers.
!
!> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1)
!> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1)
!---------------------------------------------------------------------------
SUBROUTINE AddPartSendBuffers(p, lsendnbparts, rsendnbparts, lsendpartbuff, rsendpartbuff)
!
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
TYPE(particle), INTENT(OUT) :: rsendpartbuff(:), lsendpartbuff(:)
INTEGER:: partpos, k
INTEGER:: lsendpos, rsendpos
lsendpos=0
rsendpos=0
! Loop over the outgoing particles and fill the correct send buffer
Do k=lsendnbparts+rsendnbparts,1,-1
partpos=abs(p%sendhole(k))
IF(p%sendhole(k) .GT. 0) THEN
rsendpos=rsendpos+1
CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos)
ELSE IF(p%sendhole(k) .LT. 0) THEN
lsendpos=lsendpos+1
CALL Insertsentpart(p, lsendpartbuff, lsendpos, partpos)
END IF
END DO
!
!
END SUBROUTINE AddPartSendBuffers
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Add the particles stored in the buffer to the main particle storage p in particles form
!> @param[in] p particles structure to add particles to
!> @param[in] buffer memory containing the particles to be added
!> @param[in] nb_ins number of particles stored in buffer
!---------------------------------------------------------------------------
SUBROUTINE add_list_created_part(p, buffer,nb_ins)
IMPLICIT NONE
TYPE(particles), INTENT(INOUT):: p
TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:)
INTEGER, OPTIONAL:: nb_ins
INTEGER:: i, nptotinit, parts_size_increase, nb_added
nptotinit=p%Nploc+1
if(present(nb_ins)) THEN
nb_added=nb_ins
ELSE
nb_added=size(buffer,1)
end if
IF(nb_added .le. 0) RETURN ! No particles to add
! if there is not enough space in the parts simulation buffer, increase the parst size
IF(p%Nploc + nb_added .gt. size(p%pos,2)) THEN
parts_size_increase=Max(floor(0.1*size(p%pos,2)),nb_added)
CALL change_parts_allocation(p, parts_size_increase)
END IF
DO i=1,nb_added
CALL add_created_particle(p,buffer(i))
END DO
nb_added=p%Nploc-nptotinit+1
if(p%is_field) then
IF(allocated(p%addedlist)) then
call change_array_size_int(p%addedlist,2)
else
allocate(p%addedlist(2))
end if
p%addedlist(size(p%addedlist)-1)=nptotinit
p%addedlist(size(p%addedlist))=nb_added
end if
END SUBROUTINE add_list_created_part
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Add the particles stored in the linked buffer to the main particle storage p in particles form
!> @param[in] p particles structure to add particles to
!> @param[in] linked_buffer memory containing the particles to be added in linked list format
!> @param[in] destroy Indicates if the memory of the linked buffer must be freed after copy to p
!> @param[in] zerovelocity Define if the velocity of the particles in p is set to 0 or copied from the buffer
!---------------------------------------------------------------------------
SUBROUTINE add_linked_created_part(p, linked_buffer, destroy, zerovelocity)
IMPLICIT NONE
TYPE(particles), INTENT(INOUT):: p
TYPE(linked_part_row), INTENT(in) :: linked_buffer
LOGICAL:: destroy, zerovelocity
TYPE(linked_part), POINTER:: part
INTEGER:: i, nptotinit, parts_size_increase, nb_added
nptotinit=p%Nploc+1
nb_added=linked_buffer%n
IF(nb_added .le. 0) RETURN ! No particles to add
! if there is not enough space in the parts simulation buffer, increase the parst size
IF(p%Nploc + nb_added .gt. size(p%pos,2)) THEN
parts_size_increase=Max(floor(0.2*size(p%pos,2)),nb_added)
CALL change_parts_allocation(p, parts_size_increase)
END IF
part=>linked_buffer%start
DO i=1,nb_added
CALL add_created_particle(p,part%p)
part=>part%next
END DO
nb_added=p%Nploc-nptotinit+1
if(p%is_field) then
IF(allocated(p%addedlist)) then
call change_array_size_int(p%addedlist,2)
else
allocate(p%addedlist(2))
end if
p%addedlist(size(p%addedlist)-1)=nptotinit
p%addedlist(size(p%addedlist))=nb_added
end if
if(zerovelocity)then
p%U(:,nptotinit:p%Nploc)=0
!p%UTHET(nptotinit:p%Nploc)=0
!p%UZ(nptotinit:p%Nploc)=0
end if
if (destroy) call destroy_linked_parts(linked_buffer%start)
if (p%is_field) then
! we keep track of energy by removing the ionisation energy
! with conversion from electronvolt to joules
loc_etot0=loc_etot0-sum(p%pot(nptotinit:p%Nploc)*elchar)
end if
END SUBROUTINE add_linked_created_part
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Add created particles from a buffer of type particle to the main species storages.
!
!> @param [in] p specie memory where we want to add particles
!> @param [in] part particle buffer storing the data we want to add to p
!---------------------------------------------------------------------------
SUBROUTINE add_created_particle(p,part)
USE geometry
TYPE(particles):: p
TYPE(particle):: part
p%Nploc=p%Nploc+1
p%newindex=p%newindex+1
! add the data to the p structure
CALL Insertincomingpart(p, part, p%Nploc)
p%partindex(p%Nploc)=p%newindex
! calculate the new domain weight
CALL dom_weight(p%pos(1,p%Nploc),p%pos(3,p%Nploc),p%geomweight(0,p%Nploc))
! delete the particle if it is outside of the computational domain
if( .not. is_inside(p,p%Nploc) ) then
p%Nploc=p%Nploc-1
p%newindex=p%newindex-1
RETURN
end if
! Calculate the geometric weight for the Poisson solver and the grid indices
CALL geom_weight(p%pos(1,p%Nploc),p%pos(3,p%Nploc),p%geomweight(:,p%Nploc))
call p_calc_cellindex(p,p%Nploc)
END SUBROUTINE add_created_particle
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Checks if the particle id in p is inside of the simulation domain
!
!> @param [in] p specie memory
!> @param [in] id index of the particle we want to test
!---------------------------------------------------------------------------
function is_inside(p,id)
Use basic, ONLY: rgrid,zgrid, nr, nz
IMPLICIT NONE
logical :: is_inside
type(particles) :: p
integer :: id
is_inside=.true.
! Check if the particle is in the simulation domain
if(p%geomweight(0,id).le.0)then
is_inside=.false.
return
end if
! check if the particle is in the simulation grid
if(p%pos(1,id).ge.rgrid(nr) .or. p%pos(1,id) .le. rgrid(0))then
is_inside=.false.
return
end if
if(p%pos(3,id).ge.zgrid(nz) .or. p%pos(3,id) .le. zgrid(0))then
is_inside=.false.
return
end if
end function is_inside
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Calculate the energy added by new particles to the system for diagnostic purposes
!
!> @param [in] p specie memory
!---------------------------------------------------------------------------
SUBROUTINE calc_newparts_energy(p)
USE basic, ONLY: phinorm, nlclassical
USE omp_lib
type(particles)::p
integer::i,n,nptotinit,nbadded, nptotend, j
Real(kind=db):: addedE
! exit if these particles dont participate in the Poisson equation
if( allocated(p%addedlist)) then
addedE=0.0_db
n=size(p%addedlist)
!$OMP SINGLE
nbadded=0
Do i=1,n,2
nbadded=nbadded+p%addedlist(i+1)
end do
p%nbadded=p%nbadded+nbadded
!$OMP END SINGLE
!$OMP BARRIER
if(.not. p%is_field) then
!$OMP MASTER
deallocate(p%addedlist)
!$OMP END MASTER
return
end if
! For each set of added particles
Do i=1,n,2
nptotinit=p%addedlist(i)
nbadded=p%addedlist(i+1)
nptotend=nptotinit+nbadded-1
! Potential energy\
!$OMP DO
DO j=nptotinit,nptotend,1
addedE=addedE+p%q*p%weight*p%pot(j)*phinorm
! Kinetic energy
IF(.not. nlclassical) THEN
addedE=addedE+p%m*p%weight*vlight**2*(0.5*(p%Gamma(j)+p%Gammaold(j))-1)
ELSE
addedE=addedE+0.5*p%m*p%weight*vlight**2*(p%U(1,j)*p%Uold(1,j) &
& +p%U(3,j)*p%Uold(3,j) &
& +p%U(2,j)*p%Uold(2,j))
END IF
END DO
!$OMP END DO
end do
!$OMP ATOMIC UPDATE
loc_etot0=loc_etot0+addedE
!$OMP END ATOMIC
!$OMP BARRIER
!$OMP MASTER
deallocate(p%addedlist)
!$OMP END MASTER
end if
end subroutine calc_newparts_energy
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Delete particle at given index removing its energy from the diagnosed quantities
!
!> @param [in] index index of particle to be deleted
!---------------------------------------------------------------------------
SUBROUTINE delete_part(p, index, replace)
!! This will destroy particle at the given index
USE constants, ONLY: vlight
USE bsplines
USE geometry
USE basic, ONLY: phinorm, nlclassical
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(IN) :: index
LOGICAL, OPTIONAL :: replace
LOGICAL:: repl
IF(present(replace)) THEN
repl=replace
ELSE
repl=.true.
END IF
!Computes the potential at the particle position with phi_ext+phi_s
IF(index .le. p%Nploc) THEN
IF(p%is_field) THEN
loc_etot0=loc_etot0-p%q*p%weight*(p%pot(index))*phinorm
IF(.not. nlclassical) THEN
loc_etot0=loc_etot0-p%m*p%weight*vlight**2*(p%Gamma(index)-1)
ELSE
loc_etot0=loc_etot0-0.5*p%m*p%weight*vlight**2*(p%U(1,index)**2+p%U(3,index)**2+p%U(2,index)**2)
END IF
END IF
IF(repl) THEN
! We fill the gap
CALL move_part(p, p%Nploc, index)
p%partindex(p%Nploc)=-1
! Reduce the total number of simulated parts
p%Nploc=p%Nploc-1
END IF
END IF
END SUBROUTINE delete_part
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Loads a uniform density of particles on a rectangular annulus qith maxwellian velocities
!
!> @param [inout] p particle memory to load into
!> @param [inout] VR array of radial velocity for the particles
!> @param [inout] VTHET array of azimuthal velocity for the particles
!> @param [inout] VZ array of axial velocity for the particles
!---------------------------------------------------------------------------
SUBROUTINE loaduniformRZ(p, VR,VZ,VTHET)
USE basic, ONLY: plasmadim, rnorm, temp, qsim, msim
USE constants, ONLY: me, kb, elchar
REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:)
TYPE(particles), INTENT(INOUT):: p
CALL creat_parts(p, size(VR,1))
p%Nploc=size(VR,1)
p%Nptot=size(VR,1)
p%q=sign(elchar,qsim)
p%weight=msim/me
p%m=me
p%qmRatio=qsim/msim
! Initial distribution in z with normalisation
CALL loduni(1,p%pos(3,1:p%Nploc))
p%pos(3,1:p%Nploc)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*p%pos(3,1:p%Nploc))/rnorm
! Initial distribution in r with normalisation
CALL lodlinr(2,p%pos(1,1:p%Nploc),plasmadim(3),plasmadim(4))
p%pos(1,1:p%Nploc)=p%pos(1,1:p%Nploc)/rnorm
! Initial velocities distribution
CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp)
END SUBROUTINE loaduniformRZ
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Loads a cloud of electrons trapped in a magnetic mirror according to Davidsons equilibrium
!> p117 of physics of non-neutral plasma book
!
!> @param [inout] p particle memory to load into
!> @param [inout] VR array of radial velocity for the particles
!> @param [inout] VTHET array of azimuthal velocity for the particles
!> @param [inout] VZ array of axial velocity for the particles
!---------------------------------------------------------------------------
SUBROUTINE loadDavidson(p, VR,VZ,VTHET, lodr)
USE constants, ONLY: me, kb, elchar
USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, width, qsim, msim, &
& omegac, zgrid, nz, rnorm, n0, nblock, temp
procedure(rloader)::lodr
TYPE(particles), INTENT(INOUT):: p
REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:)
REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb, z
REAL(kind=db) :: r0, deltar2, halfLz, Mirrorratio, Le, VOL
INTEGER :: j, n, blockstart, blockend, addedpart, remainparts
INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize
CALL creat_parts(p, size(VR,1))
p%Nploc=size(VR,1)
p%Nptot=p%Nploc
Allocate(ra(nblock),rb(nblock), z(0:nblock))
!r0=(plasmadim(4)+plasmadim(3))/2
r0=sqrt(4*H0/(me*omegac**2))
halfLz=(zgrid(nz)+zgrid(0))/2
MirrorRatio=(Rcurv-1)/(Rcurv+1)
z(0)=plasmadim(1)
DO n=1,nblock
! Compute limits in radius and load radii for each part
Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1)
z(n)=z(0)+n*(plasmadim(2)-plasmadim(1))/nblock
deltar2=1-MirrorRatio*cos(2*pi*Le/width)
rb(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2))
ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/H0*deltar2))
END DO
VOL=SUM(2*pi*MINVAL(ra)*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock)
qsim=VOL*n0*elchar/nplasma
msim=abs(qsim)/elchar*me
p%weight=abs(qsim)/elchar
p%m=me
p%q=sign(elchar,qsim)
p%qmRatio=p%q/p%m
blockstart=1
blockend=0
ALLOCATE(blocksize(nblock))
WRITE(*,*) "blocksize: ", size(blocksize), nblock
DO n=1,nblock
blocksize(n)=nplasma/VOL*2*pi*MINVAL(ra)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock
END DO
remainparts=p%Nploc-SUM(blocksize)
addedpart=1
n=nblock/2
j=1
DO WHILE(remainparts .GT. 0)
blocksize(n)=blocksize(n)+addedpart
remainparts=remainparts-addedpart
n=n+j
j=-1*(j+SIGN(1,j))
END DO
CALL loadPartSlices(p, lodr, ra, rb, z, blocksize)
IF(distribtype .eq. 5) THEN
CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp)
VZ=VZ/4
VR=VR*8
VTHET=VTHET*8
ELSE
Call loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0)
END IF
END SUBROUTINE loadDavidson
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Computes the velocities for a cloud of electrons trapped in a magnetic mirror according to Davidsons equilibrium
!> p117 of physics of non-neutral plasma book. This equilibrium assume mono energy and mono canonical angular momentum
!
!> @param [inout] p particle memory to load into
!> @param [inout] VR array of radial velocity for the particles
!> @param [inout] VTHET array of azimuthal velocity for the particles
!> @param [inout] VZ array of axial velocity for the particles
!> @param [in] H0 Total energy of each particle
!> @param [in] P0 Initial canonical angular momentum of each particle
!---------------------------------------------------------------------------
SUBROUTINE loadDavidsonVelocities(p, VR,VZ,VTHET, H0, P0)
USE constants, ONLY: me, kb, elchar
USE basic, ONLY: rnorm, Rcurv, B0, width, vnorm, zgrid, nz
TYPE(particles), INTENT(INOUT):: p
REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:)
REAL(kind=db), INTENT(IN):: H0, P0
REAL(kind=db) :: athetpos, rg, zg, halfLz, Mirrorratio, Pcomp, Acomp
INTEGER :: i
MirrorRatio=(Rcurv-1)/(Rcurv+1)
halfLz=(zgrid(nz)+zgrid(0))/2
! Load velocities theta velocity
! Loading of r and z velocity is done in adapt_vinit to have
! access to parts%pot
DO i=1,p%Nploc
! Interpolation for Magnetic potential
rg=p%pos(1,i)*rnorm
zg=(p%pos(3,i)-halfLz)*rnorm
Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width))
Pcomp=P0/rg/p%m
Acomp=-p%qmRatio*Athetpos
VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/p%m)),Pcomp+Acomp)
!VTHET(i)=Pcomp+Acomp
END DO
VTHET=VTHET/vnorm
VZ=0._db
VR=0._db
p%Davidson=.true.
p%H0=H0
p%P0=P0
END SUBROUTINE loadDavidsonvelocities
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Computes the particles velocities according to a maxwellian distribution of temperature temperature [K]
!
!> @param [inout] p particle memory to load into
!> @param [inout] VR array of radial velocity for the particles
!> @param [inout] VTHET array of azimuthal velocity for the particles
!> @param [inout] VZ array of axial velocity for the particles
!> @param [in] temperature temperature in [k] of the distribution function
!---------------------------------------------------------------------------
SUBROUTINE loadGaussianVelocities(p, VR,VZ,VTHET, temperature)
USE basic, ONLY: vnorm
USE constants, ONLY: kb
REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:)
TYPE(particles), INTENT(INOUT):: p
REAL(kind=db), INTENT(IN):: temperature
REAL(kind=db):: vth
! Initial velocities distribution
vth=sqrt(2.0/3.0*kb*temperature/p%m)/vnorm !thermal velocity
CALL lodgaus(3,VZ)
CALL lodgaus(5,VR)
CALL lodgaus(7,VTHET)
VZ=VZ*vth
VR=VR*vth
VTHET=VTHET*vth
p%temperature=temperature
p%Davidson=.false.
END SUBROUTINE loadGaussianVelocities
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Computes the particles velocities with a uniform distribution centered in meanv and limited by meanv+spanv and meanv-spanv
!
!> @param [inout] p particle memory to load into
!> @param [inout] VR array of radial velocity for the particles
!> @param [inout] VTHET array of azimuthal velocity for the particles
!> @param [inout] VZ array of axial velocity for the particles
!> @param [in] meanv mean velocity in each direction [m/s]
!> @param [in] spanv extent of the velocity in each direction above and below the mean velocity [m/s]
!---------------------------------------------------------------------------
SUBROUTINE loadFlatTopVelocities(p, VR,VZ,VTHET, meanv, spanv)
USE basic, ONLY: vnorm
USE constants, ONLY: kb
REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:)
TYPE(particles), INTENT(INOUT):: p
REAL(kind=db), INTENT(INOUT):: meanv(3), spanv(3)
! Initial velocities distribution
meanv=meanv/vnorm !thermal velocity
spanv=spanv/vnorm
CALL loduni(3,VZ)
CALL loduni(5,VR)
CALL loduni(7,VTHET)
VR=(VR*2-1)*spanv(1)+meanv(1)
VTHET=(VTHET*2-1)*spanv(2)+meanv(2)
VZ=(VZ*2-1)*spanv(3)+meanv(3)
p%Davidson=.false.
END SUBROUTINE loadFlatTopVelocities
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Load slices of particles defined by axial and radial limits
!
!> @param [inout] p particle memory to load into
!> @param [in] lodr sampling function definig the particle distribution in r
!> @param [in] ra lower radial limit of the slice
!> @param [in] rb upper radial limit of the slice
!> @param [in] z array giving the axial limits of each slice (slice i is betwwen z(i-1) and z(i))
!> @param [in] blocksize array containing the number of particles for each slice
!---------------------------------------------------------------------------
SUBROUTINE loadPartslices(p, lodr, ra, rb, z, blocksize)
USE constants
USE basic, ONLY: rnorm
TYPE(particles), INTENT(INOUT):: p
REAL(kind=db), INTENT(IN)::ra(:), rb(:), z(0:)
INTEGER, DIMENSION(:), INTENT(IN) :: blocksize
procedure(rloader)::lodr
INTEGER :: n, i, blockstart, blockend, nblock
nblock=size(blocksize,1)
blockstart=1
blockend=0
DO n=1,nblock
blockstart=blockend+1
blockend=MIN(blockstart+blocksize(n)-1,p%Nploc)
! Initial distribution in r, theta, z
CALL lodr(2,p%pos(1,blockstart:blockend), ra(n), rb(n))
p%pos(1,blockstart:blockend)=p%pos(1,blockstart:blockend)/rnorm
DO i=blockstart,blockend,2
CALL loduni(1,p%pos(2,i:i+1))
END DO
p%pos(2,blockstart:blockend)=p%pos(2,blockstart:blockend)*2*pi
CALL loduni(1,p%pos(3,blockstart:blockend))
p%pos(3,blockstart:blockend)=(z(n-1)+p%pos(3,blockstart:blockend)*(z(n)-z(n-1)))/rnorm
END DO
END SUBROUTINE loadPartslices
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Read a particle file format to load a simulated specie in the simulation
!
!---------------------------------------------------------------------------
SUBROUTINE read_part_file(p, partfilename, VR, VZ, VTHET)
USE basic, ONLY: lu_partfile, rnorm, vnorm
implicit None
TYPE(particles), INTENT(INOUT):: p
REAL(kind=db), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VR, VZ, VTHET
CHARACTER(len=*)::partfilename
INTEGER:: nblock = 0
REAL(kind=db), Dimension(:), ALLOCATABLE:: ra, rb, z
INTEGER, Dimension(:), ALLOCATABLE:: npartsslice
INTEGER:: velocitytype=1 !< 1) gaussian with temp 2) Davidson with H0, P0
INTEGER:: radialtype=1 !< 1) 1/R 2) uniform 3) 1/R^2 4) gauss
INTEGER:: npartsalloc !< initial size of particles arrays
INTEGER:: see_id !< index of species to add particles to for SEE
INTEGER:: iiee_id !< index of species to add particles to for IIEE
INTEGER:: neuttype_id !< index of neutral gas producing ions
INTEGER:: material_id !< index determining the type of electrode material
LOGICAL:: zero_vel !< logical to chose wether or not el. are gen. with non 0 init. vel.
REAL(kind=db):: mass=me
REAL(kind=db):: charge=-elchar
REAL(kind=db):: weight=1.0
REAL(kind=db):: qmratioscale
REAL(kind=db):: meanv(3) !< mean velocity in each direction for velocitytype 3 [m/s]
REAL(kind=db):: spanv(3) !< pos/neg extent of velocity in each direction for velocitytype 3 [m/s]
CHARACTER(len=256) :: header=' ' !< header of csv file section
REAL(kind=db):: H0=3.2e-14 !< Total energy [J]
REAL(kind=db):: P0=8.66e-25 !< Canonical angular momentum [kg m^2/s]
REAL(kind=db):: temperature=10000 !< temperature in kelvins
real(kind=db):: n0 !< density factor
LOGICAL :: is_test !< Defines if particle are saved on ittracer or not
LOGICAL :: is_field !< Defines if particle contributes to Poisson solver
LOGICAL :: calc_moments !< Defines if moments matrix must be calculated each it2d
CHARACTER(len=16) :: partformat = 'slices'
INTEGER:: i, ierr, openerr
NAMELIST /partsload/ nblock, mass, charge, weight, npartsalloc, velocitytype, &
& radialtype, temperature, H0, P0, is_test, n0, partformat, meanv, spanv, &
& calc_moments, qmratioscale, is_field, see_id, iiee_id, neuttype_id, material_id, zero_vel
! Set defaults
qmratioscale=1.0
weight=1.0
meanv=0
spanv=0
mass=me
charge=-elchar
calc_moments=.false.
is_test=.false.
is_field=.true.
see_id=-1
iiee_id=-1
neuttype_id=1
material_id=1
zero_vel=.false.
! Open the paticle file
OPEN(UNIT=lu_partfile,FILE=trim(partfilename),ACTION='READ',IOSTAT=openerr)
header=' '
IF(openerr .ne. 0) THEN
CLOSE(unit=lu_partfile)
RETURN
END IF
READ(lu_partfile,partsload)
IF(mpirank .eq.0) THEN
WRITE(*,'(a,a)')"reading partfile: ", trim(partfilename)
WRITE(*,partsload)
END IF
! The plasma cloud is defined as a set of slices
IF(trim(partformat).eq.'slices') THEN
IF( nblock .ge. 1) THEN
ALLOCATE(z(0:nblock),ra(nblock),rb(nblock), npartsslice(nblock))
DO WHILE(header(1:8) .ne. '//slices')
READ(lu_partfile,'(a)') header
END DO
DO i=1,nblock
READ(lu_partfile,*) z(i-1),ra(i),rb(i),npartsslice(i)
END DO
READ(lu_partfile,*) z(nblock)
CALL creat_parts(p,max(npartsalloc,sum(npartsslice)))
p%Nploc=sum(npartsslice)
p%Nptot=p%Nploc
IF( allocated(VR) ) THEN
DEALLOCATE(VR,VZ,VTHET)
end if
if(.not. allocated(VR)) THEN
ALLOCATE(VR(p%Nploc))
ALLOCATE(VZ(p%Nploc))
ALLOCATE(VTHET(p%Nploc))
END IF
p%m=mass
p%q=charge
p%weight=weight
p%qmRatio=charge/mass*qmratioscale
p%is_test=is_test
p%is_field=is_field
p%calc_moments=calc_moments
p%Newindex=sum(npartsslice)
p%see_id = see_id
p%iiee_id = iiee_id
p%neuttype_id = neuttype_id
p%material_id = material_id
p%zero_vel = zero_vel
SELECT CASE(radialtype)
CASE(1) ! 1/R distribution in R
CALL loadPartslices(p, lodunir, ra, rb, z, npartsslice)
CASE(2) ! flat top distribution in R
CALL loadPartslices(p, lodlinr, ra, rb, z, npartsslice)
CASE(3) ! 1/R^2 distribution in R
CALL loadPartslices(p, lodinvr, ra, rb, z, npartsslice)
CASE(4) ! gaussian distribution in R
CALL loadPartslices(p, lodgausr, ra, rb, z, npartsslice)
CASE DEFAULT
IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END SELECT
SELECT CASE(velocitytype)
CASE(1) ! Gaussian with temperature
CALL loadGaussianVelocities(p, VR, VZ, VTHET, temperature)
CASE(2) ! Davidson magnetic mirror high wr equilibrium
CALL loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0)
CASE(3) ! flat top velocity
CALL loadFlatTopVelocities(p, VR, VZ, VTHET, meanv, spanv)
CASE DEFAULT
IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of velocity distribution:", velocitytype
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END SELECT
END IF
END IF
! The plasma cloud is defined as a set individual particles
IF( trim(partformat) .eq. 'parts' ) THEN
IF( nblock .ge. 1) THEN
!Allocate necessary memory
CALL creat_parts(p,max(npartsalloc,nblock))
IF( allocated(VR) ) THEN
DEALLOCATE(VR,VZ,VTHET)
end if
if(.not. allocated(VR)) THEN
ALLOCATE(VR(nblock))
ALLOCATE(VZ(nblock))
ALLOCATE(VTHET(nblock))
END IF
! Read the particles from the file
DO WHILE(header(1:8) .ne. '//parts')
READ(lu_partfile,'(a)') header
END DO
DO i=1,nblock
READ(lu_partfile,*) p%pos(1,i),p%pos(2,i),p%pos(3,i), VR(i), VTHET(i), VZ(i)
END DO
p%Nploc=nblock
p%Nptot=p%Nploc
p%m=mass
p%q=charge
p%Newindex=nblock
p%weight=weight
p%qmRatio=charge/mass*qmratioscale
p%is_test=is_test
p%is_field=is_field
p%calc_moments=calc_moments
p%see_id=see_id
p%iiee_id=iiee_id
p%neuttype_id=neuttype_id
p%material_id=material_id
p%zero_vel=zero_vel
!normalizations
p%pos(1,:)=p%pos(1,:)/rnorm
p%pos(3,:)=p%pos(3,:)/rnorm
!p%z=p%z/rnorm
VR=VR/vnorm
VTHET=VTHET/vnorm
VZ=VZ/vnorm
END IF
END IF
CLOSE(unit=lu_partfile)
END SUBROUTINE
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Increase the number of macroparticles by separating each previous macroparticles into
!> samplefactor new macroparticles of equally divided weight. The new sub particles are distributed
!> uniformly in space to maintain the density and other moments.
!
!> @param [in] samplefactor multiplicator of the number of macroparticles.
!> @param [in] p particles type to increase.
!---------------------------------------------------------------------------
SUBROUTINE upsample(p, samplefactor)
USE basic, ONLY : nplasma, dr, dz
INTEGER, INTENT(IN) ::samplefactor
TYPE(particles), INTENT(INOUT):: p
INTEGER:: i, j, currentindex
REAL(kind=db), DIMENSION(p%Nploc) :: spreaddir ! random direction for the spread of each initial macro particle
REAL(kind=db) :: dir ! Direction in which the particle is moved
REAL(kind=db) :: dl ! Particle displacement used for
! Load and scale the direction angle for spreading the new particles
CALL loduni(2, spreaddir)
spreaddir=spreaddir*2*pi/samplefactor
dl=min(minval(dz,1,dz.GT.0),minval(dr,1,dr.GT.0))/10
DO i=1,p%Nploc
DO j=1,samplefactor-1
currentindex=p%Nploc+(i-1)*(samplefactor-1)+j
CALL move_part(p,i,currentindex)
p%partindex(currentindex)=currentindex
dir = spreaddir(i)+2*pi*j/samplefactor
p%pos(1,currentindex)=p%pos(1,currentindex) + dl*cos(dir)
p%pos(3,currentindex)=p%pos(3,currentindex) + dl*sin(dir)
END DO
p%partindex(i)=i
p%pos(1,i)=p%pos(1,i) + dl*cos(spreaddir(i))
p%pos(3,i)=p%pos(3,i) + dl*sin(spreaddir(i))
END DO
nplasma=nplasma*samplefactor
p%weight=p%weight/samplefactor
p%Nploc=p%Nploc*samplefactor
p%Nptot=p%Nptot*samplefactor
END SUBROUTINE upsample
! Taken from https://rosettacode.org/wiki/Sorting_algorithms/Radix_sort#Fortran
! No Copyright is exerted due to considerable prior art in the Public Domain.
! This Fortran version by Peter Kelly ~ peter.kelly@acm.org
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject to
! the following conditions:
! The above copyright notice and this permission notice shall be
! included in all copies or substantial portions of the Software.
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!
! Implementation of a classic Radix Sort LSD style :)
SUBROUTINE LSDRADIXSORT(A , N)
IMPLICIT NONE
!
! Dummy arguments
!
INTEGER :: N
INTEGER , target, DIMENSION(0:N - 1) :: A ! All arrays based off zero, one day I'll fix it
INTENT (IN) N
INTENT (INOUT) A
!
! Local variables
!
INTEGER , DIMENSION(0:9) :: counts
INTEGER :: digitplace
INTEGER :: i
INTEGER :: j
INTEGER :: largestnum
INTEGER, DIMENSION(0:N - 1) :: results
!
digitplace = 1 ! Count of the keys
largestnum = MAXVAL(A)
DO WHILE ( (largestnum/digitplace)>0 )
counts = 0 ! Init the count array
DO i = 0 , N - 1 , 1
J = (A(i)/digitplace)
J = MODULO(j , 10)
counts(j) = counts(j) + 1
END DO
! Change count(i) so that count(i) now contains actual position of this digit in result()
! Working similar to the counting sort algorithm
DO i = 1 , 9 , 1
counts(i) = counts(i) + counts(i - 1) ! Build up the prefix sum
END DO
!
DO i = N - 1 , 0 , -1 ! Move from left to right
j = (A(i)/digitplace)
j = MODULO(j, 10)
results(counts(j) - 1) = A(i) ! Need to subtract one as we are zero based but prefix sum is 1 based
counts(j) = counts(j) - 1
END DO
!
DO i = 0 , N - 1 , 1 ! Copy the semi-sorted data into the input
A(i) = results(i)
END DO
!
digitplace = digitplace*10
END DO ! While loop
RETURN
END SUBROUTINE LSDRADIXSORT
END MODULE beam
c4science · Help