diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index 60f966f..afcd010 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,1588 +1,1615 @@ !------------------------------------------------------------------------------ ! 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 IMPLICIT NONE !> Stores the particles properties for the run. TYPE particles INTEGER :: Nptot !< Local number of simulated particles 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 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: R !< radial coordinates of the particles DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Z !< longitudinal coordinates of the particles DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: THET !< azimuthal coordinates of the particles DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: WZ !< axial radial relative distances to the left grid line DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: WR !< radial relative distances to the bottom grid line DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: pot !< Electric potential DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Er !< Radial Electric field DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Ez !< Axial electric field DOUBLE PRECISION, DIMENSION(:),POINTER:: UR !< normalized radial velocity at the current time step DOUBLE PRECISION, DIMENSION(:),POINTER:: URold !< normalized radial velocity at the previous time step DOUBLE PRECISION, DIMENSION(:),POINTER:: UTHET !< normalized azimuthal velocity at the current time step DOUBLE PRECISION, DIMENSION(:),POINTER:: UTHETold !< normalized azimuthal velocity at the previous time step DOUBLE PRECISION, DIMENSION(:),POINTER:: UZ !< normalized axial velocity at the current time step DOUBLE PRECISION, DIMENSION(:),POINTER:: UZold !< normalized axial velocity at the previous time step DOUBLE PRECISION, DIMENSION(:),POINTER:: Gamma !< Lorentz factor at the current time step DOUBLE PRECISION, 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) DOUBLE PRECISION :: ekin !< Total kinetic energz (J) DOUBLE PRECISION :: epot !< Total potential energy (J) DOUBLE PRECISION :: etot !< Current total energy (J) DOUBLE PRECISION :: etot0 !< Initial total energy (J) INTEGER, DIMENSION(7), 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)) parts%URold=0 parts%UZold=0 parts%UTHETold=0 parts%rindex=0 parts%zindex=0 parts%WR=0 parts%WZ=0 parts%UR=0 parts%UZ=0 parts%UTHET=0 parts%Z=0 parts%R=0 parts%THET=0 parts%Gamma=1 parts%Er=0 parts%Ez=0 parts%pot=0 parts%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, Rcurv USE mpi INTEGER:: i DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V 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 @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 USE basic, ONLY: zgrid, nz, Zbounds, cstep, mpirank, partperiodic INTEGER :: i, rsendnbparts=0, lsendnbparts=0 INTEGER, DIMENSION(parts%Nptot) :: sendhole sendhole=0 lsendnbparts=0 rsendnbparts=0 IF (parts%Nptot .NE. 0) THEN ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot IF(cstep .ne. 0) THEN ! If the particle is to the right of the local simulation space, it is sent to the right MPI process IF (parts%Z(i) .ge. Zbounds(mpirank+1)) THEN !$OMP CRITICAL (nbparts) rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i !$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 (parts%Z(i) .lt. Zbounds(mpirank)) THEN !$OMP CRITICAL (nbparts) lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i !$OMP END CRITICAL (nbparts) END IF END IF ! The periodic boundary conditions in z are here applied DO WHILE (parts%Z(i) .GT. zgrid(nz)) parts%Z(i) = parts%Z(i) - zgrid(nz) + zgrid(0) END DO DO WHILE (parts%Z(i) .LT. zgrid(0)) parts%Z(i) = parts%Z(i) + zgrid(nz) - zgrid(0) END DO END DO !$OMP END PARALLEL DO IF(mpisize .gt. 1) THEN ! We send the particles leaving the local simulation space to the closest neighbour CALL particlescommunication(lsendnbparts, rsendnbparts, sendhole) ELSE IF (.NOT. partperiodic) THEN ! If the boundary conditions are not periodic, we delete the corresponding particles IF(lsendnbparts+rsendnbparts .gt. 0) THEN DO i=rsendnbparts+lsendnbparts,1,-1 CALL delete_part(abs(sendhole(i))) END DO END IF END IF 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 USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm INTEGER :: j,i, nblostparts=0 INTEGER, DIMENSION(parts%Nptot) :: losthole i=1 losthole=0 IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot parts%zindex(i)=(parts%Z(i)-zgrid(0))/dz IF (parts%R(i) .GT. rgrid(0) .AND. parts%R(i) .LT. rgrid(nnr(1))) THEN parts%rindex(i)=(parts%R(i)-rgrid(0))/dr(1) ELSE IF(parts%R(i) .GT. rgrid(nnr(1)) .AND. parts%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)))/dr(2)+nnr(1) ELSE IF(parts%R(i) .GT. rgrid(nnr(1)+nnr(2)) .AND. parts%R(i) .LT. rgrid(nr)) THEN parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)+nnr(2)))/dr(3)+nnr(1)+nnr(2) ELSE ! If the particle is outside of the simulation space in the r direction, it is deleted. WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in r:", parts%R(i)*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Particle removed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" !$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 j=nblostparts,1,-1 CALL delete_part(losthole(j)) END DO END IF !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot parts%WZ(i)=(parts%Z(i)-zgrid(parts%zindex(i)))/(zgrid(parts%zindex(i)+1)-zgrid(parts%zindex(i))); parts%WR(i)=(parts%R(i)-rgrid(parts%rindex(i)))/(rgrid(parts%rindex(i)+1)-rgrid(parts%rindex(i))); END DO !$OMP END PARALLEL DO 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 ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : nlclassical ! Store old Velocities CALL swappointer(parts%UZold, parts%UZ) CALL swappointer(parts%URold, parts%UR) CALL swappointer(parts%UTHETold, parts%UTHET) CALL swappointer(parts%Gammaold, parts%Gamma) IF (nlclassical) THEN CALL comp_velocity_fun(gamma_classical) ELSE CALL comp_velocity_fun(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(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) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gam end subroutine end interface DOUBLE PRECISION :: tau DOUBLE PRECISION:: 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 (parts%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,parts%Nptot J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1 J2=J1+1 J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1 J4=J3+1 ! Interpolation for magnetic field BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) & & +parts%WZ(i)*parts%WR(i)*Bz(J1) BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Br(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Br(J2) & & +parts%WZ(i)*parts%WR(i)*Br(J1) ! First half of electric pulse parts%UZ(i)=parts%UZold(i)+parts%Ez(i)*tau parts%UR(i)=parts%URold(i)+parts%ER(i)*tau CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) ! Rotation along magnetic field ZBZ=tau*BRZ/parts%Gamma(i) ZBR=tau*BRR/parts%Gamma(i) ZPZ=parts%UZ(i)-ZBR*parts%UTHETold(i) !u'_{z} ZPR=parts%UR(i)+ZBZ*parts%UTHETold(i) !u'_{r} ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZ(i)-ZBZ*parts%UR(i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR parts%UZ(i)=parts%UZ(i)-ZBR2*ZPTHET !u+_{z} parts%UR(i)=parts%UR(i)+ZBZ2*ZPTHET !u+_{r} parts%UTHET(i)=parts%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration parts%UZ(i)=parts%UZ(i)+parts%EZ(i)*tau parts%UR(i)=parts%UR(i)+parts%ER(i)*tau ! Final computation of the Lorentz factor CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) END DO !$OMP END PARALLEL DO 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 !--------------------------------------------------------------------------- SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, 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 !--------------------------------------------------------------------------- SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, 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 Use basic, ONLY: dt, tnorm DOUBLE PRECISION:: XP, YP, COSA, SINA, U1, U2, ALPHA INTEGER :: i IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2, ALPHA) DO i=1,parts%Nptot ! Local Cartesian coordinates XP=parts%R(i)+dt/tnorm*parts%UR(i)/parts%Gamma(i) YP=dt/tnorm*parts%UTHET(i)/parts%Gamma(i) ! Conversion to cylindrical coordiantes parts%Z(i)=parts%Z(i)+dt/tnorm*parts%UZ(i)/parts%Gamma(i) parts%R(i)=sqrt(XP**2+YP**2) ! Computation of the rotation angle IF (parts%R(i) .EQ. 0) THEN COSA=1 SINA=0 ALPHA=0 ELSE COSA=XP/parts%R(i) SINA=YP/parts%R(i) ALPHA=asin(SINA) END IF ! New azimuthal position parts%THET(i)=parts%THET(i)+ALPHA ! Velocity in rotated reference frame U1=COSA*parts%UR(i)+SINA*parts%UTHET(i) U2=-SINA*parts%UR(i)+COSA*parts%UTHET(i) parts%UR(i)=U1 parts%UTHET(i)=U2 END DO !$OMP END PARALLEL DO 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, partdensity, ierr, step, it0d, it2d, nlend, nr, itparts, nplasma, fluidur, fluiduz, fluiduthet, Presstens, partperiodic, newres DOUBLE PRECISION :: gammamid, vr, vz, vthet INTEGER :: i, gridindex INTEGER, DIMENSION(:) :: stat(7*MPI_STATUS_SIZE) CALL MPI_WAITALL(7, ireducerequest, stat, ierr) ! Reset the quantities partdensity=0 fluidur=0 fluiduz=0 fluiduthet=0 ekin=0 epot=0 etot=0 ! Computation of the kinetic and potential energy as well as fluid velocities and density !$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity,fluidur,fluiduz,fluiduthet) DEFAULT(SHARED) PRIVATE(i,gammamid,vr,vz,vthet,gridindex) 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 IF(modulo(step,it2d).eq. 0 .or. nlend) THEN gridindex=parts%zindex(i)+1+parts%rindex(i)*nz ! Number of particles per cell partdensity(gridindex)=partdensity(gridindex)+1 ! Individual particle velocities vr=parts%UR(i)/parts%Gamma(i) vz=parts%UZ(i)/parts%Gamma(i) vthet=parts%UTHET(i)/parts%Gamma(i) ! Fluid velocities fluidur(gridindex)=fluidur(gridindex)+vr fluiduz(gridindex)=fluiduz(gridindex)+vz fluiduthet(gridindex)=fluiduthet(gridindex)+vthet END IF END DO !$OMP END PARALLEL DO ! The computed energy is sent to the root process 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) 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) END IF ! Normalize fluid velocities by density avoiding inf values WHERE (partdensity.gt.0) fluiduz=fluiduz/partdensity fluidur=fluidur/partdensity fluiduthet=fluiduthet/partdensity END WHERE ! Computes the pressure tensor Presstens=0 IF(modulo(step,it2d).eq. 0 .or. nlend) THEN !$OMP PARALLEL DO REDUCTION(+: Presstens) DEFAULT(SHARED) PRIVATE(i,vr,vz,vthet,gridindex) 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) Presstens(1,gridindex)=Presstens(1,gridindex)+vr*vr Presstens(2,gridindex)=Presstens(2,gridindex)+vr*vthet Presstens(3,gridindex)=Presstens(3,gridindex)+vr*vz Presstens(4,gridindex)=Presstens(4,gridindex)+vthet*vthet Presstens(5,gridindex)=Presstens(5,gridindex)+vthet*vz Presstens(6,gridindex)=Presstens(6,gridindex)+vz*vz END DO !$OMP END PARALLEL DO !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,6 WHERE (partdensity.gt.0) Presstens(i,:)=Presstens(i,:)/partdensity END DO !$OMP END PARALLEL DO END IF ! Send the particle density and fluid velocities from the workers to the master IF(modulo(step,it2d).eq. 0 .or. nlend) THEN IF(mpirank .eq.0 ) THEN CALL MPI_IREDUCE(MPI_IN_PLACE, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) CALL MPI_IREDUCE(MPI_IN_PLACE, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(4), ierr) CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(5), ierr) CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduthet, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(6), ierr) CALL MPI_IREDUCE(MPI_IN_PLACE, Presstens, 6*nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(7), ierr) ELSE CALL MPI_IREDUCE(partdensity, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) CALL MPI_IREDUCE(fluidur, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(4), ierr) CALL MPI_IREDUCE(fluiduz, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(5), ierr) CALL MPI_IREDUCE(fluiduthet, fluiduthet, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(6), ierr) CALL MPI_IREDUCE(Presstens, Presstens, 6*nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(7), ierr) END IF END IF ! Wait for the communication of the energy to be finished CALL MPI_WAIT(ireducerequest(1), stat(1: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 ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ.1 .or. (newres .AND. step .EQ. 1)) THEN etot0 = etot 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 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 USE basic, ONLY: mpirank, mpisize, ierr 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(parts%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%THET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%THET, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%Rindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%Zindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%partindex, Nplocs_all(mpirank), MPI_INTEGER, parts%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, parts%Z, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%THET, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%partindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) 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 !! 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 DOUBLE PRECISION :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, ZBZ2, ZBR2, Vperp, v2 INTEGER :: J1, J2, J3, J4, i DOUBLE PRECISION, 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(parts%Nptot),VZ(parts%Nptot),VTHET(parts%Nptot)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=parts%UTHET/parts%Gamma*vnorm DO i=1,parts%Nptot Vperp=sqrt(MAX(2*H0/me+2*elchar/me*parts%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN parts%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 parts%Gamma(i)=sqrt(1/(1-v2/vnorm**2)) END IF parts%UR(i)=parts%Gamma(i)*VR(i)/vnorm parts%UZ(i)=parts%Gamma(i)*VZ(i)/vnorm parts%UTHET(i)=parts%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(parts%UZold, parts%UZ) CALL swappointer(parts%URold, parts%UR) CALL swappointer(parts%UTHETold, parts%UTHET) CALL swappointer(parts%Gammaold, parts%Gamma) IF (parts%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,parts%Nptot ! Compute the particle linear indices for the magnetic field interpolation J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1 J2=J1+1 J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1 J4=J3+1 ! Interpolation for magnetic field BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) & & +parts%WZ(i)*parts%WR(i)*Bz(J1) BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Br(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Br(J2) & & +parts%WZ(i)*parts%WR(i)*Br(J1) ! Half inverse Rotation along magnetic field ZBZ=tau*BRZ/parts%Gammaold(i) ZBR=tau*BRR/parts%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(parts%UZold(i)-ZBR*parts%UTHETold(i))/SQR !u-_{z} ZPR=(parts%URold(i)+ZBZ*parts%UTHETold(i))/SQR !u-_{r} ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZold(i)-ZBZ*parts%URold(i))/SQR !u-_{theta} parts%UZ(i)=ZPZ parts%UR(i)=ZPR parts%UTHET(i)=ZPTHET ! half of decceleration parts%UZ(i)=parts%UZ(i)+parts%Ez(i)*tau parts%UR(i)=parts%UR(i)+parts%Er(i)*tau IF(.not. nlclassical) THEN parts%Gamma(i)=sqrt(1+parts%UZ(i)**2+parts%UR(i)**2+parts%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, nplasma, ierr, Zbounds INTEGER:: k DOUBLE PRECISION:: 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 partspercol=0 idealnbpartsperproc = FLOOR(REAL(parts%Nptot)/REAL(mpisize)) Zmin=MINVAL(parts%Zindex) Zmax=MAXVAL(parts%Zindex) Zperproc=(Zmax-Zmin)/mpisize 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 ! 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) THEN DO k= partindexstart, partindexlast CALL move_part(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(lsendnbparts, rsendnbparts, sendholes) USE mpihelper, ONLY: particle_type USE basic, ONLY: rightproc, leftproc, ierr INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) 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(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(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(rrecvnbparts, lrecvnbparts, sendnbparts, sendholes) ! USE mpihelper INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) 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 parts%Nptot=parts%Nptot+1 partpos=parts%Nptot END IF CALL Insertincomingpart(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 parts%Nptot=parts%Nptot+1 partpos=parts%Nptot END IF CALL Insertincomingpart(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%Nptot,abs(sendholes(k))) parts%Nptot = parts%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(buffer, bufferindex, partsindex) USE mpihelper INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer parts%partindex(partsindex) = buffer(bufferindex)%partindex parts%Rindex(partsindex) = buffer(bufferindex)%Rindex parts%Zindex(partsindex) = buffer(bufferindex)%Zindex parts%R(partsindex) = buffer(bufferindex)%R parts%Z(partsindex) = buffer(bufferindex)%Z parts%THET(partsindex) = buffer(bufferindex)%THET parts%UZ(partsindex) = buffer(bufferindex)%UZ parts%UR(partsindex) = buffer(bufferindex)%UR parts%UTHET(partsindex) = buffer(bufferindex)%UTHET parts%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(buffer, bufferindex, partsindex) USE mpihelper INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer buffer(bufferindex)%partindex = parts%partindex(partsindex) buffer(bufferindex)%Rindex = parts%Rindex(partsindex) buffer(bufferindex)%Zindex = parts%Zindex(partsindex) buffer(bufferindex)%R = parts%R(partsindex) buffer(bufferindex)%Z = parts%Z(partsindex) buffer(bufferindex)%THET = parts%THET(partsindex) buffer(bufferindex)%UZ = parts%UZ(partsindex) buffer(bufferindex)%UR = parts%UR(partsindex) buffer(bufferindex)%UTHET = parts%UTHET(partsindex) buffer(bufferindex)%Gamma = parts%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(lsendnbparts, rsendnbparts, sendholes) ! USE mpihelper INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) 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(rsendpartbuff, rsendpos, partpos) ELSE IF(sendholes(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(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(index1, index2) INTEGER, INTENT(IN) :: index1, index2 DOUBLE PRECISION:: 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= parts%partindex(index1) Gamma = parts%Gamma(index1) R = parts%R(index1) Z = parts%Z(index1) THET = parts%THET(index1) UR = parts%UR(index1) UTHET = parts%UTHET(index1) UZ = parts%UZ(index1) Rindex = parts%Rindex(index1) Zindex = parts%Zindex(index1) ! Move part at index2 in part at index 1 parts%partindex(index1)= parts%partindex(index2) parts%Gamma(index1) = parts%Gamma(index2) parts%R(index1) = parts%R(index2) parts%Z(index1) = parts%Z(index2) parts%THET(index1) = parts%THET(index2) parts%UR(index1) = parts%UR(index2) parts%UTHET(index1) = parts%UTHET(index2) parts%UZ(index1) = parts%UZ(index2) parts%Rindex(index1) = parts%Rindex(index2) parts%Zindex(index1) = parts%Zindex(index2) ! Move temporary values from part(index1) to part(index2) parts%partindex(index2)= partindex parts%Gamma(index2) = Gamma parts%R(index2) = R parts%Z(index2) = Z parts%THET(index2) = THET parts%UR(index2) = UR parts%UTHET(index2) = UTHET parts%UZ(index2) = UZ parts%Rindex(index2) = Rindex parts%Zindex(index2) = Zindex END SUBROUTINE exchange_parts +!--------------------------------------------------------------------------- +!> @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(index) + !! This will destroy particle at index + INTEGER, INTENT(IN) :: index + + etot0=etot0-qsim*parts%pot(index)*phinorm + IF(.not. nlclassical) THEN + etot0=etot0-msim*vlight**2*(parts%Gamma(index)-1) + ELSE + etot0=etot0-0.5*msim*vlight**2*(abs(parts%URold(index)*parts%UR(index))+abs(parts%UZold(index)*parts%UZ(index))+abs(parts%UTHETold(index)*parts%UTHET(index))) + END IF + + ! We fill the gap + CALL move_part(parts%Nptot, index) + + ! Reduce the total number of simulated parts + parts%Nptot=parts%Nptot-1 + 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 destroy 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(sourceindex, destindex) - !! !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex ! Move part at sourceindex in part at destindex parts%partindex(destindex)= parts%partindex(sourceindex) parts%Gamma(destindex) = parts%Gamma(sourceindex) parts%R(destindex) = parts%R(sourceindex) parts%Z(destindex) = parts%Z(sourceindex) parts%THET(destindex) = parts%THET(sourceindex) parts%UR(destindex) = parts%UR(sourceindex) parts%UTHET(destindex) = parts%UTHET(sourceindex) parts%UZ(destindex) = parts%UZ(sourceindex) parts%Rindex(destindex) = parts%Rindex(sourceindex) parts%Zindex(destindex) = parts%Zindex(sourceindex) END SUBROUTINE move_part !--------------------------------------------------------------------------- !> @author !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load a uniform distribution using the Hammersley's sequence (0<=y<=1). !> (nbase = 0 => Random sampling) ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE loduni(nbase, y) ! ! Load a uniform distribution using the Hammersley's sequence. ! (nbase = 0 => Random sampling) ! INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER :: n, i ! n = SIZE(y) ! SELECT CASE (nbase) CASE(0) CALL RANDOM_NUMBER(y) CASE(1) DO i=1,n y(i) = (i-0.5_db)/n END DO CASE(2:) DO i=1,n y(i) = rev(nbase,i) END DO CASE default WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase END SELECT ! END SUBROUTINE loduni !--------------------------------------------------------------------------- !> @author !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load a gaussian distribution using a uniform distribution according to the Hammersley's sequence. !> (nbase = 0 => Random sampling) ! !> @param [in] nbase Base of the Hammersley's sequence. (0 => Random) !> @param [out] y array containing the random variables. !--------------------------------------------------------------------------- SUBROUTINE lodgaus(nbase, y) ! ! Sample y from the Gauss distributrion ! DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase ! INTEGER :: n, i, iflag DOUBLE PRECISION :: r(SIZE(y)) DOUBLE PRECISION :: t ! DOUBLE PRECISION :: c0, c1, c2 DOUBLE PRECISION :: d1, d2, d3 DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/ DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/ ! n = SIZE(y) CALL loduni(nbase, r) r = 1.0E-5_db + 0.99998_db*r ! DO i=1,n iflag = 1 IF (r(i) .GT. 0.5_db) THEN r(i) = 1.0_db - r(i) iflag = -1 END IF t = SQRT(LOG(1.0_db/(r(i)*r(i)))) y(i) = t - (c0+c1*t+c2*t**2) / (1.0_db+d1*t+d2*t**2+d3*t**3) y(i) = y(i) * iflag END DO y = y - SUM(y)/REAL(n,db) END SUBROUTINE lodgaus !________________________________________________________________________________ SUBROUTINE lodinvr(nbase, y, ra, rb) ! ! Sample y from the distribution f=1/(r)H(r-ra)H(rb-r) for where H is the heavyside function ! DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra ! DOUBLE PRECISION :: r(SIZE(y)) ! CALL loduni(nbase, r) r=r*log(rb/ra) y=ra*exp(r) END SUBROUTINE lodinvr !________________________________________________________________________________ SUBROUTINE lodlinr(nbase, y, ra, rb) DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra ! DOUBLE PRECISION :: r(SIZE(y)) ! CALL loduni(nbase, r) y=ra+sqrt(r)*(rb-ra) END SUBROUTINE lodlinr !________________________________________________________________________________ SUBROUTINE lodunir(nbase, y, ra, rb) DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra ! DOUBLE PRECISION :: r(SIZE(y)) ! CALL loduni(nbase, r) y=ra+(rb-ra)*r END SUBROUTINE lodunir !________________________________________________________________________________ SUBROUTINE lodgausr(nbase, y, ra, rb) DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra ! DOUBLE PRECISION :: r(SIZE(y)) ! CALL lodgaus(nbase, r) y=0.5*(rb+ra)+ 0.1*(rb-ra)*r END SUBROUTINE lodgausr !________________________________________________________________________________ REAL(db) FUNCTION rev(nbase,i) ! ! Return an element of the Hammersley's sequence ! INTEGER, INTENT(IN) :: nbase ! Base of the Hammersley's sequence (IN) INTEGER, INTENT(IN) :: i ! Index of the sequence (IN) ! ! Local vars INTEGER :: j1, j2 REAL(db) :: xs, xsi !----------------------------------------------------------------------- xs = 0._db xsi = 1.0_db j2 = i DO xsi = xsi/nbase j1 = j2/nbase xs = xs + (j2-nbase*j1)*xsi j2 = j1 IF( j2.LE.0 ) EXIT END DO rev = xs END FUNCTION rev !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) DOUBLE PRECISION :: bessi1,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 !________________________________________________________________________________ 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 !________________________________________________________________________________ RECURSIVE SUBROUTINE quicksortparts(leftlimit, rightlimit) ! Sorts the particle according to their Z position using quicksort algorithm INTEGER,INTENT(IN):: leftlimit, rightlimit DOUBLE PRECISION:: pivot INTEGER::i, cnt, mid IF(leftlimit .ge. rightlimit) RETURN mid=(leftlimit+rightlimit)/2 IF(parts%Z(mid).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,mid) IF(parts%Z(rightlimit).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,rightlimit) IF(parts%Z(mid).lt.parts%Z(rightlimit)) CALL exchange_parts(rightlimit,mid) ! Store the pivot point for comparison pivot=parts%Z(rightlimit) cnt=leftlimit ! Move all parts with Z smaller than pivot to the left of pivot DO i=leftlimit, rightlimit IF(parts%Z(i) .le. pivot) THEN CALL exchange_parts(i,cnt) cnt=cnt+1 END IF END DO CALL quicksortparts(leftlimit,cnt-2) CALL quicksortparts(cnt,rightlimit) END SUBROUTINE quicksortparts !________________________________________________________________________________ SUBROUTINE swappointer( pointer1, pointer2) DOUBLE PRECISION, DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2 DOUBLE PRECISION, 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 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET DOUBLE PRECISION:: vth ! Initial distribution in z with normalisation CALL loduni(1,parts%Z(:)) parts%Z(:)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(:))/rnorm ! Initial distribution in r with normalisation CALL loduni(2,parts%R(:)) parts%R=(plasmadim(3)+parts%R(:)*(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) DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra end subroutine end interface DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE::ra, rb DOUBLE PRECISION :: 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)) 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 END MODULE beam