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

fennecs/src/beam_mod.f90
6af4dc7fb7b6geom_testpardiso

Last Change
Raw File

beam_mod.f90
View Options

<<<<<<< HEAD
!------------------------------------------------------------------------------
! 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.
!------------------------------------------------------------------------------
MODULE beam
!
USE constants
USE mpi
USE mpihelper
USE basic, ONLY: mpirank, mpisize
USE distrib
IMPLICIT NONE
!> Stores the particles properties for the run.
TYPE particles
INTEGER :: Nptot !< Local number of simulated particles
REAL(kind=db) :: m !< Macroparticle mass
REAL(kind=db) :: q !< Macroparticle charge
REAL(kind=db) :: qmRatio !< Charge over mass ratio
INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex !< Index in the electric potential grid for the R direction
INTEGER, DIMENSION(:), ALLOCATABLE :: Zindex !< Index in the electric potential grid for the Z direction
INTEGER, DIMENSION(:), ALLOCATABLE :: partindex !< Index of the particle to be able to follow it when it goes from one MPI host to the other
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: R !< radial coordinates of the particles
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Z !< longitudinal coordinates of the particles
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: THET !< azimuthal coordinates of the particles
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: WZ !< axial radial relative distances to the left grid line
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: WR !< radial relative distances to the bottom grid line
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: pot !< Electric potential
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Er !< Radial Electric field
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Ez !< Axial electric field
REAL(kind=db), DIMENSION(:),POINTER:: UR !< normalized radial velocity at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: URold !< normalized radial velocity at the previous time step
REAL(kind=db), DIMENSION(:),POINTER:: UTHET !< normalized azimuthal velocity at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: UTHETold !< normalized azimuthal velocity at the previous time step
REAL(kind=db), DIMENSION(:),POINTER:: UZ !< normalized axial velocity at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: UZold !< normalized axial velocity at the previous time step
REAL(kind=db), DIMENSION(:),POINTER:: Gamma !< Lorentz factor at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: Gammaold !< Lorentz factor at the previous time step
END TYPE particles
!
TYPE(particles) :: parts !< Storage for all the particles
SAVE :: parts
TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffers to send and receive particle from left and right processes
! Diagnostics (scalars)
REAL(kind=db) :: ekin=0 !< Total kinetic energz (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)
INTEGER, DIMENSION(3), SAVE :: ireducerequest=MPI_REQUEST_NULL !< MPI requests used for the IREDUCE in the diagnostic routine
!
INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process
LOGICAL:: collected=.false. !< Stores if the particles data have been collected to MPI root process during this timestep
!
CONTAINS
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Allocate the memory for the particles variable storing the particles quantities.
!
!> @param[inout] p the particles variable needing to be allocated.
!> @param[in] nparts the maximum number of particles that will be stored in this variable
!---------------------------------------------------------------------------
SUBROUTINE creat_parts(p, nparts)
TYPE(particles) :: p
INTEGER, INTENT(in) :: nparts
IF (.NOT. ALLOCATED(p%Z) ) THEN
p%Nptot = nparts
ALLOCATE(p%Z(nparts))
ALLOCATE(p%R(nparts))
ALLOCATE(p%THET(nparts))
ALLOCATE(p%WZ(nparts))
ALLOCATE(p%WR(nparts))
ALLOCATE(p%UR(nparts))
ALLOCATE(p%UZ(nparts))
ALLOCATE(p%UTHET(nparts))
ALLOCATE(p%URold(nparts))
ALLOCATE(p%UZold(nparts))
ALLOCATE(p%UTHETold(nparts))
ALLOCATE(p%Gamma(nparts))
ALLOCATE(p%Rindex(nparts))
ALLOCATE(p%Zindex(nparts))
ALLOCATE(p%partindex(nparts))
ALLOCATE(p%pot(nparts))
ALLOCATE(p%Er(nparts))
ALLOCATE(p%Ez(nparts))
ALLOCATE(p%GAMMAold(nparts))
p%partindex=-1
p%URold=0
p%UZold=0
p%UTHETold=0
p%rindex=0
p%zindex=0
p%WR=0
p%WZ=0
p%UR=0
p%UZ=0
p%UTHET=0
p%Z=0
p%R=0
p%THET=0
p%Gamma=1
p%Er=0
p%Ez=0
p%pot=0
p%gammaold=1
END IF
END SUBROUTINE creat_parts
!---------------------------------------------------------------------------
!> @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, mpisize, nlclassical
USE mpi
INTEGER:: i
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma))
parts%Nptot=nplasma
! Select case to define the type of distribution
SELECT CASE(distribtype)
CASE(1) ! Gaussian distribution in V and uniform in R
CALL loaduniformRZ(VR, VZ, VTHET)
CASE(2) !Stable distribution from Davidson 4.95 p.119
CALL loadDavidson(VR, VZ, VTHET, lodunir)
CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2
CALL loadDavidson(VR, VZ, VTHET, lodinvr)
CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R
CALL loadDavidson(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(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(VR, VZ, VTHET)
VZ = VZ/50
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
parts%partindex(i)=i
END DO
IF(nlclassical) THEN
parts%Gamma(1:nplasma)=1
ELSE
parts%Gamma(1:nplasma)=sqrt(1/(1-VR**2-VZ**2-VTHET**2))
END IF
! Normalization of the velocities
parts%UR(1:nplasma)=parts%Gamma(1:nplasma)*VR
parts%UZ(1:nplasma)=parts%Gamma(1:nplasma)*VZ
parts%UTHET(1:nplasma)=parts%Gamma(1:nplasma)*VTHET
DEALLOCATE(VZ, VR, VTHET)
CALL localisation(parts)
IF(mpisize .gt. 1) THEN
CALL distribpartsonprocs
END IF
END SUBROUTINE load_parts
!---------------------------------------------------------------------------
! 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.
!
!> @author Guillaume Le Bars EPFL/SPC
!---------------------------------------------------------------------------
SUBROUTINE bound(p)
USE basic, ONLY: zgrid, rnorm, nz, Zbounds, cstep, mpirank, partperiodic, step, leftproc, rightproc
type(particles), INTENT(INOUT):: p
INTEGER :: i, rsendnbparts, lsendnbparts, nblostparts
INTEGER, DIMENSION(p%Nptot) :: sendhole
INTEGER, DIMENSION(p%Nptot) :: losthole
LOGICAL:: leftcomm, rightcomm
rsendnbparts=0
lsendnbparts=0
nblostparts=0
losthole=0
sendhole=0
! We communicate with the left processus
leftcomm = leftproc .ne. -1 .or. partperiodic
! We communicate with the right processus
rightcomm = rightproc .ne. -1 .or. partperiodic
IF (p%Nptot .NE. 0) THEN
! Boundary condition at z direction
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,p%Nptot
IF(step .ne. 0) THEN
! If the particle is to the right of the local simulation space, it is sent to the right MPI process
IF (p%Z(i) .ge. zgrid(Zbounds(mpirank+1))) THEN
!$OMP CRITICAL (nbparts)
IF(rightcomm) THEN
rsendnbparts=rsendnbparts+1
sendhole(lsendnbparts+rsendnbparts)=i
ELSE
!WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in z"
!WRITE(*,*) "!!!!!!!!!! Particle removed !!!!!!!!!!"
nblostparts=nblostparts+1
losthole(nblostparts)=i
END IF
!$OMP END CRITICAL (nbparts)
! If the particle is to the left of the local simulation space, it is sent to the left MPI process
ELSE IF (p%Z(i) .lt. zgrid(Zbounds(mpirank))) THEN
!$OMP CRITICAL (nbparts)
IF(leftcomm) THEN
lsendnbparts=lsendnbparts+1
sendhole(lsendnbparts+rsendnbparts)=-i
ELSE
!WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in z"
!WRITE(*,*) "!!!!!!!!!! Particle removed !!!!!!!!!!"
nblostparts=nblostparts+1
losthole(nblostparts)=i
END IF
!$OMP END CRITICAL (nbparts)
END IF
END IF
! The periodic boundary conditions in z are here applied
DO WHILE (p%Z(i) .GT. zgrid(nz))
p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0)
END DO
DO WHILE (p%Z(i) .LT. zgrid(0))
p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0)
END DO
END DO
!$OMP END PARALLEL DO
END IF
IF(mpisize .gt. 1 .and. step .ne. 0) THEN
! We send the particles leaving the local simulation space to the closest neighbour
CALL particlescommunication(p, lsendnbparts, rsendnbparts, sendhole)
END IF
! If the boundary conditions are not periodic, we delete the corresponding particles
IF(nblostparts .gt. 0) THEN
DO i=nblostparts,1,-1
!WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in z: ", p%Z(losthole(i))*rnorm, zgrid(0)*rnorm, zgrid(nz)*rnorm
CALL delete_part(p, losthole(i))
END DO
WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank
END IF
END subroutine bound
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Compute the grid cell indices for each particle as well as the distance weight Wr, Wz.
!---------------------------------------------------------------------------
SUBROUTINE localisation(p)
USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm
type(particles), INTENT(INOUT):: p
INTEGER :: i, nblostparts
INTEGER, DIMENSION(p%Nptot) :: losthole
REAL(kind=db):: invdz, invdr(3)
losthole=0
nblostparts=0
invdz=1/dz
invdr=1/dr
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,p%Nptot
p%zindex(i)=(p%Z(i)-zgrid(0))*invdz
IF (p%R(i) .GT. rgrid(0) .AND. p%R(i) .LT. rgrid(nnr(1))) THEN
p%rindex(i)=(p%R(i)-rgrid(0))*invdr(1)
ELSE IF(p%R(i) .GE. rgrid(nnr(1)) .AND. p%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN
p%rindex(i)=(p%R(i)-rgrid(nnr(1)))*invdr(2)+nnr(1)
ELSE IF(p%R(i) .GE. rgrid(nnr(1)+nnr(2)) .AND. p%R(i) .LT. rgrid(nr)) THEN
p%rindex(i)=(p%R(i)-rgrid(nnr(1)+nnr(2)))*invdr(3)+nnr(1)+nnr(2)
ELSE
! If the particle is outside of the simulation space in the r direction, it is deleted.
!$OMP CRITICAL (lostparts)
nblostparts=nblostparts+1
losthole(nblostparts)=i
!$OMP END CRITICAL (lostparts)
END IF
END DO
!$OMP END PARALLEL DO
IF(nblostparts.gt.0) THEN
DO i=nblostparts,1,-1
!WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in r: ", p%R(losthole(i))*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm
CALL delete_part(p,losthole(i))
END DO
WRITE(*,'(i6.2,a,i6.2)') nblostparts, " particles lost in r on process: ", mpirank
END IF
!$OMP PARALLEL DO SIMD DEFAULT(SHARED)
DO i=1,p%Nptot
p%WZ(i)=(p%Z(i)-zgrid(p%zindex(i)))/dz;
p%WR(i)=(p%R(i)-rgrid(p%rindex(i)))/(rgrid(p%rindex(i)+1)-rgrid(p%rindex(i)));
END DO
!$OMP END PARALLEL DO SIMD
END IF
END SUBROUTINE localisation
!---------------------------------------------------------------------------
!> @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.
!
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity(p)
!
! Computes the new velocity of the particles due to Lorentz force
!
USE basic, ONLY : nlclassical
type(particles), INTENT(INOUT):: p
! Store old Velocities
CALL swappointer(p%UZold, p%UZ)
CALL swappointer(p%URold, p%UR)
CALL swappointer(p%UTHETold, p%UTHET)
CALL swappointer(p%Gammaold, p%Gamma)
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$
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity_fun(p, gamma)
!
! Computes the new velocity of the particles due to Lorentz force
!
USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz
interface
subroutine gamma(gam, UZ, UR, UTHET)
USE constants
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
REAL(kind=db), INTENT(OUT):: gam
end subroutine
end interface
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: tau
REAL(kind=db):: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2
INTEGER:: J1, J2, J3, J4
INTEGER:: i
! Normalized time increment
tau=omegac/2/omegap*dt/tnorm
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
DO i=1,p%Nptot
J1=(p%rindex(i))*(nz+1) + p%zindex(i)+1
J2=(p%rindex(i))*(nz+1) + p%zindex(i)+2
J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1
J4=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+2
! Interpolation for magnetic field
BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) &
& +p%WZ(i)*(1-p%WR(i))*Bz(J3) &
& +(1-p%WZ(i))*p%WR(i)*Bz(J2) &
& +p%WZ(i)*p%WR(i)*Bz(J1)
BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) &
& +p%WZ(i)*(1-p%WR(i))*Br(J3) &
& +(1-p%WZ(i))*p%WR(i)*Br(J2) &
& +p%WZ(i)*p%WR(i)*Br(J1)
! First half of electric pulse
p%UZ(i)=p%UZold(i)+p%Ez(i)*tau
p%UR(i)=p%URold(i)+p%ER(i)*tau
CALL gamma(p%Gamma(i), p%UZ(i), p%UR(i), p%UTHETold(i))
! Rotation along magnetic field
ZBZ=tau*BRZ/p%Gamma(i)
ZBR=tau*BRR/p%Gamma(i)
ZPZ=p%UZ(i)-ZBR*p%UTHETold(i) !u'_{z}
ZPR=p%UR(i)+ZBZ*p%UTHETold(i) !u'_{r}
ZPTHET=p%UTHETold(i)+(ZBR*p%UZ(i)-ZBZ*p%UR(i)) !u'_{theta}
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZBZ2=2*ZBZ/SQR
ZBR2=2*ZBR/SQR
p%UZ(i)=p%UZ(i)-ZBR2*ZPTHET !u+_{z}
p%UR(i)=p%UR(i)+ZBZ2*ZPTHET !u+_{r}
p%UTHET(i)=p%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta}
! Second half of acceleration
p%UZ(i)=p%UZ(i)+p%EZ(i)*tau
p%UR(i)=p%UR(i)+p%ER(i)*tau
! Final computation of the Lorentz factor
CALL gamma(p%Gamma(i), p%UZ(i), p%UR(i), p%UTHETold(i))
END DO
!$OMP END PARALLEL DO SIMD
END IF
collected=.false.
END SUBROUTINE comp_velocity_fun
!---------------------------------------------------------------------------
!> @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
!---------------------------------------------------------------------------
ELEMENTAL SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET)
!$OMP declare simd(gamma_classical)
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
REAL(kind=db), INTENT(OUT):: gamma
gamma=1.0
END SUBROUTINE 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
!---------------------------------------------------------------------------
ELEMENTAL SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET)
!$OMP declare simd(gamma_relativistic)
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
REAL(kind=db), INTENT(OUT):: gamma
gamma=sqrt(1+UZ**2+UR**2+UTHET**2)
END SUBROUTINE
!---------------------------------------------------------------------------
!> @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.
!---------------------------------------------------------------------------
SUBROUTINE push(p)
Use basic, ONLY: dt, tnorm
type(particles), INTENT(INOUT):: p
REAL(kind=db):: XP, YP, COSA, SINA, U1, U2, ALPHA
INTEGER :: i
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(XP, YP, COSA, SINA, U1, U2, ALPHA)
DO i=1,p%Nptot
! Local Cartesian coordinates
XP=p%R(i)+dt/tnorm*p%UR(i)/p%Gamma(i)
YP=dt/tnorm*p%UTHET(i)/p%Gamma(i)
! Conversion to cylindrical coordiantes
p%Z(i)=p%Z(i)+dt/tnorm*p%UZ(i)/p%Gamma(i)
p%R(i)=sqrt(XP**2+YP**2)
! Computation of the rotation angle
IF (p%R(i) .EQ. 0) THEN
COSA=1
SINA=0
ALPHA=0
ELSE
COSA=XP/p%R(i)
SINA=YP/p%R(i)
ALPHA=asin(SINA)
END IF
! New azimuthal position
p%THET(i)=MOD(p%THET(i)+ALPHA,2*pi)
! Velocity in rotated reference frame
U1=COSA*p%UR(i)+SINA*p%UTHET(i)
U2=-SINA*p%UR(i)+COSA*p%UTHET(i)
p%UR(i)=U1
p%UTHET(i)=U2
END DO
!$OMP END PARALLEL DO SIMD
END IF
collected=.false.
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 also computes the first and second order
!> moments of the distribution function
!
!---------------------------------------------------------------------------
SUBROUTINE diagnostics
!
! Compute energies
!
USE constants, ONLY: vlight
USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, ierr, step, it0d, it2d, nlend, nr, itparts, nplasma, partperiodic, newres, Zbounds
INTEGER :: i, gridindex
INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE)
CALL MPI_WAITALL(3, ireducerequest, stat, ierr)
! Reset the quantities
ekin=0
epot=0
etot=0
! Computation of the kinetic and potential energy as well as fluid velocities and density
!$OMP PARALLEL DO SIMD REDUCTION(+:epot, ekin) DEFAULT(SHARED)
DO i=1,parts%Nptot
! Potential energy
epot=epot+0.5*qsim*parts%pot(i)*phinorm
! Kinetic energy
IF(.not. nlclassical) THEN
ekin=ekin+msim*vlight**2*(parts%Gamma(i)-1)
ELSE
ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))
END IF
END DO
!$OMP END PARALLEL DO SIMD
! Shift to Etot at cstep=1 (not valable yet at cstep=0!)
IF(cstep.EQ.1 .or. (newres .AND. step .EQ. 1)) THEN
! Compute the local total energy
loc_etot0 = epot+ekin
END IF
etot=loc_etot0
! The computed energy is sent to the root process
IF(mpisize .gt.1) THEN
IF(mpirank .eq.0 ) THEN
CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
CALL MPI_IREDUCE(etot, etot0, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
ELSE
CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
CALL MPI_IREDUCE(etot, etot0, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
END IF
ELSE
etot0=loc_etot0
END IF
! Wait for the communication of the energy to be finished
CALL MPI_WAITALL(3, ireducerequest(1:3), stat(1:3*MPI_STATUS_SIZE), ierr)
!CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr)
! Compute the total energy
etot=epot+ekin
! Send the local number of particles on each node to the root process
IF(modulo(step,it0d).eq. 0 .and. mpisize .gt. 1) THEN
Nplocs_all(mpirank)=parts%Nptot
IF(mpirank .eq.0 ) THEN
CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
ELSE
CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
END IF
IF(mpirank .eq. 0 .AND. partperiodic .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN
!WRITE(*,*) "Error particle lost on some processus"
!CALL MPI_abort(MPI_COMM_WORLD,-1,ierr)
END IF
END IF
!Calls the communications to send the particles data to the root process for diagnostic file every itparts time steps
IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN
CALL collectparts(parts)
END IF
end subroutine diagnostics
!---------------------------------------------------------------------------
!> @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)
USE basic, ONLY: mpirank, mpisize, ierr
type(particles), INTENT(INOUT):: p
INTEGER, DIMENSION(:), ALLOCATABLE :: displs
INTEGER:: i
IF(collected) RETURN ! exit subroutine if particles have already been collected during this time step
IF(mpirank .ne. 0) THEN
ALLOCATE(displs(0:mpisize-1))
displs=0
! Send Particles informations to root process
CALL MPI_Gatherv(p%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%Z, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%R, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%THET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%THET, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UR, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UZ, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UTHET, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%pot, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%Rindex, Nplocs_all(mpirank), MPI_INTEGER, p%Rindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%Zindex, Nplocs_all(mpirank), MPI_INTEGER, p%Zindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%partindex, Nplocs_all(mpirank), MPI_INTEGER, p%partindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
ELSE
ALLOCATE(displs(0:mpisize-1))
displs(0)=0
Do i=1,mpisize-1
displs(i)=displs(i-1)+Nplocs_all(i-1)
END DO
! Receive particle information from all processes
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%Z, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%R, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%THET, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UR, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UZ, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UTHET, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%pot, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%Rindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%Zindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%partindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
p%partindex(sum(Nplocs_all)+1:)=-1
END IF
collected=.TRUE.
END SUBROUTINE collectparts
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Computes the velocities at time t-1/2 to keep the second order precision in time on the velocity.
!
!---------------------------------------------------------------------------
SUBROUTINE adapt_vinit(p)
!! Computes the velocity at time -dt/2 from velocities computed at time 0
!
USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, &
& SQR, ZBZ2, ZBR2, Vperp, v2
INTEGER :: J1, J2, J3, J4, 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) THEN
ALLOCATE(VR(p%Nptot),VZ(p%Nptot),VTHET(p%Nptot))
CALL loduni(7,VZ)
VZ=VZ*2*pi
VTHET=p%UTHET/p%Gamma*vnorm
DO i=1,p%Nptot
Vperp=sqrt(MAX(2*H0/me+2*elchar/me*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%UR(i)=p%Gamma(i)*VR(i)/vnorm
p%UZ(i)=p%Gamma(i)*VZ(i)/vnorm
p%UTHET(i)=p%Gamma(i)*VTHET(i)/vnorm
END DO
DEALLOCATE(VR,VZ,VTHET)
END IF
! Normalised time increment
tau=-omegac/2/omegap*dt/tnorm
! Store old Velocities
CALL swappointer(p%UZold, p%UZ)
CALL swappointer(p%URold, p%UR)
CALL swappointer(p%UTHETold, p%UTHET)
CALL swappointer(p%Gammaold, p%Gamma)
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
DO i=1,p%Nptot
! Compute the particle linear indices for the magnetic field interpolation
J1=(p%rindex(i))*(nz+1)+p%zindex(i)+1
J2=J1+1
J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1
J4=J3+1
! Interpolation for magnetic field
BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) &
& +p%WZ(i)*(1-p%WR(i))*Bz(J3) &
& +(1-p%WZ(i))*p%WR(i)*Bz(J2) &
& +p%WZ(i)*p%WR(i)*Bz(J1)
BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) &
& +p%WZ(i)*(1-p%WR(i))*Br(J3) &
& +(1-p%WZ(i))*p%WR(i)*Br(J2) &
& +p%WZ(i)*p%WR(i)*Br(J1)
! Half inverse Rotation along magnetic field
ZBZ=tau*BRZ/p%Gammaold(i)
ZBR=tau*BRR/p%Gammaold(i)
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZPZ=(p%UZold(i)-ZBR*p%UTHETold(i))/SQR !u-_{z}
ZPR=(p%URold(i)+ZBZ*p%UTHETold(i))/SQR !u-_{r}
ZPTHET=p%UTHETold(i)+(ZBR*p%UZold(i)-ZBZ*p%URold(i))/SQR !u-_{theta}
p%UZ(i)=ZPZ
p%UR(i)=ZPR
p%UTHET(i)=ZPTHET
! half of decceleration
p%UZ(i)=p%UZ(i)+p%Ez(i)*tau
p%UR(i)=p%UR(i)+p%Er(i)*tau
IF(.not. nlclassical) THEN
p%Gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2)
END IF
END DO
!$OMP END PARALLEL DO
END IF
END SUBROUTINE adapt_vinit
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief In the case of MPI parallelism, computes the indices of the particle assigned to the current process.
!---------------------------------------------------------------------------
SUBROUTINE distribpartsonprocs
! 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, ierr, Zbounds
INTEGER:: k, i, nbparts
REAL(kind=db):: idealnbpartsperproc
INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k)
INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector
INTEGER:: partindexlast ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector
INTEGER, DIMENSION(0:nz):: 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, zindex
partspercol=0
idealnbpartsperproc = FLOOR(REAL(parts%Nptot)/REAL(mpisize))
Zmin=MINVAL(parts%Zindex(1:parts%Nptot))
Zmax=MAXVAL(parts%Zindex(1:parts%Nptot))
Zperproc=(Zmax-Zmin)/mpisize
IF(Zmax .eq. 0) Zmax=nz
DO k=1,parts%Nptot
zindex=parts%Zindex(k)
partspercol(zindex)=partspercol(zindex)+1
END DO
IF (Zperproc .eq. 0) THEN
!! No particles are present initially
Zperproc=nz/mpisize
Zmin=0
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
DO k=1,mpisize-1
nbparts=0
DO WHILE(nbparts<0.99*idealnbpartsperproc .and. i .lt. Zmax)
nbparts=nbparts+partspercol(i)
i=i+1
END DO
Zbounds(k)=i
END DO
END IF
partindexstart=1
totparts=0
partindexlast=parts%Nptot
!DO k=0,mpisize-2
! Nplocs_all(k)=idealnbpartsperproc
!END DO
! NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize)
!DO k=1, mpisize-1
! Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc))
!END DO
DO k=0,mpisize-1
Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1))
END DO
! Store local end of particle vector
partindexlast=SUM(Nplocs_all(0:mpirank))
! Store local start of particle vector
IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1
IF(mpirank .ne. 0 .and. Nplocs_all(mpirank) .gt. 0) THEN
DO k= partindexstart, partindexlast
CALL move_part(parts,k,k-partindexstart+1)
END DO
END IF
IF(INT(SUM(Nplocs_all)) .ne. parts%Nptot) THEN
WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " Nptot:", parts%Nptot
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
parts%Nptot=Nplocs_all(mpirank)
WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot
END SUBROUTINE distribpartsonprocs
!_______________________________________________________________________________
!---------------------------------------------------------------------------
!> @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)
!> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative
!---------------------------------------------------------------------------
SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, sendholes)
USE mpihelper, ONLY: particle_type
USE basic, ONLY: rightproc, leftproc, ierr
type(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(in) :: sendholes(:)
INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0
INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest
INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE)
INTEGER :: lsentnbparts, rsentnbparts
INTEGER :: lreceivednbparts, rreceivednbparts
lsentnbparts=lsendnbparts
rsentnbparts=rsendnbparts
sendrequest=MPI_REQUEST_NULL
recvrequest=MPI_REQUEST_NULL
! Send and receive the number of particles to exchange
CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, recvrequest(1), ierr)
CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, recvrequest(2), ierr)
CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, sendrequest(1), ierr)
CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, sendrequest(2), ierr)
CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr)
CALL MPI_Wait(recvrequest(2), recvstatus(2), ierr)
recvrequest=MPI_REQUEST_NULL
lreceivednbparts=lrecvnbparts
rreceivednbparts=rrecvnbparts
! Re/allocate enough memory to store the incoming particles
IF( ALLOCATED(rrecvpartbuff) .AND. (size(rrecvpartbuff) .lt. rrecvnbparts) ) DEALLOCATE(rrecvpartbuff)
IF(.not. ALLOCATED(rrecvpartbuff) .and. rrecvnbparts .gt. 0) ALLOCATE(rrecvpartbuff(INT(rreceivednbparts*1.3)))
IF( ALLOCATED(lrecvpartbuff) .AND. (size(lrecvpartbuff) .lt. lrecvnbparts) ) DEALLOCATE(lrecvpartbuff)
IF((.not. ALLOCATED(lrecvpartbuff) ).and. lrecvnbparts .gt. 0) ALLOCATE(lrecvpartbuff(INT(lreceivednbparts*1.3)))
! Receive particles from left and right processes to the corresponding buffers
IF ( lrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, recvrequest(1), ierr)
END IF
IF( rrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, recvrequest(2), ierr)
END IF
! Copy the leaving particles to the corresponding send buffers
IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN
CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes)
END IF
CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr)
CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr)
! Send the particles to the left and right neighbours
IF( lsendnbparts .gt. 0) THEN
CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, sendrequest(1), ierr)
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, sendrequest(2), ierr)
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
END IF
IF ( rreceivednbparts .gt. 0) THEN
CALL MPI_Wait(recvrequest(2), recvstatus(MPI_STATUS_SIZE+1), 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
END IF
! Copy the incoming particles from the receive buffers to the simulation parts variable
CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes)
! 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)
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_Wait(sendrequest(2), sendstatus(2), ierr)
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
!> @param [in] sendholes array containing the indices of the particle having left the local domain in ascending order.
!---------------------------------------------------------------------------
SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts, sendholes)
!
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts
INTEGER, INTENT(in) :: sendholes(:)
INTEGER k,partpos, partdiff
! computes if we received less particles than we sent
partdiff=sendnbparts-rrecvnbparts-lrecvnbparts
! 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
partpos=abs(sendholes(k))
ELSE
! Add at the end of parts
p%Nptot=p%Nptot+1
partpos=p%Nptot
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
partpos=abs(sendholes(k+rrecvnbparts))
ELSE
p%Nptot=p%Nptot+1
partpos=p%Nptot
END IF
CALL Insertincomingpart(p, lrecvpartbuff, k, partpos)
END DO
END IF
! If we received less particles than we sent, we fill the remaining holes with the particles from the end of the
! parts arrays
IF(partdiff .gt. 0) THEN
DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts
CALL move_part(parts, parts%Nptot, abs(sendholes(k)))
p%Nptot = p%Nptot-1
END DO
END IF
!
END SUBROUTINE Addincomingparts
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Copy one particle from the receive buffers to the local simulation variable parts.
!
!> @param [in] buffer receive buffer containing the particles parameters to copy from
!> @param [in] bufferindex particle index in the receive buffer
!> @param [in] partsindex destination particle index in the local parts variable
!---------------------------------------------------------------------------
SUBROUTINE Insertincomingpart(p, buffer, bufferindex, partsindex)
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: bufferindex, partsindex
TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer
p%partindex(partsindex) = buffer(bufferindex)%partindex
p%Rindex(partsindex) = buffer(bufferindex)%Rindex
p%Zindex(partsindex) = buffer(bufferindex)%Zindex
p%R(partsindex) = buffer(bufferindex)%R
p%Z(partsindex) = buffer(bufferindex)%Z
p%THET(partsindex) = buffer(bufferindex)%THET
p%UZ(partsindex) = buffer(bufferindex)%UZ
p%UR(partsindex) = buffer(bufferindex)%UR
p%UTHET(partsindex) = buffer(bufferindex)%UTHET
p%Gamma(partsindex) = buffer(bufferindex)%Gamma
!
!
END SUBROUTINE Insertincomingpart
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Copy one particle from the local parts variable to the send buffer.
!
!> @param [in] buffer send buffer to copy to
!> @param [in] bufferindex particle index in the send buffer
!> @param [in] partsindex origin particle index in the local parts variable
!---------------------------------------------------------------------------
SUBROUTINE Insertsentpart(p, buffer, bufferindex, partsindex)
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: bufferindex, partsindex
TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer
buffer(bufferindex)%partindex = p%partindex(partsindex)
buffer(bufferindex)%Rindex = p%Rindex(partsindex)
buffer(bufferindex)%Zindex = p%Zindex(partsindex)
buffer(bufferindex)%R = p%R(partsindex)
buffer(bufferindex)%Z = p%Z(partsindex)
buffer(bufferindex)%THET = p%THET(partsindex)
buffer(bufferindex)%UZ = p%UZ(partsindex)
buffer(bufferindex)%UR = p%UR(partsindex)
buffer(bufferindex)%UTHET = p%UTHET(partsindex)
buffer(bufferindex)%Gamma = p%Gamma(partsindex)
!
END SUBROUTINE Insertsentpart
!---------------------------------------------------------------------------
!> @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)
!> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative
!---------------------------------------------------------------------------
SUBROUTINE AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes)
!
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(in) :: sendholes(:)
INTEGER:: partpos, k
INTEGER:: lsendpos, rsendpos
lsendpos=0
rsendpos=0
! Verify if the buffers are big enough to store the outgoing particles
IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff)
IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff)
! Allocate the buffers if needed
IF(.not. ALLOCATED(rsendpartbuff) .and. rsendnbparts .gt. 0) ALLOCATE(rsendpartbuff(INT(rsendnbparts*1.3)))
IF(.not. ALLOCATED(lsendpartbuff) .and. lsendnbparts .gt. 0) ALLOCATE(lsendpartbuff(INT(lsendnbparts*1.3)))
! Loop over the outgoing particles and fill the correct send buffer
Do k=lsendnbparts+rsendnbparts,1,-1
partpos=abs(sendholes(k))
IF(sendholes(k) .GT. 0) THEN
rsendpos=rsendpos+1
CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos)
ELSE IF(sendholes(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 Exchange two particles in the parts variable.
!
!> @param [in] index1 index in parts of the first particle to exchange.
!> @param [in] index2 index in parts of the second particle to exchange.
!---------------------------------------------------------------------------
SUBROUTINE exchange_parts(p, index1, index2)
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(IN) :: index1, index2
REAL(kind=db):: R, Z, THET, UR, UZ, UTHET, Gamma
INTEGER :: Rindex, Zindex, partindex
!! Exchange particle at index1 with particle at index2
! Store part at index1 in temporary value
partindex= p%partindex(index1)
Gamma = p%Gamma(index1)
R = p%R(index1)
Z = p%Z(index1)
THET = p%THET(index1)
UR = p%UR(index1)
UTHET = p%UTHET(index1)
UZ = p%UZ(index1)
Rindex = p%Rindex(index1)
Zindex = p%Zindex(index1)
! Move part at index2 in part at index 1
p%partindex(index1)= p%partindex(index2)
p%Gamma(index1) = p%Gamma(index2)
p%R(index1) = p%R(index2)
p%Z(index1) = p%Z(index2)
p%THET(index1) = p%THET(index2)
p%UR(index1) = p%UR(index2)
p%UTHET(index1) = p%UTHET(index2)
p%UZ(index1) = p%UZ(index2)
p%Rindex(index1) = p%Rindex(index2)
p%Zindex(index1) = p%Zindex(index2)
! Move temporary values from part(index1) to part(index2)
p%partindex(index2)= partindex
p%Gamma(index2) = Gamma
p%R(index2) = R
p%Z(index2) = Z
p%THET(index2) = THET
p%UR(index2) = UR
p%UTHET(index2) = UTHET
p%UZ(index2) = UZ
p%Rindex(index2) = Rindex
p%Zindex(index2) = Zindex
END SUBROUTINE exchange_parts
SUBROUTINE change_parts_allocation(p, sizedifference)
implicit none
TYPE(particles), INTENT(INOUT):: p
INTEGER,INTENT(IN) :: sizedifference
CALL change_array_size_int(p%Rindex, sizedifference)
CALL change_array_size_int(p%Zindex, sizedifference)
CALL change_array_size_int(p%partindex, sizedifference)
CALL change_array_size_dp(p%ER,sizedifference)
CALL change_array_size_dp(p%EZ,sizedifference)
CALL change_array_size_dp(p%pot,sizedifference)
CALL change_array_size_dp(p%R,sizedifference)
CALL change_array_size_dp(p%Z,sizedifference)
CALL change_array_size_dp(p%THET,sizedifference)
CALL change_array_size_dp(p%WR,sizedifference)
CALL change_array_size_dp(p%WZ,sizedifference)
CALL change_array_size_dp_ptr(p%UR,sizedifference)
CALL change_array_size_dp_ptr(p%URold,sizedifference)
CALL change_array_size_dp_ptr(p%UZ,sizedifference)
CALL change_array_size_dp_ptr(p%UZold,sizedifference)
CALL change_array_size_dp_ptr(p%UTHET,sizedifference)
CALL change_array_size_dp_ptr(p%UTHETold,sizedifference)
CALL change_array_size_dp_ptr(p%Gamma,sizedifference)
CALL change_array_size_dp_ptr(p%Gammaold,sizedifference)
p%Nptot=MIN(p%Nptot,size(p%R)+sizedifference)
END SUBROUTINE change_parts_allocation
SUBROUTINE change_array_size_dp(arr, sizedifference)
implicit none
REAL(kind=db), ALLOCATABLE, INTENT(INOUT):: arr(:)
INTEGER, INTENT(IN):: sizedifference
REAL(kind=db), ALLOCATABLE:: temp(:)
INTEGER:: current_size, new_size
if(allocated(arr)) THEN
current_size=size(arr)
new_size=current_size+sizedifference
ALLOCATE(temp(new_size))
temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size))
DEALLOCATE(arr)
CALL move_alloc(temp, arr)
END IF
END SUBROUTINE change_array_size_dp
SUBROUTINE change_array_size_dp_ptr(arr, sizedifference)
implicit none
REAL(kind=db), POINTER, INTENT(INOUT):: arr(:)
INTEGER, INTENT(IN):: sizedifference
REAL(kind=db), POINTER:: temp(:)
INTEGER:: current_size, new_size
if(associated(arr)) THEN
current_size=size(arr)
new_size=current_size+sizedifference
ALLOCATE(temp(new_size))
temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size))
DEALLOCATE(arr)
arr=> temp
END IF
END SUBROUTINE change_array_size_dp_ptr
SUBROUTINE change_array_size_int(arr, sizedifference)
implicit none
INTEGER, ALLOCATABLE, INTENT(INOUT):: arr(:)
INTEGER, INTENT(IN):: sizedifference
INTEGER, ALLOCATABLE:: temp(:)
INTEGER:: current_size, new_size
if(allocated(arr)) THEN
current_size=size(arr)
new_size=current_size+sizedifference
ALLOCATE(temp(new_size))
temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size))
DEALLOCATE(arr)
CALL move_alloc(temp,arr)
END IF
END SUBROUTINE change_array_size_int
SUBROUTINE add_created_part(p, buffer)
USE bsplines
USE basic, ONLY: splrz, qsim, phinorm, nlclassical, msim
IMPLICIT NONE
TYPE(particles), INTENT(INOUT):: p
TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:)
INTEGER:: i, nptotinit, parts_size_increase, nb_added
nptotinit=p%Nptot+1
nb_added=size(buffer,1)
! if there is not enough space in the parts simulation buffer, increase the parst size
IF(p%Nptot + nb_added .gt. size(p%Z,1)) THEN
parts_size_increase=Max(floor(0.1*size(p%Z,1)),nb_added)
CALL change_parts_allocation(p, parts_size_increase)
END IF
DO i=1,nb_added
p%Nptot=p%Nptot+1
CALL Insertincomingpart(p, buffer,i,p%Nptot)
END DO
! Compute the energy added by these new particles for energy conservation diagnostic
IF(nb_added .gt. 0) THEN
CALL getgrad(splrz, p%Z(nptotinit:p%Nptot), p%R(nptotinit:p%Nptot), &
& p%pot(nptotinit:p%Nptot), p%EZ(nptotinit:p%Nptot), p%ER(nptotinit:p%Nptot))
p%EZ(nptotinit:p%Nptot)=-p%Ez(nptotinit:p%Nptot)
p%ER(nptotinit:p%Nptot)=-p%ER(nptotinit:p%Nptot)
loc_etot0=loc_etot0+qsim*sum(p%pot(nptotinit:p%Nptot))*phinorm
IF(.not. nlclassical) THEN
loc_etot0=loc_etot0+sum(msim*vlight**2*(p%Gamma(nptotinit:p%Nptot)-1))
ELSE
loc_etot0=loc_etot0+0.5*msim*vlight**2*sum(p%UR(nptotinit:p%Nptot)**2 &
& +p%UZ(nptotinit:p%Nptot)**2 &
& +p%UTHET(nptotinit:p%Nptot)**2)
END IF
END IF
END SUBROUTINE
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Delete particle at given index removing its energy from the system
!
!> @param [in] index index of particle to be deleted
!---------------------------------------------------------------------------
SUBROUTINE delete_part(p, index)
!! This will destroy particle at index
USE constants, ONLY: vlight
USE basic, ONLY: qsim, phinorm, msim, nlclassical
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(IN) :: index
IF(index .le. p%Nptot) THEN
loc_etot0=loc_etot0-qsim*p%pot(index)*phinorm
IF(.not. nlclassical) THEN
loc_etot0=loc_etot0-msim*vlight**2*(p%Gamma(index)-1)
ELSE
loc_etot0=loc_etot0-0.5*msim*vlight**2*(abs(p%URold(index)*p%UR(index))+abs(p%UZold(index)*p%UZ(index))+abs(p%UTHETold(index)*p%UTHET(index)))
END IF
! We fill the gap
CALL move_part(p, p%Nptot, index)
p%partindex(p%Nptot)=-1
! Reduce the total number of simulated parts
p%Nptot=p%Nptot-1
END IF
END SUBROUTINE delete_part
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Move particle with index sourceindex to particle with index destindex.
!> !WARNING! This will overwrite particle at destindex.
!
!> @param [in] sourceindex index in parts of the particle to move.
!> @param [in] destindex index in parts of the moved particle destination.
!---------------------------------------------------------------------------
SUBROUTINE move_part(p, sourceindex, destindex)
!! This will destroy particle at destindex
INTEGER, INTENT(IN) :: destindex, sourceindex
TYPE(particles), INTENT(INOUT)::p
! Move part at sourceindex in part at destindex
p%partindex(destindex)= p%partindex(sourceindex)
p%Gamma(destindex) = p%Gamma(sourceindex)
p%R(destindex) = p%R(sourceindex)
p%Z(destindex) = p%Z(sourceindex)
p%THET(destindex) = p%THET(sourceindex)
p%UR(destindex) = p%UR(sourceindex)
p%UTHET(destindex) = p%UTHET(sourceindex)
p%UZ(destindex) = p%UZ(sourceindex)
p%Rindex(destindex) = p%Rindex(sourceindex)
p%Zindex(destindex) = p%Zindex(sourceindex)
END SUBROUTINE move_part
!________________________________________________________________________________
SUBROUTINE destroy_parts(p)
TYPE(particles) :: p
p%nptot=0
IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z)
IF(ALLOCATED(p%R)) DEALLOCATE(p%R)
IF(ALLOCATED(p%THET)) DEALLOCATE(p%THET)
IF(ALLOCATED(p%WZ)) DEALLOCATE(p%WZ)
IF(ALLOCATED(p%WR)) DEALLOCATE(p%WR)
IF(ASSOCIATED(p%UR)) DEALLOCATE(p%UR)
IF(Associated(p%URold)) DEALLOCATE(p%URold)
IF(Associated(p%UZ)) DEALLOCATE(p%UZ)
IF(Associated(p%UZold)) DEALLOCATE(p%UZold)
IF(Associated(p%UTHET)) DEALLOCATE(p%UTHET)
IF(Associated(p%UTHETold)) DEALLOCATE(p%UTHETold)
IF(Associated(p%Gamma)) DEALLOCATE(p%Gamma)
IF(Associated(p%Gammaold)) DEALLOCATE(p%Gammaold)
IF(ALLOCATED(p%Rindex)) DEALLOCATE(p%Rindex)
IF(ALLOCATED(p%Zindex)) DEALLOCATE(p%Zindex)
IF(ALLOCATED(p%partindex)) DEALLOCATE(p%partindex)
END SUBROUTINE
!________________________________________________________________________________
SUBROUTINE clean_beam
!
CALL destroy_parts(parts)
!
END SUBROUTINE clean_beam
!________________________________________________________________________________
SUBROUTINE swappointer( pointer1, pointer2)
REAL(kind=db), DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2
REAL(kind=db), DIMENSION(:), POINTER:: temppointer
temppointer=>pointer1
pointer1=>pointer2
pointer2=>temppointer
END SUBROUTINE swappointer
!_______________________________________________________________________________
SUBROUTINE loaduniformRZ(VR,VZ,VTHET)
USE basic, ONLY: plasmadim, rnorm, temp, vnorm
USE constants, ONLY: me, kb, elchar
REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:)
REAL(kind=db):: vth
! Initial distribution in z with normalisation
CALL loduni(1,parts%Z(1:parts%Nptot))
parts%Z(1:parts%Nptot)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(1:parts%Nptot))/rnorm
! Initial distribution in r with normalisation
CALL loduni(2,parts%R(1:parts%Nptot))
parts%R(1:parts%Nptot)=(plasmadim(3)+parts%R(1:parts%Nptot)*(plasmadim(4)-plasmadim(3)))/rnorm
! Initial velocities distribution
vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity
CALL lodgaus(3,VZ)
CALL lodgaus(5,VR)
CALL lodgaus(7,VTHET)
VZ=VZ*vth
VR=VR*vth
VTHET=VTHET*vth
END SUBROUTINE loaduniformRZ
!_______________________________________________________________________________
SUBROUTINE loadDavidson(VR,VZ,VTHET,lodr)
USE constants, ONLY: me, kb, elchar
USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, B0, width, qsim, msim, vnorm, &
& omegac, zgrid, nz, rnorm, V, n0, nblock, temp
interface
subroutine lodr(nbase,y,ra,rb)
USE constants
REAL(kind=db), INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
REAL(kind=db), INTENT(in) :: rb, ra
end subroutine
end interface
REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:)
REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb
REAL(kind=db) :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL, vth
INTEGER :: i, j, n, blockstart, blockend, addedpart, remainparts
INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize
Allocate(ra(nblock),rb(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)
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)
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*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock)
qsim=VOL*n0*elchar/nplasma
msim=abs(qsim)/elchar*me
V=VOL/rnorm**3
!H0=me*omegac**2*r0**2*0.5
gamma=1/sqrt(1-omegac**2*r0**2/vlight**2)
!P0=H0/omegac
blockstart=1
blockend=0
ALLOCATE(blocksize(nblock))
WRITE(*,*) "blocksize: ", size(blocksize), nblock
DO n=1,nblock
blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock
END DO
remainparts=parts%Nptot-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
DO n=1,nblock
blockstart=blockend+1
blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot)
! Initial distribution in z with normalisation between magnetic mirrors
CALL loduni(1,parts%Z(blockstart:blockend))
parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm
CALL lodr(2,parts%R(blockstart:blockend),ra(n), rb(n))
parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm
END DO
IF(distribtype .eq. 5) THEN
! Initial velocities distribution
vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity
CALL lodgaus(3,VZ)
CALL lodgaus(5,VR)
CALL lodgaus(7,VTHET)
VZ=VZ*vth/4
VR=VR*8*vth
VTHET=VTHET*8*vth
ELSE
! Load velocities theta velocity
! Loading of r and z velocity is done in adapt_vinit to have
! access to parts%pot
DO i=1,parts%Nptot
! Interpolation for Magnetic potential
rg=parts%R(i)*rnorm
zg=(parts%Z(i)-halfLz)*rnorm
Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width))
Pcomp=P0/rg/me
Acomp=-SIGN(elchar/me*Athetpos,qsim)
VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp)
!VTHET(i)=Pcomp+Acomp
END DO
VTHET=VTHET/vnorm
VZ=0._db
VR=0._db
END IF
END SUBROUTINE loadDavidson
!===============================================================================
SUBROUTINE Temp_rescale
USE basic, ONLY: temprescale, nlclassical, nz
! INTEGER:: i, gridindex
! REAL(kind=db) :: vr, vz, vthet
! DO i=1,parts%Nptot
! gridindex=parts%zindex(i)+1+parts%rindex(i)*nz
! vr=parts%UR(i)/parts%Gamma(i)-fluidur(gridindex)
! vz=parts%UZ(i)/parts%Gamma(i)-fluiduz(gridindex)
! vthet=parts%UTHET(i)/parts%Gamma(i)-fluiduthet(gridindex)
! parts%UR(i)=vr*temprescale+fluidur(gridindex)
! parts%UTHET(i)=vthet*temprescale+fluiduthet(gridindex)
! parts%UZ(i)=vz*temprescale+fluiduz(gridindex)
! IF(nlclassical) THEN
! parts%Gamma(i)=1.0
! ELSE
! parts%Gamma(i)=1/sqrt(1-(parts%UR(i)**2+parts%UTHET(i)**2+parts%UZ(i)**2))
! END IF
! parts%UR(i)=parts%UR(i)*parts%Gamma(i)
! parts%UTHET(i)=parts%UTHET(i)*parts%Gamma(i)
! parts%UZ(i)=parts%UZ(i)*parts%Gamma(i)
! END DO
END SUBROUTINE Temp_rescale
!---------------------------------------------------------------------------
!> @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, qsim, msim, dt, dr, dz
INTEGER, INTENT(IN) ::samplefactor
TYPE(particles), INTENT(INOUT):: p
INTEGER:: i, j, currentindex
REAL(kind=db), DIMENSION(p%Nptot) :: 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(dz,minval(dr))/100
DO i=1,p%Nptot
DO j=1,samplefactor-1
currentindex=p%Nptot+(i-1)*(samplefactor-1)+j
CALL move_part(p,i,currentindex)
p%partindex(currentindex)=currentindex
dir = spreaddir(i)+2*pi*j/samplefactor
p%R(currentindex)=p%R(currentindex) + dl*cos(dir)
p%Z(currentindex)=p%Z(currentindex) + dl*sin(dir)
END DO
p%partindex(i)=i
p%R(i)=p%R(i) + dl*cos(spreaddir(i))
p%Z(i)=p%Z(i) + dl*sin(spreaddir(i))
END DO
nplasma=nplasma*samplefactor
qsim=qsim/samplefactor
msim=msim/samplefactor
p%Nptot=p%Nptot*samplefactor
END SUBROUTINE upsample
END 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.
!------------------------------------------------------------------------------
MODULE beam
!
USE constants
USE mpi
USE mpihelper
USE basic, ONLY: mpirank, mpisize
USE distrib
IMPLICIT NONE
!> Stores the particles properties for the run.
TYPE particles
INTEGER :: Nptot !< Local number of simulated particles
REAL(kind=db) :: m !< Macroparticle mass
REAL(kind=db) :: q !< Macroparticle charge
REAL(kind=db) :: qmRatio !< Charge over mass ratio
INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex !< Index in the electric potential grid for the R direction
INTEGER, DIMENSION(:), ALLOCATABLE :: Zindex !< Index in the electric potential grid for the Z direction
INTEGER, DIMENSION(:), ALLOCATABLE :: partindex !< Index of the particle to be able to follow it when it goes from one MPI host to the other
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: R !< radial coordinates of the particles
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Z !< longitudinal coordinates of the particles
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: THET !< azimuthal coordinates of the particles
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: WZ !< axial radial relative distances to the left grid line
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: WR !< radial relative distances to the bottom grid line
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: pot !< Electric potential
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Er !< Radial Electric field
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Ez !< Axial electric field
REAL(kind=db), DIMENSION(:),POINTER:: UR !< normalized radial velocity at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: URold !< normalized radial velocity at the previous time step
REAL(kind=db), DIMENSION(:),POINTER:: UTHET !< normalized azimuthal velocity at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: UTHETold !< normalized azimuthal velocity at the previous time step
REAL(kind=db), DIMENSION(:),POINTER:: UZ !< normalized axial velocity at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: UZold !< normalized axial velocity at the previous time step
REAL(kind=db), DIMENSION(:),POINTER:: Gamma !< Lorentz factor at the current time step
REAL(kind=db), DIMENSION(:),POINTER:: Gammaold !< Lorentz factor at the previous time step
END TYPE particles
!
TYPE(particles) :: parts !< Storage for all the particles
SAVE :: parts
TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffers to send and receive particle from left and right processes
! Diagnostics (scalars)
REAL(kind=db) :: ekin=0 !< Total kinetic energz (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)
INTEGER, DIMENSION(3), SAVE :: ireducerequest=MPI_REQUEST_NULL !< MPI requests used for the IREDUCE in the diagnostic routine
!
INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process
LOGICAL:: collected=.false. !< Stores if the particles data have been collected to MPI root process during this timestep
!
CONTAINS
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Allocate the memory for the particles variable storing the particles quantities.
!
!> @param[inout] p the particles variable needing to be allocated.
!> @param[in] nparts the maximum number of particles that will be stored in this variable
!---------------------------------------------------------------------------
SUBROUTINE creat_parts(p, nparts)
TYPE(particles) :: p
INTEGER, INTENT(in) :: nparts
IF (.NOT. ALLOCATED(p%Z) ) THEN
p%Nptot = nparts
ALLOCATE(p%Z(nparts))
ALLOCATE(p%R(nparts))
ALLOCATE(p%THET(nparts))
ALLOCATE(p%WZ(nparts))
ALLOCATE(p%WR(nparts))
ALLOCATE(p%UR(nparts))
ALLOCATE(p%UZ(nparts))
ALLOCATE(p%UTHET(nparts))
ALLOCATE(p%URold(nparts))
ALLOCATE(p%UZold(nparts))
ALLOCATE(p%UTHETold(nparts))
ALLOCATE(p%Gamma(nparts))
ALLOCATE(p%Rindex(nparts))
ALLOCATE(p%Zindex(nparts))
ALLOCATE(p%partindex(nparts))
ALLOCATE(p%pot(nparts))
ALLOCATE(p%Er(nparts))
ALLOCATE(p%Ez(nparts))
ALLOCATE(p%GAMMAold(nparts))
p%partindex=-1
p%URold=0
p%UZold=0
p%UTHETold=0
p%rindex=0
p%zindex=0
p%WR=0
p%WZ=0
p%UR=0
p%UZ=0
p%UTHET=0
p%Z=0
p%R=0
p%THET=0
p%Gamma=1
p%Er=0
p%Ez=0
p%pot=0
p%gammaold=1
END IF
END SUBROUTINE creat_parts
!---------------------------------------------------------------------------
!> @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, mpisize, nlclassical
USE mpi
INTEGER:: i
REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET
ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma))
parts%Nptot=nplasma
! Select case to define the type of distribution
SELECT CASE(distribtype)
CASE(1) ! Gaussian distribution in V and uniform in R
CALL loaduniformRZ(VR, VZ, VTHET)
CASE(2) !Stable distribution from Davidson 4.95 p.119
CALL loadDavidson(VR, VZ, VTHET, lodunir)
CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2
CALL loadDavidson(VR, VZ, VTHET, lodinvr)
CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R
CALL loadDavidson(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(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(VR, VZ, VTHET)
VZ = VZ/50
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
parts%partindex(i)=i
END DO
IF(nlclassical) THEN
parts%Gamma(1:nplasma)=1
ELSE
parts%Gamma(1:nplasma)=sqrt(1/(1-VR**2-VZ**2-VTHET**2))
END IF
! Normalization of the velocities
parts%UR(1:nplasma)=parts%Gamma(1:nplasma)*VR
parts%UZ(1:nplasma)=parts%Gamma(1:nplasma)*VZ
parts%UTHET(1:nplasma)=parts%Gamma(1:nplasma)*VTHET
DEALLOCATE(VZ, VR, VTHET)
CALL localisation(parts)
IF(mpisize .gt. 1) THEN
CALL distribpartsonprocs
END IF
END SUBROUTINE load_parts
!---------------------------------------------------------------------------
! 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.
!
!> @author Guillaume Le Bars EPFL/SPC
!---------------------------------------------------------------------------
SUBROUTINE bound(p)
USE basic, ONLY: zgrid, rnorm, nz, Zbounds, cstep, mpirank, partperiodic, step, leftproc, rightproc
type(particles), INTENT(INOUT):: p
INTEGER :: i, rsendnbparts, lsendnbparts, nblostparts
INTEGER, DIMENSION(p%Nptot) :: sendhole
INTEGER, DIMENSION(p%Nptot) :: losthole
LOGICAL:: leftcomm, rightcomm
rsendnbparts=0
lsendnbparts=0
nblostparts=0
losthole=0
sendhole=0
! We communicate with the left processus
leftcomm = leftproc .ne. -1 .or. partperiodic
! We communicate with the right processus
rightcomm = rightproc .ne. -1 .or. partperiodic
IF (p%Nptot .NE. 0) THEN
! Boundary condition at z direction
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,p%Nptot
IF(step .ne. 0) THEN
! If the particle is to the right of the local simulation space, it is sent to the right MPI process
IF (p%Z(i) .ge. zgrid(Zbounds(mpirank+1))) THEN
!$OMP CRITICAL (nbparts)
IF(rightcomm) THEN
rsendnbparts=rsendnbparts+1
sendhole(lsendnbparts+rsendnbparts)=i
ELSE
!WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in z"
!WRITE(*,*) "!!!!!!!!!! Particle removed !!!!!!!!!!"
nblostparts=nblostparts+1
losthole(nblostparts)=i
END IF
!$OMP END CRITICAL (nbparts)
! If the particle is to the left of the local simulation space, it is sent to the left MPI process
ELSE IF (p%Z(i) .lt. zgrid(Zbounds(mpirank))) THEN
!$OMP CRITICAL (nbparts)
IF(leftcomm) THEN
lsendnbparts=lsendnbparts+1
sendhole(lsendnbparts+rsendnbparts)=-i
ELSE
!WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in z"
!WRITE(*,*) "!!!!!!!!!! Particle removed !!!!!!!!!!"
nblostparts=nblostparts+1
losthole(nblostparts)=i
END IF
!$OMP END CRITICAL (nbparts)
END IF
END IF
! The periodic boundary conditions in z are here applied
DO WHILE (p%Z(i) .GT. zgrid(nz))
p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0)
END DO
DO WHILE (p%Z(i) .LT. zgrid(0))
p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0)
END DO
END DO
!$OMP END PARALLEL DO
END IF
IF(mpisize .gt. 1 .and. step .ne. 0) THEN
! We send the particles leaving the local simulation space to the closest neighbour
CALL particlescommunication(p, lsendnbparts, rsendnbparts, sendhole)
END IF
! If the boundary conditions are not periodic, we delete the corresponding particles
IF(nblostparts .gt. 0) THEN
DO i=nblostparts,1,-1
!WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in z: ", p%Z(losthole(i))*rnorm, zgrid(0)*rnorm, zgrid(nz)*rnorm
CALL delete_part(p, losthole(i))
END DO
WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank
END IF
END subroutine bound
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Compute the grid cell indices for each particle as well as the distance weight Wr, Wz.
!---------------------------------------------------------------------------
SUBROUTINE localisation(p)
USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm
type(particles), INTENT(INOUT):: p
INTEGER :: i, nblostparts
INTEGER, DIMENSION(p%Nptot) :: losthole
REAL(kind=db):: invdz, invdr(3)
losthole=0
nblostparts=0
invdz=1/dz
invdr=1/dr
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,p%Nptot
p%zindex(i)=(p%Z(i)-zgrid(0))*invdz
IF (p%R(i) .GT. rgrid(0) .AND. p%R(i) .LT. rgrid(nnr(1))) THEN
p%rindex(i)=(p%R(i)-rgrid(0))*invdr(1)
ELSE IF(p%R(i) .GE. rgrid(nnr(1)) .AND. p%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN
p%rindex(i)=(p%R(i)-rgrid(nnr(1)))*invdr(2)+nnr(1)
ELSE IF(p%R(i) .GE. rgrid(nnr(1)+nnr(2)) .AND. p%R(i) .LT. rgrid(nr)) THEN
p%rindex(i)=(p%R(i)-rgrid(nnr(1)+nnr(2)))*invdr(3)+nnr(1)+nnr(2)
ELSE
! If the particle is outside of the simulation space in the r direction, it is deleted.
!$OMP CRITICAL (lostparts)
nblostparts=nblostparts+1
losthole(nblostparts)=i
!$OMP END CRITICAL (lostparts)
END IF
END DO
!$OMP END PARALLEL DO
IF(nblostparts.gt.0) THEN
DO i=nblostparts,1,-1
!WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in r: ", p%R(losthole(i))*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm
CALL delete_part(p,losthole(i))
END DO
WRITE(*,'(i6.2,a,i6.2)') nblostparts, " particles lost in r on process: ", mpirank
END IF
!$OMP PARALLEL DO SIMD DEFAULT(SHARED)
DO i=1,p%Nptot
p%WZ(i)=(p%Z(i)-zgrid(p%zindex(i)))/dz;
p%WR(i)=(p%R(i)-rgrid(p%rindex(i)))/(rgrid(p%rindex(i)+1)-rgrid(p%rindex(i)));
END DO
!$OMP END PARALLEL DO SIMD
END IF
END SUBROUTINE localisation
!---------------------------------------------------------------------------
!> @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.
!
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity(p)
!
! Computes the new velocity of the particles due to Lorentz force
!
USE basic, ONLY : nlclassical
type(particles), INTENT(INOUT):: p
! Store old Velocities
CALL swappointer(p%UZold, p%UZ)
CALL swappointer(p%URold, p%UR)
CALL swappointer(p%UTHETold, p%UTHET)
CALL swappointer(p%Gammaold, p%Gamma)
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$
!---------------------------------------------------------------------------
SUBROUTINE comp_velocity_fun(p, gamma)
!
! Computes the new velocity of the particles due to Lorentz force
!
USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz
interface
subroutine gamma(gam, UZ, UR, UTHET)
USE constants
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
REAL(kind=db), INTENT(OUT):: gam
end subroutine
end interface
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: tau
REAL(kind=db):: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2
INTEGER:: J1, J2, J3, J4
INTEGER:: i
! Normalized time increment
tau=omegac/2/omegap*dt/tnorm
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
DO i=1,p%Nptot
J1=(p%rindex(i))*(nz+1) + p%zindex(i)+1
J2=(p%rindex(i))*(nz+1) + p%zindex(i)+2
J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1
J4=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+2
! Interpolation for magnetic field
BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) &
& +p%WZ(i)*(1-p%WR(i))*Bz(J3) &
& +(1-p%WZ(i))*p%WR(i)*Bz(J2) &
& +p%WZ(i)*p%WR(i)*Bz(J1)
BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) &
& +p%WZ(i)*(1-p%WR(i))*Br(J3) &
& +(1-p%WZ(i))*p%WR(i)*Br(J2) &
& +p%WZ(i)*p%WR(i)*Br(J1)
! First half of electric pulse
p%UZ(i)=p%UZold(i)+p%Ez(i)*tau
p%UR(i)=p%URold(i)+p%ER(i)*tau
CALL gamma(p%Gamma(i), p%UZ(i), p%UR(i), p%UTHETold(i))
! Rotation along magnetic field
ZBZ=tau*BRZ/p%Gamma(i)
ZBR=tau*BRR/p%Gamma(i)
ZPZ=p%UZ(i)-ZBR*p%UTHETold(i) !u'_{z}
ZPR=p%UR(i)+ZBZ*p%UTHETold(i) !u'_{r}
ZPTHET=p%UTHETold(i)+(ZBR*p%UZ(i)-ZBZ*p%UR(i)) !u'_{theta}
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZBZ2=2*ZBZ/SQR
ZBR2=2*ZBR/SQR
p%UZ(i)=p%UZ(i)-ZBR2*ZPTHET !u+_{z}
p%UR(i)=p%UR(i)+ZBZ2*ZPTHET !u+_{r}
p%UTHET(i)=p%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta}
! Second half of acceleration
p%UZ(i)=p%UZ(i)+p%EZ(i)*tau
p%UR(i)=p%UR(i)+p%ER(i)*tau
! Final computation of the Lorentz factor
CALL gamma(p%Gamma(i), p%UZ(i), p%UR(i), p%UTHETold(i))
END DO
!$OMP END PARALLEL DO SIMD
END IF
collected=.false.
END SUBROUTINE comp_velocity_fun
!---------------------------------------------------------------------------
!> @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
!---------------------------------------------------------------------------
ELEMENTAL SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET)
!$OMP declare simd(gamma_classical)
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
REAL(kind=db), INTENT(OUT):: gamma
gamma=1.0
END SUBROUTINE 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
!---------------------------------------------------------------------------
ELEMENTAL SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET)
!$OMP declare simd(gamma_relativistic)
REAL(kind=db), INTENT(IN):: UR,UZ,UTHET
REAL(kind=db), INTENT(OUT):: gamma
gamma=sqrt(1+UZ**2+UR**2+UTHET**2)
END SUBROUTINE
!---------------------------------------------------------------------------
!> @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.
!---------------------------------------------------------------------------
SUBROUTINE push(p)
Use basic, ONLY: dt, tnorm
type(particles), INTENT(INOUT):: p
REAL(kind=db):: XP, YP, COSA, SINA, U1, U2, ALPHA
INTEGER :: i
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(XP, YP, COSA, SINA, U1, U2, ALPHA)
DO i=1,p%Nptot
! Local Cartesian coordinates
XP=p%R(i)+dt/tnorm*p%UR(i)/p%Gamma(i)
YP=dt/tnorm*p%UTHET(i)/p%Gamma(i)
! Conversion to cylindrical coordiantes
p%Z(i)=p%Z(i)+dt/tnorm*p%UZ(i)/p%Gamma(i)
p%R(i)=sqrt(XP**2+YP**2)
! Computation of the rotation angle
IF (p%R(i) .EQ. 0) THEN
COSA=1
SINA=0
ALPHA=0
ELSE
COSA=XP/p%R(i)
SINA=YP/p%R(i)
ALPHA=asin(SINA)
END IF
! New azimuthal position
p%THET(i)=MOD(p%THET(i)+ALPHA,2*pi)
! Velocity in rotated reference frame
U1=COSA*p%UR(i)+SINA*p%UTHET(i)
U2=-SINA*p%UR(i)+COSA*p%UTHET(i)
p%UR(i)=U1
p%UTHET(i)=U2
END DO
!$OMP END PARALLEL DO SIMD
END IF
collected=.false.
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 also computes the first and second order
!> moments of the distribution function
!
!---------------------------------------------------------------------------
SUBROUTINE diagnostics
!
! Compute energies
!
USE constants, ONLY: vlight
USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, ierr, step, it0d, it2d, nlend, nr, itparts, nplasma, partperiodic, newres, Zbounds
INTEGER :: i, gridindex
INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE)
CALL MPI_WAITALL(3, ireducerequest, stat, ierr)
! Reset the quantities
ekin=0
epot=0
etot=0
! Computation of the kinetic and potential energy as well as fluid velocities and density
!$OMP PARALLEL DO SIMD REDUCTION(+:epot, ekin) DEFAULT(SHARED)
DO i=1,parts%Nptot
! Potential energy
epot=epot+0.5*qsim*parts%pot(i)*phinorm
! Kinetic energy
IF(.not. nlclassical) THEN
ekin=ekin+msim*vlight**2*(parts%Gamma(i)-1)
ELSE
ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))
END IF
END DO
!$OMP END PARALLEL DO SIMD
! Shift to Etot at cstep=1 (not valable yet at cstep=0!)
IF(cstep.EQ.1 .or. (newres .AND. step .EQ. 1)) THEN
! Compute the local total energy
loc_etot0 = epot+ekin
END IF
etot=loc_etot0
! The computed energy is sent to the root process
IF(mpisize .gt.1) THEN
IF(mpirank .eq.0 ) THEN
CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
CALL MPI_IREDUCE(etot, etot0, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
ELSE
CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(1), ierr)
CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(2), ierr)
CALL MPI_IREDUCE(etot, etot0, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
& 0, MPI_COMM_WORLD, ireducerequest(3), ierr)
END IF
ELSE
etot0=loc_etot0
END IF
! Wait for the communication of the energy to be finished
CALL MPI_WAITALL(3, ireducerequest(1:3), stat(1:3*MPI_STATUS_SIZE), ierr)
!CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr)
! Compute the total energy
etot=epot+ekin
! Send the local number of particles on each node to the root process
IF(modulo(step,it0d).eq. 0 .and. mpisize .gt. 1) THEN
Nplocs_all(mpirank)=parts%Nptot
IF(mpirank .eq.0 ) THEN
CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
ELSE
CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,&
& 0, MPI_COMM_WORLD, ierr)
END IF
IF(mpirank .eq. 0 .AND. partperiodic .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN
!WRITE(*,*) "Error particle lost on some processus"
!CALL MPI_abort(MPI_COMM_WORLD,-1,ierr)
END IF
END IF
!Calls the communications to send the particles data to the root process for diagnostic file every itparts time steps
IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN
CALL collectparts(parts)
END IF
end subroutine diagnostics
!---------------------------------------------------------------------------
!> @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)
USE basic, ONLY: mpirank, mpisize, ierr
type(particles), INTENT(INOUT):: p
INTEGER, DIMENSION(:), ALLOCATABLE :: displs
INTEGER:: i
IF(collected) RETURN ! exit subroutine if particles have already been collected during this time step
IF(mpirank .ne. 0) THEN
ALLOCATE(displs(0:mpisize-1))
displs=0
! Send Particles informations to root process
CALL MPI_Gatherv(p%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%Z, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%R, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%THET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%THET, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UR, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UZ, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UTHET, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%pot, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%Rindex, Nplocs_all(mpirank), MPI_INTEGER, p%Rindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%Zindex, Nplocs_all(mpirank), MPI_INTEGER, p%Zindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(p%partindex, Nplocs_all(mpirank), MPI_INTEGER, p%partindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
ELSE
ALLOCATE(displs(0:mpisize-1))
displs(0)=0
Do i=1,mpisize-1
displs(i)=displs(i-1)+Nplocs_all(i-1)
END DO
! Receive particle information from all processes
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%Z, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%R, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%THET, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UR, Nplocs_all, displs,&
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UZ, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UTHET, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%pot, Nplocs_all, displs, &
& MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%Rindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%Zindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%partindex, Nplocs_all, displs, &
& MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
p%partindex(sum(Nplocs_all)+1:)=-1
END IF
collected=.TRUE.
END SUBROUTINE collectparts
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief Computes the velocities at time t-1/2 to keep the second order precision in time on the velocity.
!
!---------------------------------------------------------------------------
SUBROUTINE adapt_vinit(p)
!! Computes the velocity at time -dt/2 from velocities computed at time 0
!
USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm
type(particles), INTENT(INOUT):: p
REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, &
& SQR, ZBZ2, ZBR2, Vperp, v2
INTEGER :: J1, J2, J3, J4, 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) THEN
ALLOCATE(VR(p%Nptot),VZ(p%Nptot),VTHET(p%Nptot))
CALL loduni(7,VZ)
VZ=VZ*2*pi
VTHET=p%UTHET/p%Gamma*vnorm
DO i=1,p%Nptot
Vperp=sqrt(MAX(2*H0/me+2*elchar/me*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%UR(i)=p%Gamma(i)*VR(i)/vnorm
p%UZ(i)=p%Gamma(i)*VZ(i)/vnorm
p%UTHET(i)=p%Gamma(i)*VTHET(i)/vnorm
END DO
DEALLOCATE(VR,VZ,VTHET)
END IF
! Normalised time increment
tau=-omegac/2/omegap*dt/tnorm
! Store old Velocities
CALL swappointer(p%UZold, p%UZ)
CALL swappointer(p%URold, p%UR)
CALL swappointer(p%UTHETold, p%UTHET)
CALL swappointer(p%Gammaold, p%Gamma)
IF (p%Nptot .NE. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2)
DO i=1,p%Nptot
! Compute the particle linear indices for the magnetic field interpolation
J1=(p%rindex(i))*(nz+1)+p%zindex(i)+1
J2=J1+1
J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1
J4=J3+1
! Interpolation for magnetic field
BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) &
& +p%WZ(i)*(1-p%WR(i))*Bz(J3) &
& +(1-p%WZ(i))*p%WR(i)*Bz(J2) &
& +p%WZ(i)*p%WR(i)*Bz(J1)
BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) &
& +p%WZ(i)*(1-p%WR(i))*Br(J3) &
& +(1-p%WZ(i))*p%WR(i)*Br(J2) &
& +p%WZ(i)*p%WR(i)*Br(J1)
! Half inverse Rotation along magnetic field
ZBZ=tau*BRZ/p%Gammaold(i)
ZBR=tau*BRR/p%Gammaold(i)
SQR=1+ZBZ*ZBZ+ZBR*ZBR
ZPZ=(p%UZold(i)-ZBR*p%UTHETold(i))/SQR !u-_{z}
ZPR=(p%URold(i)+ZBZ*p%UTHETold(i))/SQR !u-_{r}
ZPTHET=p%UTHETold(i)+(ZBR*p%UZold(i)-ZBZ*p%URold(i))/SQR !u-_{theta}
p%UZ(i)=ZPZ
p%UR(i)=ZPR
p%UTHET(i)=ZPTHET
! half of decceleration
p%UZ(i)=p%UZ(i)+p%Ez(i)*tau
p%UR(i)=p%UR(i)+p%Er(i)*tau
IF(.not. nlclassical) THEN
p%Gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2)
END IF
END DO
!$OMP END PARALLEL DO
END IF
END SUBROUTINE adapt_vinit
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief In the case of MPI parallelism, computes the indices of the particle assigned to the current process.
!---------------------------------------------------------------------------
SUBROUTINE distribpartsonprocs
! 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, ierr, Zbounds
INTEGER:: k, i, nbparts
REAL(kind=db):: idealnbpartsperproc
INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k)
INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector
INTEGER:: partindexlast ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector
INTEGER, DIMENSION(0:nz):: 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, zindex
partspercol=0
idealnbpartsperproc = FLOOR(REAL(parts%Nptot)/REAL(mpisize))
Zmin=MINVAL(parts%Zindex(1:parts%Nptot))
Zmax=MAXVAL(parts%Zindex(1:parts%Nptot))
Zperproc=(Zmax-Zmin)/mpisize
IF(Zmax .eq. 0) Zmax=nz
DO k=1,parts%Nptot
zindex=parts%Zindex(k)
partspercol(zindex)=partspercol(zindex)+1
END DO
IF (Zperproc .eq. 0) THEN
!! No particles are present initially
Zperproc=nz/mpisize
Zmin=0
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
DO k=1,mpisize-1
nbparts=0
DO WHILE(nbparts<0.99*idealnbpartsperproc .and. i .lt. Zmax)
nbparts=nbparts+partspercol(i)
i=i+1
END DO
Zbounds(k)=i
END DO
END IF
partindexstart=1
totparts=0
partindexlast=parts%Nptot
!DO k=0,mpisize-2
! Nplocs_all(k)=idealnbpartsperproc
!END DO
! NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize)
!DO k=1, mpisize-1
! Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc))
!END DO
DO k=0,mpisize-1
Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1))
END DO
! Store local end of particle vector
partindexlast=SUM(Nplocs_all(0:mpirank))
! Store local start of particle vector
IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1
IF(mpirank .ne. 0 .and. Nplocs_all(mpirank) .gt. 0) THEN
DO k= partindexstart, partindexlast
CALL move_part(parts,k,k-partindexstart+1)
END DO
END IF
IF(INT(SUM(Nplocs_all)) .ne. parts%Nptot) THEN
WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " Nptot:", parts%Nptot
CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr)
END IF
parts%Nptot=Nplocs_all(mpirank)
WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot
END SUBROUTINE distribpartsonprocs
!_______________________________________________________________________________
!---------------------------------------------------------------------------
!> @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)
!> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative
!---------------------------------------------------------------------------
SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, sendholes)
USE mpihelper, ONLY: particle_type
USE basic, ONLY: rightproc, leftproc, ierr
type(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(in) :: sendholes(:)
INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0
INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest
INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE)
INTEGER :: lsentnbparts, rsentnbparts
INTEGER :: lreceivednbparts, rreceivednbparts
lsentnbparts=lsendnbparts
rsentnbparts=rsendnbparts
sendrequest=MPI_REQUEST_NULL
recvrequest=MPI_REQUEST_NULL
! Send and receive the number of particles to exchange
CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, recvrequest(1), ierr)
CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, recvrequest(2), ierr)
CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, sendrequest(1), ierr)
CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, sendrequest(2), ierr)
CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr)
CALL MPI_Wait(recvrequest(2), recvstatus(2), ierr)
recvrequest=MPI_REQUEST_NULL
lreceivednbparts=lrecvnbparts
rreceivednbparts=rrecvnbparts
! Re/allocate enough memory to store the incoming particles
IF( ALLOCATED(rrecvpartbuff) .AND. (size(rrecvpartbuff) .lt. rrecvnbparts) ) DEALLOCATE(rrecvpartbuff)
IF(.not. ALLOCATED(rrecvpartbuff) .and. rrecvnbparts .gt. 0) ALLOCATE(rrecvpartbuff(INT(rreceivednbparts*1.3)))
IF( ALLOCATED(lrecvpartbuff) .AND. (size(lrecvpartbuff) .lt. lrecvnbparts) ) DEALLOCATE(lrecvpartbuff)
IF((.not. ALLOCATED(lrecvpartbuff) ).and. lrecvnbparts .gt. 0) ALLOCATE(lrecvpartbuff(INT(lreceivednbparts*1.3)))
! Receive particles from left and right processes to the corresponding buffers
IF ( lrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, recvrequest(1), ierr)
END IF
IF( rrecvnbparts .gt. 0) THEN
CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, recvrequest(2), ierr)
END IF
! Copy the leaving particles to the corresponding send buffers
IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN
CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes)
END IF
CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr)
CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr)
! Send the particles to the left and right neighbours
IF( lsendnbparts .gt. 0) THEN
CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, sendrequest(1), ierr)
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, sendrequest(2), ierr)
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
END IF
IF ( rreceivednbparts .gt. 0) THEN
CALL MPI_Wait(recvrequest(2), recvstatus(MPI_STATUS_SIZE+1), 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
END IF
! Copy the incoming particles from the receive buffers to the simulation parts variable
CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes)
! 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)
END IF
IF( rsendnbparts .gt. 0) THEN
CALL MPI_Wait(sendrequest(2), sendstatus(2), ierr)
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
!> @param [in] sendholes array containing the indices of the particle having left the local domain in ascending order.
!---------------------------------------------------------------------------
SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts, sendholes)
!
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts
INTEGER, INTENT(in) :: sendholes(:)
INTEGER k,partpos, partdiff
! computes if we received less particles than we sent
partdiff=sendnbparts-rrecvnbparts-lrecvnbparts
! 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
partpos=abs(sendholes(k))
ELSE
! Add at the end of parts
p%Nptot=p%Nptot+1
partpos=p%Nptot
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
partpos=abs(sendholes(k+rrecvnbparts))
ELSE
p%Nptot=p%Nptot+1
partpos=p%Nptot
END IF
CALL Insertincomingpart(p, lrecvpartbuff, k, partpos)
END DO
END IF
! If we received less particles than we sent, we fill the remaining holes with the particles from the end of the
! parts arrays
IF(partdiff .gt. 0) THEN
DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts
CALL move_part(parts, parts%Nptot, abs(sendholes(k)))
p%Nptot = p%Nptot-1
END DO
END IF
!
END SUBROUTINE Addincomingparts
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Copy one particle from the receive buffers to the local simulation variable parts.
!
!> @param [in] buffer receive buffer containing the particles parameters to copy from
!> @param [in] bufferindex particle index in the receive buffer
!> @param [in] partsindex destination particle index in the local parts variable
!---------------------------------------------------------------------------
SUBROUTINE Insertincomingpart(p, buffer, bufferindex, partsindex)
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: bufferindex, partsindex
TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer
p%partindex(partsindex) = buffer(bufferindex)%partindex
p%Rindex(partsindex) = buffer(bufferindex)%Rindex
p%Zindex(partsindex) = buffer(bufferindex)%Zindex
p%R(partsindex) = buffer(bufferindex)%R
p%Z(partsindex) = buffer(bufferindex)%Z
p%THET(partsindex) = buffer(bufferindex)%THET
p%UZ(partsindex) = buffer(bufferindex)%UZ
p%UR(partsindex) = buffer(bufferindex)%UR
p%UTHET(partsindex) = buffer(bufferindex)%UTHET
p%Gamma(partsindex) = buffer(bufferindex)%Gamma
!
!
END SUBROUTINE Insertincomingpart
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Copy one particle from the local parts variable to the send buffer.
!
!> @param [in] buffer send buffer to copy to
!> @param [in] bufferindex particle index in the send buffer
!> @param [in] partsindex origin particle index in the local parts variable
!---------------------------------------------------------------------------
SUBROUTINE Insertsentpart(p, buffer, bufferindex, partsindex)
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: bufferindex, partsindex
TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer
buffer(bufferindex)%partindex = p%partindex(partsindex)
buffer(bufferindex)%Rindex = p%Rindex(partsindex)
buffer(bufferindex)%Zindex = p%Zindex(partsindex)
buffer(bufferindex)%R = p%R(partsindex)
buffer(bufferindex)%Z = p%Z(partsindex)
buffer(bufferindex)%THET = p%THET(partsindex)
buffer(bufferindex)%UZ = p%UZ(partsindex)
buffer(bufferindex)%UR = p%UR(partsindex)
buffer(bufferindex)%UTHET = p%UTHET(partsindex)
buffer(bufferindex)%Gamma = p%Gamma(partsindex)
!
END SUBROUTINE Insertsentpart
!---------------------------------------------------------------------------
!> @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)
!> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative
!---------------------------------------------------------------------------
SUBROUTINE AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes)
!
USE mpihelper
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts
INTEGER, INTENT(in) :: sendholes(:)
INTEGER:: partpos, k
INTEGER:: lsendpos, rsendpos
lsendpos=0
rsendpos=0
! Verify if the buffers are big enough to store the outgoing particles
IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff)
IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff)
! Allocate the buffers if needed
IF(.not. ALLOCATED(rsendpartbuff) .and. rsendnbparts .gt. 0) ALLOCATE(rsendpartbuff(INT(rsendnbparts*1.3)))
IF(.not. ALLOCATED(lsendpartbuff) .and. lsendnbparts .gt. 0) ALLOCATE(lsendpartbuff(INT(lsendnbparts*1.3)))
! Loop over the outgoing particles and fill the correct send buffer
Do k=lsendnbparts+rsendnbparts,1,-1
partpos=abs(sendholes(k))
IF(sendholes(k) .GT. 0) THEN
rsendpos=rsendpos+1
CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos)
ELSE IF(sendholes(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 Exchange two particles in the parts variable.
!
!> @param [in] index1 index in parts of the first particle to exchange.
!> @param [in] index2 index in parts of the second particle to exchange.
!---------------------------------------------------------------------------
SUBROUTINE exchange_parts(p, index1, index2)
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(IN) :: index1, index2
REAL(kind=db):: R, Z, THET, UR, UZ, UTHET, Gamma
INTEGER :: Rindex, Zindex, partindex
!! Exchange particle at index1 with particle at index2
! Store part at index1 in temporary value
partindex= p%partindex(index1)
Gamma = p%Gamma(index1)
R = p%R(index1)
Z = p%Z(index1)
THET = p%THET(index1)
UR = p%UR(index1)
UTHET = p%UTHET(index1)
UZ = p%UZ(index1)
Rindex = p%Rindex(index1)
Zindex = p%Zindex(index1)
! Move part at index2 in part at index 1
p%partindex(index1)= p%partindex(index2)
p%Gamma(index1) = p%Gamma(index2)
p%R(index1) = p%R(index2)
p%Z(index1) = p%Z(index2)
p%THET(index1) = p%THET(index2)
p%UR(index1) = p%UR(index2)
p%UTHET(index1) = p%UTHET(index2)
p%UZ(index1) = p%UZ(index2)
p%Rindex(index1) = p%Rindex(index2)
p%Zindex(index1) = p%Zindex(index2)
! Move temporary values from part(index1) to part(index2)
p%partindex(index2)= partindex
p%Gamma(index2) = Gamma
p%R(index2) = R
p%Z(index2) = Z
p%THET(index2) = THET
p%UR(index2) = UR
p%UTHET(index2) = UTHET
p%UZ(index2) = UZ
p%Rindex(index2) = Rindex
p%Zindex(index2) = Zindex
END SUBROUTINE exchange_parts
SUBROUTINE change_parts_allocation(p, sizedifference)
implicit none
TYPE(particles), INTENT(INOUT):: p
INTEGER,INTENT(IN) :: sizedifference
CALL change_array_size_int(p%Rindex, sizedifference)
CALL change_array_size_int(p%Zindex, sizedifference)
CALL change_array_size_int(p%partindex, sizedifference)
CALL change_array_size_dp(p%ER,sizedifference)
CALL change_array_size_dp(p%EZ,sizedifference)
CALL change_array_size_dp(p%pot,sizedifference)
CALL change_array_size_dp(p%R,sizedifference)
CALL change_array_size_dp(p%Z,sizedifference)
CALL change_array_size_dp(p%THET,sizedifference)
CALL change_array_size_dp(p%WR,sizedifference)
CALL change_array_size_dp(p%WZ,sizedifference)
CALL change_array_size_dp_ptr(p%UR,sizedifference)
CALL change_array_size_dp_ptr(p%URold,sizedifference)
CALL change_array_size_dp_ptr(p%UZ,sizedifference)
CALL change_array_size_dp_ptr(p%UZold,sizedifference)
CALL change_array_size_dp_ptr(p%UTHET,sizedifference)
CALL change_array_size_dp_ptr(p%UTHETold,sizedifference)
CALL change_array_size_dp_ptr(p%Gamma,sizedifference)
CALL change_array_size_dp_ptr(p%Gammaold,sizedifference)
p%Nptot=MIN(p%Nptot,size(p%R)+sizedifference)
END SUBROUTINE change_parts_allocation
SUBROUTINE change_array_size_dp(arr, sizedifference)
implicit none
REAL(kind=db), ALLOCATABLE, INTENT(INOUT):: arr(:)
INTEGER, INTENT(IN):: sizedifference
REAL(kind=db), ALLOCATABLE:: temp(:)
INTEGER:: current_size, new_size
if(allocated(arr)) THEN
current_size=size(arr)
new_size=current_size+sizedifference
ALLOCATE(temp(new_size))
temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size))
DEALLOCATE(arr)
CALL move_alloc(temp, arr)
END IF
END SUBROUTINE change_array_size_dp
SUBROUTINE change_array_size_dp_ptr(arr, sizedifference)
implicit none
REAL(kind=db), POINTER, INTENT(INOUT):: arr(:)
INTEGER, INTENT(IN):: sizedifference
REAL(kind=db), POINTER:: temp(:)
INTEGER:: current_size, new_size
if(associated(arr)) THEN
current_size=size(arr)
new_size=current_size+sizedifference
ALLOCATE(temp(new_size))
temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size))
DEALLOCATE(arr)
arr=> temp
END IF
END SUBROUTINE change_array_size_dp_ptr
SUBROUTINE change_array_size_int(arr, sizedifference)
implicit none
INTEGER, ALLOCATABLE, INTENT(INOUT):: arr(:)
INTEGER, INTENT(IN):: sizedifference
INTEGER, ALLOCATABLE:: temp(:)
INTEGER:: current_size, new_size
if(allocated(arr)) THEN
current_size=size(arr)
new_size=current_size+sizedifference
ALLOCATE(temp(new_size))
temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size))
DEALLOCATE(arr)
CALL move_alloc(temp,arr)
END IF
END SUBROUTINE change_array_size_int
SUBROUTINE add_created_part(p, buffer)
USE bsplines
USE basic, ONLY: splrz, qsim, phinorm, nlclassical, msim
IMPLICIT NONE
TYPE(particles), INTENT(INOUT):: p
TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:)
INTEGER:: i, nptotinit, parts_size_increase, nb_added
nptotinit=p%Nptot+1
nb_added=size(buffer,1)
! if there is not enough space in the parts simulation buffer, increase the parst size
IF(p%Nptot + nb_added .gt. size(p%Z,1)) THEN
parts_size_increase=Max(floor(0.1*size(p%Z,1)),nb_added)
CALL change_parts_allocation(p, parts_size_increase)
END IF
DO i=1,nb_added
p%Nptot=p%Nptot+1
CALL Insertincomingpart(p, buffer,i,p%Nptot)
END DO
! Compute the energy added by these new particles for energy conservation diagnostic
IF(nb_added .gt. 0) THEN
CALL getgrad(splrz, p%Z(nptotinit:p%Nptot), p%R(nptotinit:p%Nptot), &
& p%pot(nptotinit:p%Nptot), p%EZ(nptotinit:p%Nptot), p%ER(nptotinit:p%Nptot))
p%EZ(nptotinit:p%Nptot)=-p%Ez(nptotinit:p%Nptot)
p%ER(nptotinit:p%Nptot)=-p%ER(nptotinit:p%Nptot)
loc_etot0=loc_etot0+qsim*sum(p%pot(nptotinit:p%Nptot))*phinorm
IF(.not. nlclassical) THEN
loc_etot0=loc_etot0+sum(msim*vlight**2*(p%Gamma(nptotinit:p%Nptot)-1))
ELSE
loc_etot0=loc_etot0+0.5*msim*vlight**2*sum(p%UR(nptotinit:p%Nptot)**2 &
& +p%UZ(nptotinit:p%Nptot)**2 &
& +p%UTHET(nptotinit:p%Nptot)**2)
END IF
END IF
END SUBROUTINE
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Delete particle at given index removing its energy from the system
!
!> @param [in] index index of particle to be deleted
!---------------------------------------------------------------------------
SUBROUTINE delete_part(p, index)
!! This will destroy particle at index
USE constants, ONLY: vlight
USE basic, ONLY: qsim, phinorm, msim, nlclassical
TYPE(particles), INTENT(INOUT):: p
INTEGER, INTENT(IN) :: index
IF(index .le. p%Nptot) THEN
loc_etot0=loc_etot0-qsim*p%pot(index)*phinorm
IF(.not. nlclassical) THEN
loc_etot0=loc_etot0-msim*vlight**2*(p%Gamma(index)-1)
ELSE
loc_etot0=loc_etot0-0.5*msim*vlight**2*(abs(p%URold(index)*p%UR(index))+abs(p%UZold(index)*p%UZ(index))+abs(p%UTHETold(index)*p%UTHET(index)))
END IF
! We fill the gap
CALL move_part(p, p%Nptot, index)
p%partindex(p%Nptot)=-1
! Reduce the total number of simulated parts
p%Nptot=p%Nptot-1
END IF
END SUBROUTINE delete_part
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Move particle with index sourceindex to particle with index destindex.
!> !WARNING! This will overwrite particle at destindex.
!
!> @param [in] sourceindex index in parts of the particle to move.
!> @param [in] destindex index in parts of the moved particle destination.
!---------------------------------------------------------------------------
SUBROUTINE move_part(p, sourceindex, destindex)
!! This will destroy particle at destindex
INTEGER, INTENT(IN) :: destindex, sourceindex
TYPE(particles), INTENT(INOUT)::p
! Move part at sourceindex in part at destindex
p%partindex(destindex)= p%partindex(sourceindex)
p%Gamma(destindex) = p%Gamma(sourceindex)
p%R(destindex) = p%R(sourceindex)
p%Z(destindex) = p%Z(sourceindex)
p%THET(destindex) = p%THET(sourceindex)
p%UR(destindex) = p%UR(sourceindex)
p%UTHET(destindex) = p%UTHET(sourceindex)
p%UZ(destindex) = p%UZ(sourceindex)
p%Rindex(destindex) = p%Rindex(sourceindex)
p%Zindex(destindex) = p%Zindex(sourceindex)
END SUBROUTINE move_part
!________________________________________________________________________________
SUBROUTINE destroy_parts(p)
TYPE(particles) :: p
p%nptot=0
IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z)
IF(ALLOCATED(p%R)) DEALLOCATE(p%R)
IF(ALLOCATED(p%THET)) DEALLOCATE(p%THET)
IF(ALLOCATED(p%WZ)) DEALLOCATE(p%WZ)
IF(ALLOCATED(p%WR)) DEALLOCATE(p%WR)
IF(ASSOCIATED(p%UR)) DEALLOCATE(p%UR)
IF(Associated(p%URold)) DEALLOCATE(p%URold)
IF(Associated(p%UZ)) DEALLOCATE(p%UZ)
IF(Associated(p%UZold)) DEALLOCATE(p%UZold)
IF(Associated(p%UTHET)) DEALLOCATE(p%UTHET)
IF(Associated(p%UTHETold)) DEALLOCATE(p%UTHETold)
IF(Associated(p%Gamma)) DEALLOCATE(p%Gamma)
IF(Associated(p%Gammaold)) DEALLOCATE(p%Gammaold)
IF(ALLOCATED(p%Rindex)) DEALLOCATE(p%Rindex)
IF(ALLOCATED(p%Zindex)) DEALLOCATE(p%Zindex)
IF(ALLOCATED(p%partindex)) DEALLOCATE(p%partindex)
END SUBROUTINE
!________________________________________________________________________________
SUBROUTINE clean_beam
!
CALL destroy_parts(parts)
!
END SUBROUTINE clean_beam
!________________________________________________________________________________
SUBROUTINE swappointer( pointer1, pointer2)
REAL(kind=db), DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2
REAL(kind=db), DIMENSION(:), POINTER:: temppointer
temppointer=>pointer1
pointer1=>pointer2
pointer2=>temppointer
END SUBROUTINE swappointer
!_______________________________________________________________________________
SUBROUTINE loaduniformRZ(VR,VZ,VTHET)
USE basic, ONLY: plasmadim, rnorm, temp, vnorm
USE constants, ONLY: me, kb, elchar
REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:)
REAL(kind=db):: vth
! Initial distribution in z with normalisation
CALL loduni(1,parts%Z(1:parts%Nptot))
parts%Z(1:parts%Nptot)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(1:parts%Nptot))/rnorm
! Initial distribution in r with normalisation
CALL loduni(2,parts%R(1:parts%Nptot))
parts%R(1:parts%Nptot)=(plasmadim(3)+parts%R(1:parts%Nptot)*(plasmadim(4)-plasmadim(3)))/rnorm
! Initial velocities distribution
vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity
CALL lodgaus(3,VZ)
CALL lodgaus(5,VR)
CALL lodgaus(7,VTHET)
VZ=VZ*vth
VR=VR*vth
VTHET=VTHET*vth
END SUBROUTINE loaduniformRZ
!_______________________________________________________________________________
SUBROUTINE loadDavidson(VR,VZ,VTHET,lodr)
USE constants, ONLY: me, kb, elchar
USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, B0, width, qsim, msim, vnorm, &
& omegac, zgrid, nz, rnorm, V, n0, nblock, temp
interface
subroutine lodr(nbase,y,ra,rb)
USE constants
REAL(kind=db), INTENT(out) :: y(:)
INTEGER, INTENT(in) :: nbase
REAL(kind=db), INTENT(in) :: rb, ra
end subroutine
end interface
REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:)
REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb
REAL(kind=db) :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL, vth
INTEGER :: i, j, n, blockstart, blockend, addedpart, remainparts
INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize
Allocate(ra(nblock),rb(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)
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)
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*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock)
qsim=VOL*n0*elchar/nplasma
msim=abs(qsim)/elchar*me
V=VOL/rnorm**3
!H0=me*omegac**2*r0**2*0.5
gamma=1/sqrt(1-omegac**2*r0**2/vlight**2)
!P0=H0/omegac
blockstart=1
blockend=0
ALLOCATE(blocksize(nblock))
WRITE(*,*) "blocksize: ", size(blocksize), nblock
DO n=1,nblock
blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock
END DO
remainparts=parts%Nptot-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
DO n=1,nblock
blockstart=blockend+1
blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot)
! Initial distribution in z with normalisation between magnetic mirrors
CALL loduni(1,parts%Z(blockstart:blockend))
parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm
CALL lodr(2,parts%R(blockstart:blockend),ra(n), rb(n))
parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm
END DO
IF(distribtype .eq. 5) THEN
! Initial velocities distribution
vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity
CALL lodgaus(3,VZ)
CALL lodgaus(5,VR)
CALL lodgaus(7,VTHET)
VZ=VZ*vth/4
VR=VR*8*vth
VTHET=VTHET*8*vth
ELSE
! Load velocities theta velocity
! Loading of r and z velocity is done in adapt_vinit to have
! access to parts%pot
DO i=1,parts%Nptot
! Interpolation for Magnetic potential
rg=parts%R(i)*rnorm
zg=(parts%Z(i)-halfLz)*rnorm
Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width))
Pcomp=P0/rg/me
Acomp=-SIGN(elchar/me*Athetpos,qsim)
VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp)
!VTHET(i)=Pcomp+Acomp
END DO
VTHET=VTHET/vnorm
VZ=0._db
VR=0._db
END IF
END SUBROUTINE loadDavidson
!===============================================================================
SUBROUTINE Temp_rescale
USE basic, ONLY: temprescale, nlclassical, nz
! INTEGER:: i, gridindex
! REAL(kind=db) :: vr, vz, vthet
! DO i=1,parts%Nptot
! gridindex=parts%zindex(i)+1+parts%rindex(i)*nz
! vr=parts%UR(i)/parts%Gamma(i)-fluidur(gridindex)
! vz=parts%UZ(i)/parts%Gamma(i)-fluiduz(gridindex)
! vthet=parts%UTHET(i)/parts%Gamma(i)-fluiduthet(gridindex)
! parts%UR(i)=vr*temprescale+fluidur(gridindex)
! parts%UTHET(i)=vthet*temprescale+fluiduthet(gridindex)
! parts%UZ(i)=vz*temprescale+fluiduz(gridindex)
! IF(nlclassical) THEN
! parts%Gamma(i)=1.0
! ELSE
! parts%Gamma(i)=1/sqrt(1-(parts%UR(i)**2+parts%UTHET(i)**2+parts%UZ(i)**2))
! END IF
! parts%UR(i)=parts%UR(i)*parts%Gamma(i)
! parts%UTHET(i)=parts%UTHET(i)*parts%Gamma(i)
! parts%UZ(i)=parts%UZ(i)*parts%Gamma(i)
! END DO
END SUBROUTINE Temp_rescale
!---------------------------------------------------------------------------
!> @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, qsim, msim, dt, dr, dz
INTEGER, INTENT(IN) ::samplefactor
TYPE(particles), INTENT(INOUT):: p
INTEGER:: i, j, currentindex
REAL(kind=db), DIMENSION(p%Nptot) :: 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(dz,minval(dr))/100
DO i=1,p%Nptot
DO j=1,samplefactor-1
currentindex=p%Nptot+(i-1)*(samplefactor-1)+j
CALL move_part(p,i,currentindex)
p%partindex(currentindex)=currentindex
dir = spreaddir(i)+2*pi*j/samplefactor
p%R(currentindex)=p%R(currentindex) + dl*cos(dir)
p%Z(currentindex)=p%Z(currentindex) + dl*sin(dir)
END DO
p%partindex(i)=i
p%R(i)=p%R(i) + dl*cos(spreaddir(i))
p%Z(i)=p%Z(i) + dl*sin(spreaddir(i))
END DO
nplasma=nplasma*samplefactor
qsim=qsim/samplefactor
msim=msim/samplefactor
p%Nptot=p%Nptot*samplefactor
END SUBROUTINE upsample
END MODULE beam
>>>>>>> 580212fe518b24d0c099fe2f77d86387df7612cc
c4science · Help