diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index 954ddf2..0eba63e 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,1965 +1,1936 @@ MODULE beam !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Guillaume Le Bars EPFL/SPC !> Patryk Kaminski EPFL/SPC !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> Module responsible for loading, advancing and computing the necessary diagnostics for the simulated particles. !------------------------------------------------------------------------------ ! USE constants use mpi USE mpihelper USE basic, ONLY: mpirank, mpisize USE distrib USE particletypes USE weighttypes IMPLICIT NONE ! !TYPE(particles) :: parts !< Storage for all the particles !SAVE :: parts TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist ! Diagnostics (scalars) REAL(kind=db) :: ekin=0 !< Total kinetic energy (J) REAL(kind=db) :: epot=0 !< Total potential energy (J) REAL(kind=db) :: etot=0 !< Current total energy (J) REAL(kind=db) :: etot0=0 !< Initial total energy (J) REAL(kind=db) :: loc_etot0=0 !< theoretical local total energy (J) REAL(kind=db) :: Energies(4) !< (1) kinetic energy, (2) potential energy, (3) total energy and (4) gained/lossed energy due to gain or loss of particles (J) ! INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process INTERFACE add_created_part MODULE PROCEDURE add_linked_created_part, add_list_created_part END INTERFACE add_created_part ! abstract interface subroutine rloader(nbase,y,rminus,rplus) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rplus, rminus end subroutine REAL(kind=db) FUNCTION gamma(UZ, UR, UTHET) USE constants REAL(kind=db), INTENT(IN):: UR,UZ,UTHET end FUNCTION end interface CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Loads the particles at the beginning of the simulation and create the parts variable if necessary !--------------------------------------------------------------------------- SUBROUTINE load_parts USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical, nbspecies, Zbounds, partfile use mpi INTEGER:: i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V, uniform in Z and 1/R in R CALL loaduniformRZ(partslist(1), VR, VZ, VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(3) !Stable distribution from Davidson 4.95 p.119 but with constant distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodinvr) CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodgausr) CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz @brief Checks for each particle if the z position is outside of the local/global simulation space. !> Depending on the boundary conditions, the leaving particles are sent to the correct neighbouring MPI process !> or deleted. ! !> @param[in] p particles structure ! !> @author Guillaume Le Bars EPFL/SPC !--------------------------------------------------------------------------- SUBROUTINE bound(p) USE basic, ONLY: zgrid, nz, Zbounds, mpirank, step, leftproc, rightproc, partperiodic IMPLICIT NONE type(particles), INTENT(INOUT):: p INTEGER :: i, rsendnbparts, lsendnbparts, nblostparts INTEGER :: receivednbparts, partdiff INTEGER, DIMENSION(p%Nploc) :: sendhole INTEGER, DIMENSION(p%Nploc) :: losthole LOGICAL:: leftcomm, rightcomm INTEGER, ALLOCATABLE:: partstoremove(:) receivednbparts=0 nblostparts=0 rsendnbparts=0 lsendnbparts=0 IF (p%Nploc .gt. 0) THEN losthole=0 sendhole=0 ! We communicate with the left processus leftcomm = leftproc .ne. -1 ! We communicate with the right processus rightcomm = rightproc .ne. -1 ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) DO i=1,p%Nploc ! 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 IF(partperiodic) THEN DO WHILE (p%Z(i) .GT. zgrid(nz)) p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0) END DO END IF !$OMP CRITICAL (nbparts) IF(rightcomm) THEN rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i ELSE nblostparts=nblostparts+1 losthole(nblostparts)=i p%nblost(2)=p%nblost(2)+1 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 IF(partperiodic) THEN DO WHILE (p%Z(i) .LT. zgrid(0)) p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0) END DO END IF !$OMP CRITICAL (nbparts) IF(leftcomm) THEN ! We send the particle to the left process lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i ELSE ! we destroy the particle nblostparts=nblostparts+1 losthole(nblostparts)=i p%nblost(1)=p%nblost(1)+1 END IF !$OMP END CRITICAL (nbparts) END IF END DO !$OMP END PARALLEL DO END IF IF(mpisize .gt. 1) THEN ! We send the particles leaving the local simulation space to the closest neighbour CALL particlescommunication(p, lsendnbparts, rsendnbparts, sendhole, receivednbparts, (/leftproc,rightproc/)) END IF ! If the boundary conditions are not periodic, we delete the corresponding particles IF(nblostparts .gt. 0 .and. step .ne. 0) THEN DO i=1,nblostparts CALL delete_part(p, losthole(i), .false. ) END DO !WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank END IF ! computes if we received less particles than we sent partdiff=max(lsendnbparts+rsendnbparts-receivednbparts,0) IF(nblostparts + partdiff .gt. 0) THEN ALLOCATE(partstoremove(nblostparts+partdiff)) partstoremove(1:partdiff)=abs(sendhole(receivednbparts+1:receivednbparts+partdiff)) partstoremove(partdiff+1:partdiff+nblostparts)=abs(losthole(1:nblostparts)) call LSDRADIXSORT(partstoremove,size(partstoremove)) ! If we received less particles than we sent, or lost particles we fill the remaining holes with the particles from the end of the parts arrays DO i=nblostparts+partdiff,1,-1 CALL move_part(p, p%Nploc, partstoremove(i)) p%partindex(p%Nploc)=-1 p%Nploc = p%Nploc-1 END DO END IF END subroutine bound !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Check if a particle is outside the simulation domain and remove it if needed !> @param[in] p particles structure !--------------------------------------------------------------------------- SUBROUTINE boundary_loss(p) USE basic, ONLY: rgrid, nr Use geometry, ONLY: r_a, geom_weight, dom_weight type(particles), INTENT(INOUT):: p INTEGER :: i,j,isup, nblostparts, iend,nbunch INTEGER, DIMENSION(p%Nploc) :: losthole INTEGER, DIMENSION(16)::idwall INTEGER :: nblost(size(p%nblost,1)) nblostparts=0 nblost=0 nbunch=16 IF (p%Nploc .le. 0) return losthole=0 !$OMP PARALLEL DEFAULT(SHARED), private(i,iend,j,isup,idwall) !$OMP DO reduction(+:nblost) DO i=1,p%Nploc,nbunch ! Avoid segmentation fault caused by accessing non relevant data iend=min(i+nbunch-1,p%Nploc) ! calculate the weight do determine if a particle is inside the simulation domain. call dom_weight(p%Z(i:iend), p%r(i:iend), p%geomweight(i:iend,0),idwall) do j=i,iend if(p%geomweight(j,0).le.0 .or. p%R(j) .ge. rgrid(nr) .or. p%R(j) .le. rgrid(0)) then ! If the particle is outside of the simulation space in the r direction, or if it is outside of the vacuum region it is deleted. !$OMP CRITICAL (lostparts) nblostparts=nblostparts+1 losthole(nblostparts)=j !$OMP END CRITICAL (lostparts) isup=0 if(p%R(j) .ge. rgrid(nr) .or. idwall(j-i+1) .gt.0) then isup=1 end if nblost(3+isup+idwall(j-i+1))=nblost(3+isup+idwall(j-i+1))+1 + else + call p_calc_rzindex(p,j) + call geom_weight(p%Z(j), p%r(j), p%geomweight(j,:)) end if end do END DO !$OMP END DO !$OMP END PARALLEL IF(nblostparts.gt.0) THEN p%nblost=nblost+p%nblost !call qsort(losthole,p%Nploc,sizeof(losthole(1)),compare_int) call LSDRADIXSORT(losthole(1:nblostparts),nblostparts) !Write(*,'(a,60i)') "losthole: ", losthole(nblostparts:nblostparts+1) DO i=nblostparts,1,-1 CALL delete_part(p,losthole(i)) END DO END IF END SUBROUTINE boundary_loss - !--------------------------------------------------------------------------- - !> @author - !> Guillaume Le Bars EPFL/SPC - ! - ! DESCRIPTION: - !> @brief Compute the grid cell indices for each particle as well as the distance weight Wr, Wz. - !> @param[in] p particles structure - !--------------------------------------------------------------------------- -SUBROUTINE localisation(p) - USE basic, ONLY: rgrid, nr - Use geometry, ONLY: r_a, geom_weight, dom_weight - type(particles), INTENT(INOUT):: p - INTEGER :: i, j, iend,nbunch - nbunch=16 - - IF (p%Nploc .le. 0) return - !$OMP PARALLEL DEFAULT(SHARED), private(i,iend,j) - !$OMP DO - DO i=1,p%Nploc,nbunch - ! Avoid segmentation fault by accessing non relevant data - iend=min(i+nbunch-1,p%Nploc) - do j=i,iend - call p_calc_rzindex(p,j) - end do - call geom_weight(p%Z(i:iend), p%r(i:iend), p%geomweight(i:iend,:)) - END DO - !$OMP END DO - !$OMP END PARALLEL -END SUBROUTINE localisation - subroutine p_calc_rzindex(p,i) use basic, only: rgrid,zgrid,invdz,invdr, nnr, nr, nsubr integer::i,j,k type(particles)::p k=0 do j=1,nsubr IF (p%R(i) .GT. rgrid(k) .AND. p%R(i) .LT. rgrid(k+nnr(j))) THEN p%rindex(i)=floor((p%R(i)-rgrid(k))*invdr(j))+k exit end if k=k+nnr(j) end do !ELSE IF(p%R(i) .GE. rgrid(nnr(1)) .AND. p%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN ! p%rindex(i)=floor((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)=floor((p%R(i)-rgrid(nnr(1)+nnr(2)))*invdr(3))+nnr(1)+nnr(2) !End if p%zindex(i)=floor((p%Z(i)-zgrid(0))*invdz) end subroutine p_calc_rzindex SUBROUTINE comp_mag_p(p) USE basic, ONLY: zgrid, rgrid, BZ, BR, nz, invdz type(particles), INTENT(INOUT):: p INTEGER :: i Real(kind=db):: WZ,WR INTEGER:: j1,j2,j3,j4 !$OMP PARALLEL DO SIMD DEFAULT(SHARED) Private(J1,J2,J3,J4,WZ,WR) DO i=1,p%Nploc WZ=(p%Z(i)-zgrid(p%zindex(i)))*invdz; WR=(p%R(i)-rgrid(p%rindex(i)))/(rgrid(p%rindex(i)+1)-rgrid(p%rindex(i))); 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 p%BZ(i)=(1-WZ)*(1-WR)*Bz(J4) & & +WZ*(1-WR)*Bz(J3) & & +(1-WZ)*WR*Bz(J2) & & +WZ*WR*Bz(J1) p%BR(i)=(1-WZ)*(1-WR)*Br(J4) & & +WZ*(1-WR)*Br(J3) & & +(1-WZ)*WR*Br(J2) & & +WZ*WR*Br(J1) END DO !$OMP END PARALLEL DO SIMD end subroutine comp_mag_p !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the classical simulations. !> This routine systematically returns 1.0 to treat the system according to classical dynamic. ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- REAL(kind=db) FUNCTION gamma_classical(UZ, UR, UTHET) #if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_classical) #endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_classical=1.0 END FUNCTION gamma_classical !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the relativistic simulations. !> This routine computes the Lorentz factor \f$\gamma=\sqrt{1+\mathbf{\gamma\beta}^2}\f$ ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- REAL(kind=db) FUNCTION gamma_relativistic(UZ, UR, UTHET) #if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_relativistic) #endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_relativistic=sqrt(1+UZ**2+UR**2+UTHET**2) END FUNCTION gamma_relativistic !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief General routine to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity(p) ! ! 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$ !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity_fun(p, gammafun) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : bnorm, dt, tnorm procedure(gamma)::gammafun 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=p%qmRatio*bnorm*tnorm*0.5*dt/tnorm tau=p%qmRatio*bnorm*0.5*dt*tnorm IF (p%Nploc .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%Nploc ! 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 p%Gamma(i)=gammafun(p%UZ(i), p%UR(i), p%UTHETold(i)) ! Rotation along magnetic field ZBZ=tau*p%BZ(i)/p%Gamma(i) ZBR=tau*p%BR(i)/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 !p%ur(i)=0.001 ! Final computation of the Lorentz factor p%Gamma(i)=gammafun(p%UZ(i), p%UR(i), p%UTHET(i)) END DO !$OMP END PARALLEL DO SIMD END IF p%collected=.false. END SUBROUTINE comp_velocity_fun !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles position at time t+1 !> This routine computes the particles position at time t+1 according to the Bunemann algorithm. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE push(p) Use basic, ONLY: dt, tnorm type(particles), INTENT(INOUT):: p REAL(kind=db):: XP, YP, COSA, SINA, U1, U2, ALPHA INTEGER :: i IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(XP, YP, COSA, SINA, U1, U2, ALPHA) DO i=1,p%Nploc ! Local Cartesian coordinates XP=p%R(i)+dt*p%UR(i)/p%Gamma(i) YP=dt*p%UTHET(i)/p%Gamma(i) ! Conversion to cylindrical coordiantes p%Z(i)=p%Z(i)+dt*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 p%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 keeps track of the reference energy and the number of particle per mpi node. ! !--------------------------------------------------------------------------- SUBROUTINE partdiagnostics ! ! Compute energies ! USE constants, ONLY: vlight USE basic, ONLY: phinorm, cstep, nlclassical, ierr, step, nlend,& & itparts, nbspecies INTEGER:: i,j ! 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 REDUCTION(+:epot, ekin) DEFAULT(SHARED), PRIVATE(i,j) Do j=1,nbspecies if(.not. partslist(j)%is_field) CYCLE DO i=1,partslist(j)%Nploc ! Potential energy epot=epot+(partslist(j)%pot(i)+partslist(j)%potxt(i))*partslist(j)%q*partslist(j)%weight ! Kinetic energy IF(.not. nlclassical) THEN ekin=ekin+(0.5*(partslist(j)%Gammaold(i)+partslist(j)%Gamma(i))-1)*partslist(j)%m*partslist(j)%weight ELSE ekin=ekin+0.5*( partslist(j)%UR(i)*partslist(j)%URold(i) & & + partslist(j)%UZ(i)*partslist(j)%UZold(i) & & + partslist(j)%UTHET(i)*partslist(j)%UTHETold(i) )*partslist(j)%m*partslist(j)%weight END IF END DO END DO !$OMP END PARALLEL DO epot=epot*phinorm*0.5 ekin=ekin*vlight**2 ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ. 1) THEN ! Compute the local total energy loc_etot0 = epot+ekin etot0=0 END IF !etot=loc_etot0 ! Compute the total energy etot=epot+ekin Energies=(/ekin,epot,etot,loc_etot0/) ! The computed energy is sent to the root process IF(mpisize .gt.1) THEN IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) etot0=etot0+Energies(4) ekin=Energies(1) epot=Energies(2) etot=Energies(3) ELSE CALL MPI_REDUCE(Energies, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) END IF ELSE etot0=etot0+loc_etot0 END IF loc_etot0=0 ! Send the local number of particles on each node to the root process IF(mpisize .gt. 1) THEN Nplocs_all(mpirank)=partslist(1)%Nploc IF(mpirank .eq.0 ) THEN CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) !CALL MPI_REDUCE(MPI_IN_PLACE,partslist(1)%nudcol,3,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) partslist(1)%Nptot=sum(Nplocs_all) !partslist(1)%nudcol=partslist(1)%nudcol/partslist(1)%Nptot ELSE CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) !CALL MPI_REDUCE(partslist(1)%nudcol,partslist(1)%nudcol,3,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) END IF ELSE partslist(1)%Nptot=partslist(1)%Nploc END IF end subroutine partdiagnostics !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Collect the particles positions and velocities on the root process. !> If the collection has already been performed at this time step, the routine does nothing. ! !--------------------------------------------------------------------------- SUBROUTINE collectparts(p) USE basic, ONLY: mpirank, mpisize, ierr type(particles), INTENT(INOUT):: p INTEGER, DIMENSION(:), ALLOCATABLE :: displs, Nploc INTEGER:: i INTEGER:: particles_type(mpisize-1) !< Stores the MPI data type used for particles gathering on node 0 and broadcast from node 0 INTEGER :: part_requests(mpisize-1) INTEGER:: stats(MPI_STATUS_SIZE,mpisize-1) part_requests=MPI_REQUEST_NULL particles_type=MPI_DATATYPE_NULL IF(p%collected) RETURN ! exit subroutine if particles have already been collected during this time step ALLOCATE(Nploc(0:mpisize-1)) ALLOCATE(displs(0:mpisize-1)) displs=0 Nploc(mpirank)=p%Nploc CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, Nploc, 1, MPI_INTEGER,& & MPI_COMM_WORLD, ierr) p%Nptot=sum(Nploc) IF(p%Nptot .eq. 0 ) THEN p%collected=.true. RETURN END IF Do i=1,mpisize-1 displs(i)=displs(i-1)+Nploc(i-1) END DO IF(mpirank.eq.0 .and. p%Nptot .gt. size(p%R,1)) THEN CALL change_parts_allocation(p,max(p%Nptot-size(P%R,1),floor(0.5*size(P%R,1)))) END IF IF(mpirank .ne. 0) THEN if(Nploc(mpirank) .gt. 0) THEN Call init_particles_gather_mpi(p,1,Nploc(mpirank),particles_type(mpirank)) ! Send Particles informations to root process CALL MPI_SEND(p, 1, particles_type(mpirank), 0, partsgather_tag, MPI_COMM_WORLD, ierr) CALL MPI_TYPE_FREE(particles_type(mpirank),ierr) END IF ELSE ! Receive particle information from all processes DO i=1,mpisize-1 if(Nploc(i) .lt. 1) cycle Call init_particles_gather_mpi(p,displs(i)+1,Nploc(i),particles_type(i)) CALL MPI_IRECV(p,1,particles_type(i),i,partsgather_tag,MPI_COMM_WORLD, part_requests(i), ierr) END DO CALL MPI_WAITALL(mpisize-1,part_requests, stats, ierr) p%partindex(sum(Nploc)+1:)=-1 Do i=1,mpisize-1 if(Nploc(i) .lt. 1) cycle CALL MPI_TYPE_FREE(particles_type(i),ierr) END DO END IF p%collected=.TRUE. END SUBROUTINE collectparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Computes the velocities at time t-1/2 delta t to keep the second order precision in time on the velocity. !> This should only be used at particle initialisation time, ot in the case of a restart. ! !--------------------------------------------------------------------------- SUBROUTINE adapt_vinit(p) !! Computes the velocity at time -dt/2 from velocities computed at time 0 ! USE basic, ONLY : bnorm, dt, tnorm, nlclassical, phinorm, distribtype, vnorm type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, Vperp, v2 INTEGER :: 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 .or. p%Davidson) THEN ALLOCATE(VR(p%Nploc),VZ(p%Nploc),VTHET(p%Nploc)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=p%UTHET/p%Gamma*vnorm DO i=1,p%Nploc Vperp=sqrt(MAX(2*p%H0/p%m-2*p%qmRatio*p%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN p%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 p%Gamma(i)=sqrt(1/(1-v2/vnorm**2)) END IF p%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 tau=-p%qmRatio*bnorm*0.5*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%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR) DO i=1,p%Nploc ! Half inverse Rotation along magnetic field ZBZ=tau*p%BZ(i)/p%Gammaold(i) ZBR=tau*p%BR(i)/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 SIMD END IF END SUBROUTINE adapt_vinit !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Calculates the number of paticles per column of the spatial grid ( at fixed axial cell position) !> This facilitate the computation of the axial grid limits for each MPI worker ! !--------------------------------------------------------------------------- SUBROUTINE calcnbperz(p,nbperz) USE basic, only: nz IMPLICIT NONE type(particles):: p INTEGER, INTENT(INOUT):: nbperz(0:) Integer::i, zindex nbperz=0 !$OMP PARALLEL DO DEFAULT(SHARED) reduction(+:nbperz), private(zindex,i) Do i=1,p%Nploc ! we make sure zindex is in [0, nz-1] to avoid segmentation faults zindex=min(nz-1,max(p%zindex(i),0)) nbperz(zindex)=nbperz(zindex)+1 END DO !$OMP END PARALLEL DO END SUBROUTINE calcnbperz !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief In the case of MPI parallelism, computes the axial limits of the local domain. !--------------------------------------------------------------------------- SUBROUTINE calc_Zbounds(p, Zbounds, norder) ! Computes the start and end indices for the Z boundaries on local processus ! Computes the particle indices from initial particle loading vector, that stay in current process USE basic, ONLY: nz, cstep, mpirank, mpisize,step USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER:: Zbounds(0:) INTEGER:: norder(2) INTEGER:: old_Zbounds(0:size(Zbounds,1)-1) INTEGER:: k, i, nbparts REAL(kind=db):: idealnbpartsperproc INTEGER, DIMENSION(0:nz-1):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1) INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction INTEGER:: Zperproc, ierr, remparts CHARACTER(12)::fmt ! calculatese the axial disstibution integrated along the radial direction call calcnbperz(p,partspercol) ! gather this data on all nodes if(step .gt. 0 .and. mpisize .gt. 1) THEN old_Zbounds=Zbounds CALL MPI_ALLREDUCE(MPI_IN_PLACE, partspercol, nz, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr) END IF ! estimate the ideal number of particle per MPI worker idealnbpartsperproc = p%Nptot/mpisize ! find the start and end indices where particles are present Zmin=0 Zmax=nz-1 Do k=0,nz-1 if(partspercol(k) .gt.0) then Zmin=k exit end if end do Do k=nz-1,0,-1 if(partspercol(k) .gt.0) then Zmax=k exit end if end do !Zmin=findloc(partspercol.gt.0,.true.,1,kind=kind(Zmin))-1 !Zmax=findloc(partspercol.gt.0,.true.,1,kind=kind(Zmax),BACK=.true.)-1 ! Find naive axial limits assuming uniform axial distribution IF(Zmax .le. 0) Zmax=nz-1 IF(Zmin .gt. nz) Zmin=0 Zperproc=(Zmax-Zmin)/mpisize IF (Zperproc .lt. 1 .or. cstep .eq. 0) THEN !! No particles are present initially Zperproc=nz/mpisize Zmin=0 ! Define boundaries using naive guess on start or restart (allow to start with 0 parts) DO k=1,mpisize-1 IF(k .lt. mpisize-1-MODULO(Zmax-Zmin,mpisize)) THEN Zbounds(k)=Zmin+k*Zperproc-1 ELSE Zbounds(k)=Zmin+k*Zperproc-1+k-mpisize+2+MODULO(Zmax-Zmin,mpisize) END IF END DO ELSE i=0 ! Define axial boundaries using the axial distribution information. ! the subdomains are not equal remparts=p%Nptot DO k=1,mpisize-1 nbparts=0 DO WHILE(nbparts<0.98*idealnbpartsperproc .and. i .lt. Zmax .and. (nbparts+partspercol(i)).lt.1.25*idealnbpartsperproc) nbparts=nbparts+partspercol(i) i=i+1 END DO remparts=remparts-nbparts Zbounds(k)=i END DO END IF IF(step .gt. 0 .and. mpirank .eq. 0) THEN Do i=1,mpisize-1 !We check that the new limits will not exceed the old limits of the left and right process ! This avoids particle communications with process >mpirank+2 and < mpirank-2 ! However this should converge over time IF(Zbounds(i) .lt. old_Zbounds(i-1)) Zbounds(i)=old_Zbounds(i-1) if(Zbounds(i) .gt. old_Zbounds(i+1))Zbounds(i)=old_Zbounds(i+1) ! If a process would have an axial domain shoter than axial norder, we revert to the old boundaries. IF((Zbounds(i)-Zbounds(i-1)).lt. norder(1) .or. (Zbounds(i+1)-Zbounds(i)).lt. norder(1)) THEN Zbounds=old_Zbounds EXIT END IF END DO END IF ! send the new boundaries to all the workers CALL MPI_Bcast(Zbounds,mpisize+1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) DO k=0,mpisize-1 Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1)) END DO if(mpirank .eq. 0) THEN WRITE(fmt,'(a,i3,a)')"(a,",mpisize+1, "i5)" WRITE(*,fmt) "Zbounds: ", Zbounds WRITE(fmt,'(a,i3,a)')"(a,",mpisize, "i8)" WRITE(*,fmt) "Nplocs: ", Nplocs_all END IF END SUBROUTINE calc_Zbounds !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief After a restart keep only the particles in the local domain of the current MPI worker !--------------------------------------------------------------------------- SUBROUTINE keep_mpi_self_parts(p,Zbounds) TYPE(particles),INTENT(INOUT):: p INTEGER,INTENT(in)::Zbounds(0:) INTEGER :: i, partstart, old_sum,ierr partstart=1 p%Nploc=0 Do i=1,p%Nptot IF(p%Zindex(i).ge.Zbounds(mpirank).and.p%Zindex(i).lt.Zbounds(mpirank+1)) THEN p%Nploc=p%Nploc+1 CALL move_part(p,i,p%Nploc) END IF END DO old_sum=p%Nptot CALL MPI_REDUCE(p%Nploc, p%Nptot,1,MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr) IF(p%Nptot .ne. old_sum) THEN WRITE(*,*) "Error in particle distribution kept: ", p%Nptot, "/",old_sum !call MPI_Abort(MPI_COMM_WORLD, -1, ierr) !stop END IF END SUBROUTINE keep_mpi_self_parts !_______________________________________________________________________________ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Manage the particle communication between neighbours. !> This routine is responsible to receive the incoming particles from the MPI neighbours and to send its outgoing !> particles to these neighbours ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @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, receivednbparts, procs) USE mpihelper, ONLY: particle_type #ifdef _DEBUG USE basic, ONLY: step #endif type(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(out) :: receivednbparts INTEGER, INTENT(in) :: sendholes(:) INTEGER, INTENT(in) :: procs(2) INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0 INTEGER, ASYNCHRONOUS :: sendrequest(2), recvrequest(2) INTEGER, ASYNCHRONOUS :: sendstatus(MPI_STATUS_SIZE,2), recvstatus(MPI_STATUS_SIZE,2) TYPE(particle), ALLOCATABLE :: rrecvpartbuff(:), lrecvpartbuff(:), rsendpartbuff(:), lsendpartbuff(:) ! buffers to send and receive particle from left and right processes INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts, ierr lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL lrecvnbparts=0 rrecvnbparts=0 ! Send and receive the number of particles to exchange CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, procs(1), nbpartsexchange_tag, MPI_COMM_WORLD, recvrequest(1), ierr) CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, procs(2), nbpartsexchange_tag, MPI_COMM_WORLD, recvrequest(2), ierr) CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, procs(1), nbpartsexchange_tag, MPI_COMM_WORLD, sendrequest(1), ierr) CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, procs(2), nbpartsexchange_tag, MPI_COMM_WORLD, sendrequest(2), ierr) CALL MPI_Waitall(2,recvrequest(1:2), recvstatus(:,1:2), ierr) recvrequest=MPI_REQUEST_NULL lreceivednbparts=lrecvnbparts rreceivednbparts=rrecvnbparts ! Re/allocate enough memory to store the incoming particles ALLOCATE(rrecvpartbuff(rreceivednbparts)) ALLOCATE(lrecvpartbuff(lreceivednbparts)) ! Receive particles from left and right processes to the corresponding buffers IF ( lrecvnbparts .gt. 0) THEN CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, procs(1), partsexchange_tag, MPI_COMM_WORLD, recvrequest(1), ierr) END IF IF( rrecvnbparts .gt. 0) THEN CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, procs(2), partsexchange_tag, MPI_COMM_WORLD, recvrequest(2), ierr) END IF ALLOCATE(rsendpartbuff(rsendnbparts)) ALLOCATE(lsendpartbuff(lsendnbparts)) ! Copy the leaving particles to the corresponding send buffers IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes, lsendpartbuff, rsendpartbuff) END IF CALL MPI_Waitall(2,sendrequest(1:2), sendstatus(:,1:2), ierr) ! Send the particles to the left and right neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, procs(1), partsexchange_tag, MPI_COMM_WORLD, sendrequest(1), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", lsendnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, procs(2), partsexchange_tag, MPI_COMM_WORLD, sendrequest(2), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", rsendnbparts , " to right at step: ",step #endif END IF ! Receive the incoming parts in the receive buffers IF ( lreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(1), recvstatus(:,1), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,1) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", lreceivednbparts , " from left at step: ",step #endif END IF IF ( rreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(2), recvstatus(:,2), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,2) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", rreceivednbparts , " from right at step: ",step #endif END IF receivednbparts=rreceivednbparts+lreceivednbparts IF(p%Nploc+receivednbparts-lsendnbparts-rsendnbparts .gt. size(p%R,1)) THEN CALL change_parts_allocation(p,receivednbparts) END IF ! Copy the incoming particles from the receive buffers to the simulation parts variable CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, & & sendholes, lrecvpartbuff, rrecvpartbuff) ! Wait for the outgoing particles to be fully received by the neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(1), sendstatus(:,1), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", lsentnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(2), sendstatus(:,2), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", rsentnbparts , " to right at step: ",step #endif END IF ! ! END SUBROUTINE particlescommunication !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the receive buffers to the local simulation variable parts. !> The incoming particles will first be stored in the holes left by the outgoing particles, then they !> will be added at the end of the parts variable ! !> @param [in] rrecvnbparts number of particles received from the right neighbour (mpirank+1) !> @param [in] lrecvnbparts number of particles received from the left neighbour (mpirank-1) !> @param [in] sendnbparts total number of particles having left the local domain !> @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,lrecvpartbuff, rrecvpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(:) TYPE(particle), INTENT(IN) :: rrecvpartbuff(:), lrecvpartbuff(:) INTEGER k,partpos ! First import the particles coming from the right IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(sendholes(k)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, rrecvpartbuff(k), partpos) END DO END IF ! Then import the particles coming from the left IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(sendholes(k+rrecvnbparts)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, lrecvpartbuff(k), partpos) END DO END IF ! END SUBROUTINE Addincomingparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the local parts variable to the left and right send buffers. ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @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, lsendpartbuff, rsendpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(:) TYPE(particle), INTENT(OUT) :: rsendpartbuff(:), lsendpartbuff(:) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 ! Loop over the outgoing particles and fill the correct send buffer Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(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 SUBROUTINE add_list_created_part(p, buffer,nb_ins) IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:) INTEGER, OPTIONAL:: nb_ins INTEGER:: i, nptotinit, parts_size_increase, nb_added nptotinit=p%Nploc+1 if(present(nb_ins)) THEN nb_added=nb_ins ELSE nb_added=size(buffer,1) end if IF(nb_added .le. 0) RETURN ! No particles to add ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nploc + nb_added .gt. size(p%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 CALL add_created_particle(p,buffer(i)) END DO nb_added=p%Nploc-nptotinit+1 if(p%is_field) then IF(allocated(p%addedlist)) then call change_array_size_int(p%addedlist,2) else allocate(p%addedlist(2)) end if p%addedlist(size(p%addedlist)-1)=nptotinit p%addedlist(size(p%addedlist))=nb_added end if END SUBROUTINE add_list_created_part SUBROUTINE add_linked_created_part(p, linked_buffer, destroy, zerovelocity) IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(linked_part_row), INTENT(in) :: linked_buffer LOGICAL:: destroy, zerovelocity TYPE(linked_part), POINTER:: part INTEGER:: i, nptotinit, parts_size_increase, nb_added nptotinit=p%Nploc+1 nb_added=linked_buffer%n IF(nb_added .le. 0) RETURN ! No particles to add ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nploc + nb_added .gt. size(p%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 part=>linked_buffer%start DO i=1,nb_added CALL add_created_particle(p,part%p) part=>part%next END DO nb_added=p%Nploc-nptotinit+1 if(p%is_field) then IF(allocated(p%addedlist)) then call change_array_size_int(p%addedlist,2) else allocate(p%addedlist(2)) end if p%addedlist(size(p%addedlist)-1)=nptotinit p%addedlist(size(p%addedlist))=nb_added end if if(zerovelocity)then p%UR(nptotinit:p%Nploc)=0 p%UTHET(nptotinit:p%Nploc)=0 p%UZ(nptotinit:p%Nploc)=0 end if if (destroy) call destroy_linked_parts(linked_buffer%start) if (p%is_field) then ! we keep track of energy by removing the ionisation energy ! with conversion from electronvolt to joules loc_etot0=loc_etot0-sum(p%pot(nptotinit:p%Nploc)*elchar) end if END SUBROUTINE add_linked_created_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief 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 add_created_particle(p,part) USE geometry TYPE(particles):: p TYPE(particle):: part p%Nploc=p%Nploc+1 p%newindex=p%newindex+1 ! add the data to the p structure CALL Insertincomingpart(p, part, p%Nploc) p%partindex(p%Nploc)=p%newindex ! calculate the new domain weight CALL dom_weight(p%Z(p%Nploc),p%R(p%Nploc),p%geomweight(p%Nploc,0)) ! delete the particle if it is outside of the computational domain if( .not. is_inside(p,p%Nploc) ) then p%Nploc=p%Nploc-1 p%newindex=p%newindex-1 RETURN end if ! Calculate the geometric weight for the Poisson solver and the grid indices CALL geom_weight(p%Z(p%Nploc),p%R(p%Nploc),p%geomweight(p%Nploc,:)) call p_calc_rzindex(p,p%Nploc) END SUBROUTINE add_created_particle function is_inside(p,id) Use basic, ONLY: rgrid,zgrid, nr, nz IMPLICIT NONE logical :: is_inside type(particles) :: p integer :: id is_inside=.true. if(p%geomweight(id,0).le.0)then is_inside=.false. return end if if(p%R(id).ge.rgrid(nr) .or. p%R(id) .le. rgrid(0))then is_inside=.false. return end if if(p%Z(id).ge.zgrid(nz) .or. p%Z(id) .le. zgrid(0))then is_inside=.false. return end if end function is_inside SUBROUTINE calc_newparts_energy(p) USE basic, ONLY: phinorm, nlclassical type(particles)::p integer::i,n,nptotinit,nbadded, nptotend if(.not. p%is_field) return if( allocated(p%addedlist)) then n=size(p%addedlist) !write(*,*) n, "addedlist: ", p%addedlist Do i=1,n,2 nptotinit=p%addedlist(i) nbadded=p%addedlist(i+1) p%nbadded=p%nbadded+nbadded nptotend=nptotinit+nbadded-1 ! Potential energy loc_etot0=loc_etot0+p%q*p%weight*sum(p%pot(nptotinit:nptotend))*phinorm ! Kinetic energy IF(.not. nlclassical) THEN loc_etot0=loc_etot0+p%m*p%weight*vlight**2*sum(0.5*(p%Gamma(nptotinit:nptotend)+p%Gammaold(nptotinit:nptotend))-1) ELSE loc_etot0=loc_etot0+0.5*p%m*p%weight*vlight**2*sum(p%UR(nptotinit:nptotend)*p%URold(nptotinit:nptotend) & & +p%UZ(nptotinit:nptotend)*p%UZold(nptotinit:nptotend) & & +p%UTHET(nptotinit:nptotend)*p%UTHETold(nptotinit:nptotend)) END IF end do deallocate(p%addedlist) end if end subroutine calc_newparts_energy !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Delete particle at given index removing its energy from the system ! !> @param [in] index index of particle to be deleted !--------------------------------------------------------------------------- SUBROUTINE delete_part(p, index, replace) !! This will destroy particle at the given index USE constants, ONLY: vlight USE bsplines USE geometry USE basic, ONLY: phinorm, nlclassical TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index LOGICAL, OPTIONAL :: replace LOGICAL:: repl IF(present(replace)) THEN repl=replace ELSE repl=.true. END IF !Computes the potential at the particle position with phi_ext+phi_s IF(index .le. p%Nploc) THEN IF(p%is_field) THEN loc_etot0=loc_etot0-p%q*p%weight*(p%pot(index))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0-p%m*p%weight*vlight**2*(p%Gamma(index)-1) ELSE loc_etot0=loc_etot0-0.5*p%m*p%weight*vlight**2*(p%UR(index)**2+p%UZ(index)**2+p%UTHET(index)**2) END IF END IF IF(repl) THEN ! We fill the gap CALL move_part(p, p%Nploc, index) p%partindex(p%Nploc)=-1 ! Reduce the total number of simulated parts p%Nploc=p%Nploc-1 END IF END IF END SUBROUTINE delete_part !_______________________________________________________________________________ SUBROUTINE loaduniformRZ(p, VR,VZ,VTHET) USE basic, ONLY: plasmadim, rnorm, temp, qsim, msim USE constants, ONLY: me, kb, elchar REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=size(VR,1) p%q=sign(elchar,qsim) p%weight=msim/me p%m=me p%qmRatio=qsim/msim ! Initial distribution in z with normalisation CALL loduni(1,p%Z(1:p%Nploc)) p%Z(1:p%Nploc)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*p%Z(1:p%Nploc))/rnorm ! Initial distribution in r with normalisation CALL lodlinr(2,p%R(1:p%Nploc),plasmadim(3),plasmadim(4)) p%R(1:p%Nploc)=p%R(1:p%Nploc)/rnorm ! Initial velocities distribution CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) END SUBROUTINE loaduniformRZ !_______________________________________________________________________________ SUBROUTINE loadDavidson(p, VR,VZ,VTHET, lodr) USE constants, ONLY: me, kb, elchar USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, width, qsim, msim, & & omegac, zgrid, nz, rnorm, n0, nblock, temp procedure(rloader)::lodr TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb, z REAL(kind=db) :: r0, deltar2, halfLz, Mirrorratio, Le, VOL INTEGER :: j, n, blockstart, blockend, addedpart, remainparts INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=p%Nploc Allocate(ra(nblock),rb(nblock), z(0:nblock)) !r0=(plasmadim(4)+plasmadim(3))/2 r0=sqrt(4*H0/(me*omegac**2)) halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) z(0)=plasmadim(1) DO n=1,nblock ! Compute limits in radius and load radii for each part Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1) z(n)=z(0)+n*(plasmadim(2)-plasmadim(1))/nblock deltar2=1-MirrorRatio*cos(2*pi*Le/width) rb(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2)) ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/H0*deltar2)) END DO VOL=SUM(2*pi*MINVAL(ra)*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) qsim=VOL*n0*elchar/nplasma msim=abs(qsim)/elchar*me p%weight=abs(qsim)/elchar p%m=me p%q=sign(elchar,qsim) p%qmRatio=p%q/p%m blockstart=1 blockend=0 ALLOCATE(blocksize(nblock)) WRITE(*,*) "blocksize: ", size(blocksize), nblock DO n=1,nblock blocksize(n)=nplasma/VOL*2*pi*MINVAL(ra)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock END DO remainparts=p%Nploc-SUM(blocksize) addedpart=1 n=nblock/2 j=1 DO WHILE(remainparts .GT. 0) blocksize(n)=blocksize(n)+addedpart remainparts=remainparts-addedpart n=n+j j=-1*(j+SIGN(1,j)) END DO CALL loadPartSlices(p, lodr, ra, rb, z, blocksize) IF(distribtype .eq. 5) THEN CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) VZ=VZ/4 VR=VR*8 VTHET=VTHET*8 ELSE Call loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) END IF END SUBROUTINE loadDavidson SUBROUTINE loadDavidsonVelocities(p, VR,VZ,VTHET, H0, P0) USE constants, ONLY: me, kb, elchar USE basic, ONLY: rnorm, Rcurv, B0, width, vnorm, zgrid, nz TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), INTENT(IN):: H0, P0 REAL(kind=db) :: athetpos, rg, zg, halfLz, Mirrorratio, Pcomp, Acomp INTEGER :: i MirrorRatio=(Rcurv-1)/(Rcurv+1) halfLz=(zgrid(nz)+zgrid(0))/2 ! Load velocities theta velocity ! Loading of r and z velocity is done in adapt_vinit to have ! access to parts%pot DO i=1,p%Nploc ! Interpolation for Magnetic potential rg=p%R(i)*rnorm zg=(p%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/p%m Acomp=-p%qmRatio*Athetpos VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/p%m)),Pcomp+Acomp) !VTHET(i)=Pcomp+Acomp END DO VTHET=VTHET/vnorm VZ=0._db VR=0._db p%Davidson=.true. p%H0=H0 p%P0=P0 END SUBROUTINE loadDavidsonvelocities SUBROUTINE loadGaussianVelocities(p, VR,VZ,VTHET, temperature) USE basic, ONLY: vnorm USE constants, ONLY: kb REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN):: temperature REAL(kind=db):: vth ! Initial velocities distribution vth=sqrt(2.0/3.0*kb*temperature/p%m)/vnorm !thermal velocity CALL lodgaus(3,VZ) CALL lodgaus(5,VR) CALL lodgaus(7,VTHET) VZ=VZ*vth VR=VR*vth VTHET=VTHET*vth p%temperature=temperature p%Davidson=.false. END SUBROUTINE loadGaussianVelocities SUBROUTINE loadFlatTopVelocities(p, VR,VZ,VTHET, meanv, spanv) USE basic, ONLY: vnorm USE constants, ONLY: kb REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT):: meanv(3), spanv(3) ! Initial velocities distribution meanv=meanv/vnorm !thermal velocity spanv=spanv/vnorm CALL loduni(3,VZ) CALL loduni(5,VR) CALL loduni(7,VTHET) VR=(VR*2-1)*spanv(1)+meanv(1) VTHET=(VTHET*2-1)*spanv(2)+meanv(2) VZ=(VZ*2-1)*spanv(3)+meanv(3) p%Davidson=.false. END SUBROUTINE loadFlatTopVelocities SUBROUTINE loadPartslices(p, lodr, ra, rb, z, blocksize) USE basic, ONLY: rnorm TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN)::ra(:), rb(:), z(0:) INTEGER, DIMENSION(:), INTENT(IN) :: blocksize procedure(rloader)::lodr INTEGER :: n, blockstart, blockend, nblock nblock=size(blocksize,1) blockstart=1 blockend=0 DO n=1,nblock blockstart=blockend+1 blockend=MIN(blockstart+blocksize(n)-1,p%Nploc) ! Initial distribution in z with normalisation between magnetic mirrors CALL loduni(1, p%Z(blockstart:blockend)) p%Z(blockstart:blockend)= (z(n-1)+p%Z(blockstart:blockend)*(z(n)-z(n-1)))/rnorm CALL lodr(2, p%R(blockstart:blockend), ra(n), rb(n)) p%R(blockstart:blockend)=p%R(blockstart:blockend)/rnorm END DO END SUBROUTINE loadPartslices !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Read a particle file format to load a simulated specie in the simulation ! !--------------------------------------------------------------------------- SUBROUTINE read_part_file(p, partfilename, VR, VZ, VTHET) USE basic, ONLY: lu_partfile, rnorm, vnorm implicit None TYPE(particles), INTENT(INOUT):: p REAL(kind=db), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VR, VZ, VTHET CHARACTER(len=*)::partfilename INTEGER:: nblock = 0 REAL(kind=db), Dimension(:), ALLOCATABLE:: ra, rb, z INTEGER, Dimension(:), ALLOCATABLE:: npartsslice INTEGER:: velocitytype=1 !< 1) gaussian with temp 2) Davidson with H0, P0 INTEGER:: radialtype=1 !< 1) 1/R 2) uniform 3) 1/R^2 4) gauss INTEGER:: npartsalloc !< initial size of particles arrays REAL(kind=db):: mass=me REAL(kind=db):: charge=-elchar REAL(kind=db):: weight=1.0 REAL(kind=db):: qmratioscale REAL(kind=db):: meanv(3) !< mean velocity in each direction for velocitytype 3 REAL(kind=db):: spanv(3) !< pos/neg extent of velocity in each direction for velocitytype 3 CHARACTER(len=256) :: header=' ' !< header of csv file REAL(kind=db):: H0=3.2e-14 !< Total energy REAL(kind=db):: P0=8.66e-25 !< Canonical angular momentum REAL(kind=db):: temperature=10000 !< temperature in kelvins real(kind=db):: n0 !< density factor LOGICAL :: is_test !< Defines if particle are saved on ittracer or not LOGICAL :: is_field !< Defines if particle contributes to Poisson solver LOGICAL :: calc_moments !< Defines if moments matrix must be calculated each it2d CHARACTER(len=16) :: partformat = 'slice' INTEGER:: i, ierr, openerr NAMELIST /partsload/ nblock, mass, charge, weight, npartsalloc, velocitytype, & & radialtype, temperature, H0, P0, is_test, n0, partformat, meanv, spanv, & & calc_moments, qmratioscale, is_field ! Set defaults qmratioscale=1.0 weight=1.0 meanv=0 spanv=0 mass=me charge=-elchar calc_moments=.false. is_test=.false. is_field=.true. ! Open the paticle file OPEN(UNIT=lu_partfile,FILE=trim(partfilename),ACTION='READ',IOSTAT=openerr) header=' ' IF(openerr .ne. 0) THEN CLOSE(unit=lu_partfile) RETURN END IF READ(lu_partfile,partsload) IF(mpirank .eq.0) THEN WRITE(*,'(a,a)')"reading partfile: ", trim(partfilename) WRITE(*,partsload) END IF ! The plasma cloud is defined as a set of slices IF(trim(partformat).eq.'slice') THEN IF( nblock .ge. 1) THEN ALLOCATE(z(0:nblock),ra(nblock),rb(nblock), npartsslice(nblock)) DO WHILE(header(1:8) .ne. '//slices') READ(lu_partfile,'(a)') header END DO DO i=1,nblock READ(lu_partfile,*) z(i-1),ra(i),rb(i),npartsslice(i) END DO READ(lu_partfile,*) z(nblock) CALL creat_parts(p,max(npartsalloc,sum(npartsslice))) p%Nploc=sum(npartsslice) p%Nptot=p%Nploc IF( allocated(VR) ) THEN DEALLOCATE(VR,VZ,VTHET) end if if(.not. allocated(VR)) THEN ALLOCATE(VR(p%Nploc)) ALLOCATE(VZ(p%Nploc)) ALLOCATE(VTHET(p%Nploc)) END IF p%m=mass p%q=charge p%weight=weight p%qmRatio=charge/mass*qmratioscale p%is_test=is_test p%is_field=is_field p%calc_moments=calc_moments p%Newindex=sum(npartsslice) SELECT CASE(radialtype) CASE(1) ! 1/R distribution in R CALL loadPartslices(p, lodunir, ra, rb, z, npartsslice) CASE(2) ! flat top distribution in R CALL loadPartslices(p, lodlinr, ra, rb, z, npartsslice) CASE(3) ! 1/R^2 distribution in R CALL loadPartslices(p, lodinvr, ra, rb, z, npartsslice) CASE(4) ! gaussian distribution in R CALL loadPartslices(p, lodgausr, ra, rb, z, npartsslice) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT SELECT CASE(velocitytype) CASE(1) ! Gaussian with temperature CALL loadGaussianVelocities(p, VR, VZ, VTHET, temperature) CASE(2) ! Davidson magnetic mirror high wr equilibrium CALL loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) CASE(3) ! flat top velocity CALL loadFlatTopVelocities(p, VR, VZ, VTHET, meanv, spanv) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of velocity distribution:", velocitytype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT END IF END IF ! The plasma cloud is defined as a set individual particles IF( trim(partformat) .eq. 'parts' ) THEN IF( nblock .ge. 1) THEN !Allocate necessary memory CALL creat_parts(p,max(npartsalloc,nblock)) IF( allocated(VR) ) THEN DEALLOCATE(VR,VZ,VTHET) end if if(.not. allocated(VR)) THEN ALLOCATE(VR(nblock)) ALLOCATE(VZ(nblock)) ALLOCATE(VTHET(nblock)) END IF ! Read the particles from the file DO WHILE(header(1:8) .ne. '//parts') READ(lu_partfile,'(a)') header END DO DO i=1,nblock READ(lu_partfile,*) p%R(i),p%THET(i),p%Z(i), VR(i), VTHET(i), VZ(i) END DO p%Nploc=nblock p%Nptot=p%Nploc p%m=mass p%q=charge p%Newindex=nblock p%weight=weight p%qmRatio=charge/mass*qmratioscale p%is_test=is_test p%is_field=is_field p%calc_moments=calc_moments !normalizations p%r=p%r/rnorm p%z=p%z/rnorm VR=VR/vnorm VTHET=VTHET/vnorm VZ=VZ/vnorm END IF END IF CLOSE(unit=lu_partfile) END SUBROUTINE !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Increase the number of macroparticles by separating each previous macroparticles into !> samplefactor new macroparticles of equally divided weight. The new sub particles are distributed !> uniformly in space to maintain the density and other moments. ! !> @param [in] samplefactor multiplicator of the number of macroparticles. !> @param [in] p particles type to increase. !--------------------------------------------------------------------------- SUBROUTINE upsample(p, samplefactor) USE basic, ONLY : nplasma, dr, dz INTEGER, INTENT(IN) ::samplefactor TYPE(particles), INTENT(INOUT):: p INTEGER:: i, j, currentindex REAL(kind=db), DIMENSION(p%Nploc) :: spreaddir ! random direction for the spread of each initial macro particle REAL(kind=db) :: dir ! Direction in which the particle is moved REAL(kind=db) :: dl ! Particle displacement used for ! Load and scale the direction angle for spreading the new particles CALL loduni(2, spreaddir) spreaddir=spreaddir*2*pi/samplefactor dl=min(dz,minval(dr,1,dr.GT.0))/100 DO i=1,p%Nploc DO j=1,samplefactor-1 currentindex=p%Nploc+(i-1)*(samplefactor-1)+j CALL move_part(p,i,currentindex) p%partindex(currentindex)=currentindex dir = spreaddir(i)+2*pi*j/samplefactor p%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 p%weight=p%weight/samplefactor p%Nploc=p%Nploc*samplefactor p%Nptot=p%Nptot*samplefactor END SUBROUTINE upsample ! Taken from https://rosettacode.org/wiki/Sorting_algorithms/Radix_sort#Fortran ! No Copyright is exerted due to considerable prior art in the Public Domain. ! This Fortran version by Peter Kelly ~ peter.kelly@acm.org ! ! Permission is hereby granted, free of charge, to any person obtaining ! a copy of this software and associated documentation files (the ! "Software"), to deal in the Software without restriction, including ! without limitation the rights to use, copy, modify, merge, publish, ! distribute, sublicense, and/or sell copies of the Software, and to ! permit persons to whom the Software is furnished to do so, subject to ! the following conditions: ! The above copyright notice and this permission notice shall be ! included in all copies or substantial portions of the Software. ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. ! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY ! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, ! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE ! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ! ! Implementation of a classic Radix Sort LSD style :) SUBROUTINE LSDRADIXSORT(A , N) IMPLICIT NONE ! ! Dummy arguments ! INTEGER :: N INTEGER , target, DIMENSION(0:N - 1) :: A ! All arrays based off zero, one day I'll fix it INTENT (IN) N INTENT (INOUT) A ! ! Local variables ! INTEGER , DIMENSION(0:9) :: counts INTEGER :: digitplace INTEGER :: i INTEGER :: j INTEGER :: largestnum INTEGER, DIMENSION(0:N - 1) :: results ! digitplace = 1 ! Count of the keys largestnum = MAXVAL(A) DO WHILE ( (largestnum/digitplace)>0 ) counts = 0 ! Init the count array DO i = 0 , N - 1 , 1 J = (A(i)/digitplace) J = MODULO(j , 10) counts(j) = counts(j) + 1 END DO ! Change count(i) so that count(i) now contains actual position of this digit in result() ! Working similar to the counting sort algorithm DO i = 1 , 9 , 1 counts(i) = counts(i) + counts(i - 1) ! Build up the prefix sum END DO ! DO i = N - 1 , 0 , -1 ! Move from left to right j = (A(i)/digitplace) j = MODULO(j, 10) results(counts(j) - 1) = A(i) ! Need to subtract one as we are zero based but prefix sum is 1 based counts(j) = counts(j) - 1 END DO ! DO i = 0 , N - 1 , 1 ! Copy the semi-sorted data into the input A(i) = results(i) END DO ! digitplace = digitplace*10 END DO ! While loop RETURN END SUBROUTINE LSDRADIXSORT END MODULE beam diff --git a/src/diagnose.f90 b/src/diagnose.f90 index f12d1d7..8fa4880 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,434 +1,467 @@ SUBROUTINE diagnose(kstep) ! ! Diagnostics ! USE basic USE futils USE hashtable Use maxwsrce Use neutcol USE beam, ONLY : partslist, epot, ekin, etot, etot0, Nplocs_all, collectparts #if USE_X == 1 USE xg, ONLY : initw, updt_xg_var #endif USE fields, ONLY: phi_spline, nbmoments USE celldiag use mpi Use geometry Use splinebound Use weighttypes use psupply use filemanip IMPLICIT NONE ! INTEGER, INTENT(in) :: kstep ! ! Local vars and arrays INTEGER, PARAMETER :: BUFSIZE = 20 CHARACTER(len=128) :: str, fname, grpname CHARACTER(len=12):: charspec INTEGER:: Ntotal ! Total number of simulated particles remaining in the simulation INTEGER:: partsrank, partsdim(2) INTEGER:: Nplocs_all_save(64) INTEGER:: i, nbbounds INTEGER, allocatable, save:: partnbBuffer(:,:) + REAL(kind=db), ALLOCATABLE :: magr(:), magz(:) + REAL(kind=db), ALLOCATABLE :: tempBr(:, :), tempBz(:, :), tempAthet(:, :) + INTEGER :: magn(2), magrank, magfid !________________________________________________________________________________ ! 1. Initial diagnostics IF( kstep .EQ. 0 .and. mpirank .eq. 0) THEN ! Only process 0 should save on file ! WRITE(*,'(a)') ' Initial diagnostics' ! ! 1.1 Initial run or when NEWRES set to .TRUE. ! IF( .NOT. nlres .OR. newres) THEN CALL creatf(resfile, fidres, 'Espic2d result file', real_prec='d') WRITE(*,'(3x,a,a)') TRIM(resfile), ' created' ! ! Label the run IF( LEN_TRIM(label1).GT.0 ) CALL attach(fidres, "/", "label1", TRIM(label1)) IF( LEN_TRIM(label2).GT.0 ) CALL attach(fidres, "/", "label2", TRIM(label2)) IF( LEN_TRIM(label3).GT.0 ) CALL attach(fidres, "/", "label3", TRIM(label3)) IF( LEN_TRIM(label4).GT.0 ) CALL attach(fidres, "/", "label4", TRIM(label4)) ! ! Job number jobnum = 0 ! ! Data group CALL creatg(fidres, "/data", "data") CALL creatg(fidres, "/data/var0d", "0d history arrays") CALL creatg(fidres, "/data/var1d","1d history arrays") CALL creatg(fidres, "/data/part", "part phase space") CALL creatg(fidres, "/data/fields", "Electro static potential and Er Ez fields") ! ! File group CALL creatg(fidres, "/files", "files") CALL attach(fidres, "/files", "jobnum", jobnum) CALL putarr(fidres, "/data/var1d/rgrid", rgrid) CALL putarr(fidres, "/data/var1d/zgrid", zgrid) CALL creatd(fidres, 1, (/ 64 /), "/data/var0d/Nplocs_all", "Nplocs_all") CALL creatd(fidres, 1, (/3/), "/data/var0d/nudcol", "nudcol") ! Part and fields vectors ! Initialize time-dependent particle variables CALL creatd(fidres, 1, SHAPE(partslist(1)%R), "/data/part/R", "R",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/Z", "Z",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/THET", "THET",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%UZ), "/data/part/UR", "UZ",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%UR), "/data/part/UZ", "UR",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%UTHET), "/data/part/UTHET", "UTHET",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Rindex), "/data/part/Rindex", "Rindex",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Zindex), "/data/part/Zindex", "Zindex",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%partindex), "/data/part/partindex", "partindex",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%pot), "/data/part/pot", "pot",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 0, SHAPE(time), "/data/part/time", "time" ) CALL creatd(fidres, 0, SHAPE(Ntotal), "/data/part/Nparts", "number of remaining parts in the simulation space") CALL creatd(fidres, 1, (/7/), "/data/part/nbchange", "number of added parts, lost parts zm,zp,rm,rp, and collisions per type io, ela") CALL attach(fidres,'/data/part', "q", partslist(1)%q) CALL attach(fidres,'/data/part', "m", partslist(1)%m) CALL attach(fidres,'/data/part', "weight", partslist(1)%weight) CALL creatd(fidres, 1, SHAPE(pot), "/data/fields/pot", "pot") CALL creatd(fidres, 1, SHAPE(Er), "/data/fields/Er", "Er") CALL creatd(fidres, 1, SHAPE(Ez), "/data/fields/Ez", "Ez") CALL creatd(fidres, 1, SHAPE(phi_spline), "/data/fields/phi", "spline form of Phi") CALL creatd(fidres, 2, (/nbmoments,nrank(1)*nrank(2)/), "/data/fields/moments", "moments",compress=.true.,chunking=(/1,nrank(1)*nrank(2)/)) !CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments") CALL creatd(fidres, 0, SHAPE(time), "/data/fields/time", "time" ) CALL putarr(fidres, "/data/fields/Br", Br) CALL putarr(fidres, "/data/fields/Bz", Bz) CALL putarr(fidres, "/data/fields/Athet", Athet) CALL putarr(fidres, "/data/fields/volume", Volume) + + ! We save the magnetic field as defined by the h5 file + if(len_trim(magnetfile) .gt. 1) then + CALL openf(trim(magnetfile), magfid, 'r', real_prec='d') + + CALL getdims(magfid, '/mag/Athet', magrank, magn) + + ALLOCATE (magr(magn(2)), magz(magn(1))) + ALLOCATE (tempAthet(magn(1), magn(2)), tempBr(magn(1), magn(2)), tempBz(magn(1), magn(2))) + + ! Read r and z coordinates for the definition of A_\thet, and B + CALL getarr(magfid, '/mag/r', magr) + CALL getarr(magfid, '/mag/z', magz) + CALL getarr(magfid, '/mag/Athet', tempAthet) + + IF (isdataset(magfid, '/mag/Br') .and. isdataset(magfid, '/mag/Bz')) THEN + CALL getarr(magfid, '/mag/Br', tempBr) + CALL getarr(magfid, '/mag/Bz', tempBz) + end if + + CALL creatg(fidres, '/data/inputmag') + CALL putarr(fidres, '/data/inputmag/r',magr) + CALL putarr(fidres, '/data/inputmag/z',magz) + CALL putarr(fidres, '/data/inputmag/Athet',tempAthet) + CALL putarr(fidres, '/data/inputmag/Br',tempBr) + CALL putarr(fidres, '/data/inputmag/Bz',tempBz) + + + call closef(magfid) + end if ! ! 1.2 Restart run ! ELSE CALL cp2bk(resfile) ! backup previous result file CALL openf(resfile, fidres, real_prec='d') WRITE(*,'(3x,a,a)') TRIM(resfile), ' open' CALL getatt(fidres, "/files", "jobnum", jobnum) jobnum = jobnum+1 WRITE(*,'(3x,a,i3)') "Current Job Number =", jobnum CALL attach(fidres, "/files", "jobnum", jobnum) !allocate(partnbBuffer(nbspecies,4+size(partslist(1)%nblost,1))) !partnbBuffer=0 END IF ! ! Add input namelist variables as attributes of /data/input WRITE(str,'(a,i2.2)') "/data/input.",jobnum CALL creatg(fidres, TRIM(str)) CALL attach(fidres, TRIM(str), "job_time", job_time) CALL attach(fidres, TRIM(str), "extra_time", extra_time) CALL attach(fidres, TRIM(str), "dt", dt*tnorm) CALL attach(fidres, TRIM(str), "tmax", tmax) CALL attach(fidres, TRIM(str), "nrun", nrun) CALL attach(fidres, TRIM(str), "nlres", nlres) CALL attach(fidres, TRIM(str), "nlsave", nlsave) CALL attach(fidres, TRIM(str), "newres", newres) CALL attach(fidres, TRIM(str), "nz", nz) CALL putarr(fidres, TRIM(str)//"/lz", lz) CALL attach(fidres, TRIM(str), "nplasma", nplasma) CALL attach(fidres, TRIM(str), "potinn", potinn) CALL attach(fidres, TRIM(str), "potout", potout) CALL attach(fidres, TRIM(str), "B0", B0) CALL attach(fidres, TRIM(str), "Rcurv", Rcurv) CALL attach(fidres, TRIM(str), "width", width) CALL attach(fidres, TRIM(str), "n0", n0) CALL attach(fidres, TRIM(str), "temp", partslist(1)%temperature) CALL attach(fidres, TRIM(str), "it0d", it0d) CALL attach(fidres, TRIM(str), "it2d", it2d) CALL attach(fidres, TRIM(str), "itparts", itparts) CALL attach(fidres, TRIM(str), "nlclassical", nlclassical) CALL attach(fidres, TRIM(str), "nlPhis", nlPhis) CALL attach(fidres, TRIM(str), "qsim", partslist(1)%q*partslist(1)%weight) CALL attach(fidres, TRIM(str), "msim", partslist(1)%m*partslist(1)%weight) CALL attach(fidres, TRIM(str), "startstep", cstep) CALL attach(fidres, TRIM(str), "H0", partslist(1)%H0) CALL attach(fidres, TRIM(str), "P0", partslist(1)%P0) CALL putarr(fidres, TRIM(str)//"/femorder", femorder) CALL putarr(fidres, TRIM(str)//"/ngauss", ngauss) CALL putarr(fidres, TRIM(str)//"/nnr", nnr) CALL putarr(fidres, TRIM(str)//"/radii", radii) CALL putarr(fidres, TRIM(str)//"/plasmadim", plasmadim) CALL attach(fidres, TRIM(str), "rawparts", .true.) CALL attach(fidres, TRIM(str), "nbspecies", nbspecies) CALL putarr(fidres, TRIM(str)//"/potxt", potxt) CALL putarr(fidres, TRIM(str)//"/Erxt", Erxt) CALL putarr(fidres, TRIM(str)//"/Ezxt", Ezxt) ! Save geometry parameters for non conforming boundary conditions Call geom_diag(fidres,str,rnorm) ! Save geometry parameters for non conforming boundary conditions using b-spline curves call splinebound_diag(fidres, str, the_domain) ! Save Maxwellsource parameters for the ad-hoc source Call maxwsrce_diag(fidres,str,vnorm) ! Save neutcol parameters for the electron collisions with neutrals Call neutcol_diag(fidres,str,vnorm) if(.not. isdataset(fidres,'/data/var0d/nudcol'))then CALL creatd(fidres, 1, (/3/), "/data/var0d/nudcol", "nudcol") end if ! Save psupply parameters for the simulation of realistic power supplies Call psupply_diag(fidres,str) if(.not. isdataset(fidres,'/data/var0d/biases'))then nbbounds=2 if(the_domain%nbsplines .gt. 0) nbbounds=the_domain%nbsplines CALL creatd(fidres, 1, (/nbbounds/), "/data/var0d/biases", "biases") end if ! Save STDIN of this run WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum INQUIRE(unit=lu_in, name=fname) CALL putfile(fidres, TRIM(str), TRIM(fname)) ! Prepare hdf5 file for storing test particles DO i=2,nbspecies WRITE(grpname,'(a,i2)')'/data/part/',i call create_parts_group(partslist(i),trim(grpname),time) END DO CALL attach(fidres, "/data/part", "nbspecies", nbspecies) ! ! Initialize buffers for 0d history arrays CALL htable_init(hbuf0, BUFSIZE) CALL set_htable_fileid(hbuf0, fidres, "/data/var0d") ! ! Initialize Xgrafix #if USE_X == 1 IF(nlxg) THEN CALL initw END IF #endif END IF IF(kstep .EQ. 0) THEN ! Initialize particle cell diagnostic CALL celldiag_init(lu_in, fidres) CLOSE(lu_in) allocate(partnbBuffer(nbspecies,4+size(partslist(1)%nblost,1))) partnbBuffer=0 END IF !________________________________________________________________________________ ! 2. Periodic diagnostics IF( kstep .NE. -1) THEN IF(modulo(step,ittracer) .eq. 0 .or. nlend) THEN ! We gather the traced particles on the mpi host DO i=1,nbspecies IF(partslist(i)%is_test) CALL collectparts(partslist(i)) END DO END IF IF(modulo(step,itrestart) .eq. 0 .or. modulo(step,itparts) .eq. 0 .or. nlend) THEN ! We gather the traced particles on the mpi host DO i=1,nbspecies CALL collectparts(partslist(i)) END DO END IF do i=1,nbspecies partnbBuffer(i,1)=partnbBuffer(i,1)+partslist(i)%nbadded partnbBuffer(i,2:3)=partnbBuffer(i,2:3)+partslist(i)%nbcolls partnbBuffer(i,4)=partslist(i)%Nploc partnbBuffer(i,5:)=partnbBuffer(i,5:)+partslist(i)%nblost partslist(i)%nbadded=0 partslist(i)%nblost=0 partslist(i)%nbcolls=0 end do IF(modulo(step,ittext) .eq. 0 .or. nlend) THEN ! We gather the number of gained and lost particles on the mpi host IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, partnbBuffer, nbspecies*(4+size(partslist(1)%nblost,1)), MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_REDUCE(partnbBuffer, partnbBuffer, nbspecies*(4+size(partslist(1)%nblost,1)), MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) partnbBuffer=0 END IF end if ! ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN ! IF (mpisize .gt. 1) THEN partslist(1)%Nptot=sum(Nplocs_all) END IF ! IF(modulo(step,ittext).eq. 0 .or. nlend) THEN WRITE(*,'(a,1x,i8.8,a1,i8.8,20x,a,1pe10.3)') '*** Timestep (this run/total) =', & & step, '/', cstep, 'Time =', time if( abs(etot).gt. 0) then WRITE(*,'(a,6(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Etot0, Eerr, Eerr rel, Nbparts/Ntotal', epot, ekin, etot, etot0, etot-etot0,(etot-etot0)/etot, partslist(1)%Nptot,'/',nplasma else WRITE(*,'(a,4(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Eerr, Nbparts/Ntotal', epot, ekin, etot, etot-etot0, partslist(1)%Nptot,'/',nplasma end if IF(mpisize .gt. 1 ) then WRITE(*,'(a,64i10.7)') 'Nbparts per proc', Nplocs_all end if Write(*,'(a)')"speci, added, iocoll, elacoll, tot var, tot, Losses (zmin zmax rmin rmax boundaries(i))" write(charspec,'(a,i02,a)') '(i04,',size(partnbBuffer,2)+1,'i9.7)' do i=1,nbspecies WRITE(*,charspec) i, partnbBuffer(i,1),partnbBuffer(i,2:3), partnbBuffer(i,1)-sum(partnbBuffer(i,5:)), partnbBuffer(i,4),-partnbBuffer(i,5:) partslist(i)%nptot= partnbBuffer(i,4) end do partnbBuffer=0 END IF !________________________________________________________________________________ ! ! 2.1 0d history arrays ! ! if we do a restart, we don't want to save the same data twice IF( kstep .eq. 0 .and. nlres .and. (.not. newres)) return IF(modulo(step,it0d).eq. 0 .or. nlend) THEN CALL add_record(hbuf0, "time", "simulation time", time) CALL add_record(hbuf0, "epot", "potential energy", epot) CALL add_record(hbuf0, "ekin", "kinetic energy", ekin) CALL add_record(hbuf0, "etot", "total energy", etot) CALL add_record(hbuf0, "etot0", "theoretical total energy", etot0) CALL add_record(hbuf0, "nbparts", "number of remaining parts in the simulation space", REAL(partslist(1)%Nptot,kind=db)) CALL add_record(hbuf0,"current", "unscaled current flowing between the electrodes of the power supplies [A]", the_ps%current(1)*qnorm/tnorm) CALL htable_endstep(hbuf0) Nplocs_all_save=0 Nplocs_all_save(1:mpisize)=Nplocs_all(0:mpisize-1) CALL append(fidres, "/data/var0d/Nplocs_all", REAL(Nplocs_all_save,kind=db)) CALL append(fidres, "/data/var0d/nudcol", partslist(1)%nudcol/(dt*tnorm)) CALL append(fidres, "/data/var0d/biases", the_ps%biases*phinorm) END IF ! ! 2.2 2d profiles IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL append(fidres, "/data/fields/time", time) CALL append(fidres, "/data/fields/pot", pot*phinorm) CALL append(fidres, "/data/fields/Er", Er*enorm) CALL append(fidres, "/data/fields/Ez", Ez*enorm) CALL append(fidres, "/data/fields/phi", phi_spline*phinorm) CALL append(fidres, "/data/fields/moments", partslist(1)%moments) DO i=2,nbspecies IF ( .not. partslist(i)%calc_moments) CYCLE WRITE(grpname,'(a,i2,a)')'/data/part/',i,'/' CALL append(fidres, trim(grpname) // "moments", partslist(i)%moments) end DO END IF ! ! 2.3 main specie quantities IF(modulo(step,itparts).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' CALL append(fidres, "/data/part/time", time) CALL append(fidres, "/data/part/Nparts", REAL(partslist(1)%Nptot,kind=db)) !CALL append(fidres, "/data/part/nbchange", REAL((/partslist(1)%nbadded,partslist(1)%nblost,partslist(1)%nbcolls/),kind=db)) IF ( isdataset(fidres,'/data/part/R') ) THEN CALL getdims(fidres, '/data/part/R', partsrank, partsdim) partsdim(1)=min(size(partslist(1)%R,1), partsdim(1)) CALL append(fidres, "/data/part/R", partslist(1)%R(1:partsdim(1))*rnorm) CALL append(fidres, "/data/part/Z", partslist(1)%Z(1:partsdim(1))*rnorm) CALL append(fidres, "/data/part/THET", partslist(1)%THET(1:partsdim(1))) CALL append(fidres, "/data/part/UZ", 0.5*(partslist(1)%UZ(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))+partslist(1)%UZold(1:partsdim(1))/partslist(1)%gammaold(1:partsdim(1)))) CALL append(fidres, "/data/part/UR", 0.5*(partslist(1)%UR(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))+partslist(1)%URold(1:partsdim(1))/partslist(1)%gammaold(1:partsdim(1)))) CALL append(fidres, "/data/part/UTHET", 0.5*(partslist(1)%UTHET(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))+partslist(1)%UTHETold(1:partsdim(1))/partslist(1)%gammaold(1:partsdim(1)))) CALL append(fidres, "/data/part/pot", partslist(1)%pot(1:partsdim(1))*phinorm) CALL append(fidres, "/data/part/Rindex", REAL(partslist(1)%Rindex(1:partsdim(1)),kind=db)) CALL append(fidres, "/data/part/Zindex", REAL(partslist(1)%Zindex(1:partsdim(1)),kind=db)) CALL append(fidres, "/data/part/partindex", REAL(partslist(1)%partindex(1:partsdim(1)),kind=db)) END IF END IF ! ! 2.4 Tracer quantities IF(modulo(step,ittracer).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' DO i=2,nbspecies IF ( .not. partslist(i)%is_test) CYCLE WRITE(grpname,'(a,i2,a)')'/data/part/',i,'/' CALL append(fidres, trim(grpname) // "time", time) CALL append(fidres, trim(grpname) //"Nparts", REAL(partslist(i)%Nptot,kind=db)) !CALL append(fidres, trim(grpname) //"nbchange", REAL((/partslist(i)%nbadded,partslist(i)%nblost,partslist(i)%nbcolls/),kind=db)) IF ( isdataset(fidres,trim(grpname)//'R') ) THEN CALL getdims(fidres, trim(grpname) // 'R', partsrank, partsdim) partsdim(1)=min(size(partslist(i)%R,1), partsdim(1)) CALL append(fidres, trim(grpname) // "R", partslist(i)%R(1:partsdim(1))*rnorm) CALL append(fidres, trim(grpname) // "Z", partslist(i)%Z(1:partsdim(1))*rnorm) CALL append(fidres, trim(grpname) // "THET", partslist(i)%THET(1:partsdim(1))) CALL append(fidres, trim(grpname) // "UZ", 0.5*(partslist(i)%UZ(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)) + partslist(i)%UZold(1:partsdim(1))/partslist(i)%gammaold(1:partsdim(1)))) CALL append(fidres, trim(grpname) // "UR", 0.5*(partslist(i)%UR(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)) + partslist(i)%URold(1:partsdim(1))/partslist(i)%gammaold(1:partsdim(1)))) CALL append(fidres, trim(grpname) // "UTHET", 0.5*(partslist(i)%UTHET(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)) + partslist(i)%UTHETold(1:partsdim(1))/partslist(i)%gammaold(1:partsdim(1)))) CALL append(fidres, trim(grpname) // "pot", partslist(i)%pot(1:partsdim(1))*phinorm) CALL append(fidres, trim(grpname) // "Rindex", REAL(partslist(i)%Rindex(1:partsdim(1)),kind=db)) CALL append(fidres, trim(grpname) // "Zindex", REAL(partslist(i)%Zindex(1:partsdim(1)),kind=db)) CALL append(fidres, trim(grpname) // "partindex", REAL(partslist(i)%partindex(1:partsdim(1)),kind=db)) END IF END DO ! END IF ! 2.5 3d profiles ! ! ! 2.6 Xgrafix ! #if USE_X == 1 IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0) THEN call xgevent CALL updt_xg_var CALL xgupdate END IF #endif !________________________________________________________________________________ ! 3. Final diagnostics ELSE ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN ! ! Flush 0d history array buffers CALL htable_hdf5_flush(hbuf0) ! ! Close all diagnostic files CALL closef(fidres) !________________________________________________________________________________ END IF ! CONTAINS SUBROUTINE create_parts_group(p,grpname, time) USE beam,ONLY: particles type(particles):: p real(kind=db):: time character(len=*):: grpname If(isgroup(fidres, trim(grpname))) return CALL creatg(fidres, grpname, "specific specie phase space") CALL creatd(fidres, 0, SHAPE(time), trim(grpname) // "/time", "time") CALL creatd(fidres, 0, SHAPE(time), trim(grpname) //"/Nparts", "number of remaining parts") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/R", "radial pos") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/Z", "axial pos") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/THET", "azimuthal pos") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/UZ", "axial beta*gamma") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/UR", "radial beta*gamma") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/UTHET", "azimuthal beta*gamma") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/pot", "electric potential") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/Rindex", "radial grid index") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/Zindex", "axial grid index") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/partindex", "particle index") CALL creatd(fidres, 1, (/7/), trim(grpname) // "nbchange", "number of added parts, lost parts zm,zp,rm,rp, and collisions per type io, ela") CALL attach(fidres,trim(grpname), "q", p%q) CALL attach(fidres,trim(grpname), "m", p%m) CALL attach(fidres,trim(grpname), "weight", p%weight) CALL creatd(fidres, 2, (/nbmoments,nrank(1)*nrank(2)/), trim(grpname) // "/moments", "moments",compress=.true.,chunking=(/1,nrank(1)*nrank(2)/)) END SUBROUTINE create_parts_group END SUBROUTINE diagnose diff --git a/src/neutcol_mod.f90 b/src/neutcol_mod.f90 index a686c86..b848f37 100644 --- a/src/neutcol_mod.f90 +++ b/src/neutcol_mod.f90 @@ -1,454 +1,454 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: neutcol ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for handling the electron-neutral collisions and creating electrons !> by ionisation Based on the paper by Birdsall 1991 and Sengupta et al. !------------------------------------------------------------------------------ module neutcol USE constants IMPLICIT NONE private LOGICAL :: nlcol=.false. !< Flag to activate or not electron neutral collisions LOGICAL :: nlmaxwellio=.false. !< Flag to define how ionised electrons are created (physically or according to maxwellian) INTEGER :: itcol = 1 !< number of dt between each evaluation of neutcol_step Real(kind=db) :: neutdens=2.4e16 !< Neutral particle density in m-3 Real(kind=db) :: neuttemp=300 !< Neutral particle temperature in K Real(kind=db) :: neutpressure=1e-6 !< Neutral particle pressure in mbar Real(kind=db) :: scatter_fac = 24.2 !< Energy scattering factor for the considered gas (here for Ne) real(kind=db) :: Eion = 21.56 !< Ionisation energy (eV) (here for Ne) Real(kind=db) :: E0 = 27.21 !< Atomic unit of energy used for calculation of deviation angles real(kind=db) :: collfactor !< collision factor (n_n \sigma \delta t) INTEGER :: nb_io_cross=0 Real(kind=db), ALLOCATABLE :: io_cross_sec(:,:) !< Ionisation cross-section table Real(kind=db), ALLOCATABLE :: io_growth_cross_sec(:) !< Ionisation exponential fitting factor INTEGER :: nb_ela_cross=0 Real(kind=db), ALLOCATABLE :: ela_cross_sec(:,:) !< Elastic collision cross section table Real(kind=db), ALLOCATABLE :: ela_growth_cross_sec(:) !< Elastic collision exponential fitting factor Real(kind=db) :: Escale !< Energy normalisation factor used to reduce computation costs CHARACTER(len=128) :: io_cross_sec_file='' CHARACTER(len=128) :: ela_cross_sec_file='' Real(kind=db) :: etemp=22000 !< In case of nlmaxwelio, defines the temperature of created electrons Real(kind=db) :: vth !< In case of nlmaxwelio, defines the thermal velocity of created electrons LOGICAL :: nldragio=.true. !< INTEGER :: species(2) NAMELIST /neutcolparams/ neutdens, neuttemp, neutpressure, Eion, & & scatter_fac, nlcol, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, & & nldragio, itcol, species PUBLIC:: neutcol_init, neutcol_step, neutcol_diag, itcol, neutdens CONTAINS subroutine neutcol_init(lu_in, p) use mpi Use basic, only: mpirank, dt, nlclassical,rnorm,tnorm, vnorm Use beam, only: particles Use constants implicit none INTEGER, INTENT(IN) :: lu_in TYPE(particles) :: p INTEGER:: ierr, istat character(len=1000) :: line species(1)=1 species(2)=-1 Rewind(lu_in) READ(lu_in, neutcolparams, iostat=istat) if (istat.gt.0) then backspace(lu_in) read(lu_in,fmt='(A)') line write(*,'(A)') & 'Invalid line in neutcolparams: '//trim(line) call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if IF(mpirank .eq. 0) THEN WRITE(*, neutcolparams) END IF if(.not. nlcol) return if(nlclassical)THEN Escale=0.5*p%m/elchar*vlight**2 else Escale=p%m*vlight**2/elchar end if if (nlmaxwellio) vth=sqrt(kb*etemp/p%m)/vnorm if(io_cross_sec_file .ne.'') then call read_cross_sec(io_cross_sec_file,io_cross_sec, nb_io_cross) if(nb_io_cross .gt. 0) then allocate(io_growth_cross_sec(nb_io_cross-1)) ! Normalisations io_cross_sec(:,2)=io_cross_sec(:,2)/rnorm**2 ! Precomputing of exponential fitting factor for faster execution io_growth_cross_sec=log(io_cross_sec(2:nb_io_cross,2)/io_cross_sec(1:nb_io_cross-1,2))/ & & log(io_cross_sec(2:nb_io_cross,1)/io_cross_sec(1:nb_io_cross-1,1)) end if end if if(ela_cross_sec_file .ne.'') then call read_cross_sec(ela_cross_sec_file,ela_cross_sec, nb_ela_cross) if(nb_ela_cross .gt. 0) then allocate(ela_growth_cross_sec(nb_ela_cross-1)) ! Normalisations ela_cross_sec(:,2)=ela_cross_sec(:,2)/rnorm**2 ! Precomputing of exponential fitting factor for faster execution ela_growth_cross_sec=log(ela_cross_sec(2:nb_ela_cross,2)/ela_cross_sec(1:nb_ela_cross-1,2))/ & & log(ela_cross_sec(2:nb_ela_cross,1)/ela_cross_sec(1:nb_ela_cross-1,1)) end if END IF nlcol=nlcol .and. (allocated(io_cross_sec) .or. allocated(ela_cross_sec)) ! Collision factor depending on neutral gas parameters collfactor=neutdens*dt*rnorm**3*itcol neutpressure=neutdens*kb*300/100 end subroutine neutcol_init Subroutine neutcol_diag(File_handle, str, vnorm) use mpi Use futils Integer:: File_handle Real(kind=db):: vnorm Character(len=*):: str CHARACTER(len=256):: grpname Integer:: ierr, mpirank CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0 .and. nlcol) THEN Write(grpname,'(a,a)') trim(str),"/neutcol" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "neutdens", neutdens) Call attach(File_handle, trim(grpname), "neuttemp", neuttemp) Call attach(File_handle, trim(grpname), "neutpressure", neutpressure) Call attach(File_handle, trim(grpname), "scatter_fac", scatter_fac) Call attach(File_handle, trim(grpname), "Eion", Eion) Call attach(File_handle, trim(grpname), "E0", E0) Call attach(File_handle, trim(grpname), "Escale", Escale) Call putarr(File_handle,trim(grpname)//"species", species) if (allocated(io_cross_sec)) Call putarr(File_handle, trim(grpname)//"/io_cross_sec", io_cross_sec) if (allocated(ela_cross_sec)) Call putarr(File_handle, trim(grpname)//"/ela_cross_sec", ela_cross_sec) END IF End subroutine neutcol_diag !------------------------------------------------------------- SUBROUTINE neutcol_step(plist) ! USE random USE beam USE omp_lib USE basic, ONLY: nlclassical USE distrib, ONLY: lodgaus type(particles), TARGET::plist(:) type(particles),pointer::p INTEGER:: i, omp_thread, num_threads, j, nbcolls_ela, nbcolls_io real(kind=db):: Rand(5) real(kind=db):: v2, v, ek, Everif, es, cosChi, thet, sig_io, sig_ela, vfact, xsi type(linked_part_row), ALLOCATABLE:: ins_p(:) type(linked_part), POINTER:: created, created2 real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3) p=>plist(species(1)) if(.not. nlcol .or. p%nploc .le. 0) return num_threads=omp_get_max_threads() Allocate(ins_p(num_threads)) nbcolls_ela=0 nbcolls_io=0 nucol=0 !$OMP PARALLEL DEFAULT(SHARED), private(collisionfact,i,omp_thread,Rand,v2,ek,sig_io,sig_ela,es,coschi,thet,vfact, created, v, everif,xsi,vinit,vend), reduction(+:nbcolls_ela,nbcolls_io, nucol) omp_thread=omp_get_thread_num()+1 allocate(ins_p(omp_thread)%start) ins_p(omp_thread)%n=0 created=>ins_p(omp_thread)%start !$OMP DO DO i=1,p%Nploc CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread)) v2=(p%UR(i)**2+p%UTHET(i)**2+p%UZ(i)**2) if(nlclassical) THEN ek=v2*escale v=sqrt(v2) vinit=(/p%UR(i),p%UTHET(i),p%UZ(i)/) ELSE ek=(p%gamma(i)-1)*escale v=sqrt(v2)/p%gamma(i) vinit=(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i) end if sig_io=0 sig_ela=0 ! The ionisation event can only occur if the incoming electron energy is above the binding energy if (ek .gt. Eion .and. nb_io_cross .gt. 1) then sig_io=sig_fit(io_cross_sec,io_growth_cross_sec,ek, nb_io_cross) end if if (nb_ela_cross .gt. 1) then sig_ela=sig_fit(ela_cross_sec,ela_growth_cross_sec,ek, nb_ela_cross) xsi=Ek/(0.25*E0+Ek) sig_ela=sig_ela*(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi)) end if collisionfact=1-exp(-collfactor*(sig_io+sig_ela)*v) ! If we have a collision if (Rand(1) .lt.collisionfact) THEN CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread)) IF(Rand(1).gt. sig_ela/(sig_io+sig_ela)) THEN ! An ionisation collision happened and we create the necessary electron ins_p(omp_thread)%n=ins_p(omp_thread)%n+1 allocate(created%next) created%next%prev=>created ! Fill created particle physical quantities created%p%R=p%R(i) created%p%THET=p%THET(i) created%p%Z=p%Z(i) IF( nlmaxwellio ) THEN ! the new electron is created according to a Maxwellian CALL lodgaus(0, Rand(1:3)) ! get random velocity created%p%UR=vth*(Rand(1)) created%p%UTHET=vth*(Rand(2)) created%p%UZ=vth*(Rand(3)) ELSE CALL random_array(Rand,3,ran_index(omp_thread),ran_array(:,omp_thread)) ! Compute created electron energy Es=scatter_fac*tan(Rand(1)*atan((Ek-Eion)/(2*scatter_fac))) ! Compute scattering angles for created electron cosChi=1-2*Rand(2)/(1+8*Es/E0*(1-Rand(2))) thet=Rand(3)*2*pi if(nlclassical)THEN ! new velocity factor for created particle vfact=sqrt(Es/Ek) ELSE ! new velocity factor for created particle vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale))) END IF ! Fill created particle physical quantities created%p%UR=vfact*(p%UR(i)) created%p%UTHET=vfact*(p%UTHET(i)) created%p%UZ=vfact*(p%UZ(i)) call rotate(created%p%UR,created%p%UTHET, created%p%UZ, coschi, thet) END IF vend=(/created%p%UR,created%p%UTHET,created%p%UZ/) if(nlclassical)THEN ! Lorentz factor for created particle created%p%gamma=1.0 ELSE ! Lorentz factor for created particle created%p%gamma=sqrt(1+created%p%UZ**2+created%p%UR**2+created%p%UTHET**2) vend=vend/created%p%gamma END IF ! We prepare the next created particle ins_p(omp_thread)%end=>created created=>created%next ! We keep track of what changed nbcolls_io=nbcolls_io+1 nucol=nucol-vend/vinit ! If we want the incoming electron to be scattered, we need to compute the if (nldragio) THEN ! We store the lossed energy in pot for keeping track of energy conservation created%prev%p%pot=Eion+Es CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread)) Es=Ek-Eion-Es if(nlclassical)THEN ! new velocity factor for scattered particle vfact=sqrt(Es/Ek) ELSE ! new velocity factor for scattered particle vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale))) END IF ELSE CYCLE END IF ELSE ! An elastic collision event happens CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread)) Es=Ek vfact=1.0 nbcolls_ela=nbcolls_ela+1 END IF ! We calculate the scattered velocity angle for the scattered electron cosChi=1-2*Rand(1)/(1+8*Es/E0*(1-Rand(1))) thet=Rand(2)*2*pi ! Change the incident electron velocity direction and amplitude if necessary p%UR(i)=p%UR(i)*vfact p%UTHET(i)=p%UTHET(i)*vfact p%UZ(i)=p%UZ(i)*vfact call rotate(p%UR(i),p%UTHET(i), p%UZ(i), coschi, thet) if(nlclassical) THEN vend=(/p%UR(i),p%UTHET(i),p%UZ(i)/) ELSE p%gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2) vend=(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i) END IF nucol=nucol+1-vend/vinit END IF END DO !$OMP END DO if(associated(created%prev)) then created=>created%prev ins_p(omp_thread)%end=>created deallocate(created%next) else deallocate(ins_p(omp_thread)%start) end if !$OMP END PARALLEL ! We collect all created particules into one linked list Do i=1,num_threads if(associated(ins_p(i)%start)) then created=>ins_p(i)%end Do j=i+1,num_threads created%next=>ins_p(j)%start ins_p(i)%n=ins_p(i)%n+ins_p(j)%n IF(ASSOCIATED(created%next)) then created=>ins_p(j)%end END IF End Do if(species(2).gt.0) then - CALL add_created_part(plist(species(2)), ins_p(i),.false.,.true.) + CALL add_created_part(plist(species(2)), ins_p(i), .false.,.true.) end if CALL add_created_part(p,ins_p(i),.true.,.false.) exit end if end do DEALLOCATE(ins_p) p%nbcolls=p%nbcolls+(/nbcolls_io, nbcolls_ela/) p%nudcol=nucol !Write(*,*)"mpirank: ", mpirank, " Nb colls ela, io: ",nbcolls_ela, nbcolls_io ! END SUBROUTINE neutcol_step FUNCTION sig_fit(sig_vec,growth_vec,ek,nb_cross) use distrib, ONLY: closest real(kind=db)::sig_fit, ek real(kind=db):: sig_vec(:,:), growth_vec(:) Integer:: k, nb_cross sig_fit=0 k=closest(sig_vec(:,1),ek, nb_cross-1) if(k.lt.1) return !sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1) ! Exponential fitting relevant at high energies sig_fit=sig_vec(k,2)*(ek/sig_vec(k,1))**growth_vec(k) END FUNCTION sig_fit SUBROUTINE rotate(Ur, Uthet, Uz, coschi, thet) real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet real(kind=db):: norm, perp(3), U(3), U0(3), normf real(kind=db):: sinchi, sinthet, costhet Integer :: iperp1,iperp2 U0=(/Ur,Uthet,Uz/) norm=sqrt(sum(U0**2)) U=U0/norm ! Find a vector perpendicular to U for chi rotation ! find the direction with maximum amplitude perp=(/1,1,1/) iperp1=maxloc(abs(U),1) ! find second direction with next max amplitude perp(iperp1)=0 iperp2=maxloc(abs(perp*U),1) perp=0 perp(iperp2)=U(iperp1) perp(iperp1)=-U(iperp2) ! Normalise the rotation vector perp=perp/sqrt(sum(perp**2)) ! Compute sinus and cosinus for rotation sinchi=sqrt(1-coschi**2) costhet=cos(thet) sinthet=sin(thet) ! Rotation of angle chi around perp Ur = (coschi+perp(1)**2*(1-coschi))*U0(1) + (perp(1)*perp(2)*(1-coschi)-perp(3)*sinchi)*U0(2) + (perp(1)*perp(3)*(1-coschi) + perp(2)*sinchi)*U0(3) Uthet = (perp(1)*perp(2)*(1-coschi)+perp(3)*sinchi)*U0(1) + (coschi + perp(2)**2*(1-coschi))*U0(2) +(perp(2)*perp(3)*(1-coschi)-perp(1)*sinchi)*U0(3) Uz = (perp(1)*perp(3)*(1-coschi)-perp(2)*sinchi)*U0(1) +(perp(3)*perp(2)*(1-coschi)+perp(1)*sinchi)*U0(2) +( coschi+perp(3)**2*(1-coschi))*U0(3) U0 =(/Ur,Uthet,Uz/) ! second rotation according to uniform distribution ! Rotation of angle theta around U Ur = (costhet+U(1)**2*(1-costhet))*U0(1) + (U(1)*U(2)*(1-costhet) - U(3)*sinthet)*U0(2) + (U(1)*U(3)*(1-costhet)+U(2)*sinthet)*U0(3) Uthet = (U(2)*U(1)*(1-costhet)+U(3)*sinthet)*U0(1) + (costhet + U(2)**2*(1-costhet))*U0(2) + (U(2)*U(3)*(1-costhet)-U(1)*sinthet)*U0(3) Uz = (U(3)*U(1)*(1-costhet) - U(2)*sinthet)*U0(1) + (U(3)*U(2)*(1-costhet)+U(1)*sinthet)*U0(2) + (costhet +U(3)**2*(1-costhet))*U0(3) !normf=sqrt(Ur**2+Uthet**2+Uz**2) !if(abs(norm-normf)/norm .gt. 1e-14) WRITE(*,*) "Error in rotate the norm of v changed" END SUBROUTINE rotate SUBROUTINE read_cross_sec(filename,cross_sec, nb_cross) CHARACTER(len=*) ::filename Real(kind=db), ALLOCATABLE :: cross_sec(:,:) INTEGER:: nb_cross INTEGER :: lu_cross_sec=9999 INTEGER:: i, openerr, reason CHARACTER(len=256) :: header real(kind=db):: t1,t2 nb_cross=0 OPEN(UNIT=lu_cross_sec,FILE=trim(filename),ACTION='READ',IOSTAT=openerr) header=' ' IF(openerr .ne. 0) THEN CLOSE(unit=lu_cross_sec) RETURN END IF ! The cross section table is defined as a two column energy and cross_section DO WHILE(.true.) READ(lu_cross_sec,'(a)',IOSTAT=reason) header header=adjustl(header) if(reason .lt. 0 ) exit ! We reached end of file if( header(1:1) .ne. '!') then READ(header,*) t1, t2 if(t1 .ne. 0 .and. t2.ne. 0) nb_cross=nb_cross+1 end if END DO if (allocated(cross_sec)) deallocate(cross_sec) allocate(cross_sec(nb_cross,2)) REWIND(lu_cross_sec) ! The cross section table is defined as a two column energy and cross_section i=1 DO WHILE(i .le. nb_cross) READ(lu_cross_sec,'(a)',IOSTAT=reason) header header=adjustl(header) if(reason .lt. 0 ) exit ! We reached end of file if( header(1:1) .ne. '!') then READ(header,*) cross_sec(i,1), cross_sec(i,2) if(cross_sec(i,1) .ne. 0 .and. cross_sec(i,2).ne. 0) i=i+1 end if END DO CLOSE(unit=lu_cross_sec) END subroutine read_cross_sec end module neutcol diff --git a/src/psupply_mod.f90 b/src/psupply_mod.f90 index 8727d76..13d3395 100644 --- a/src/psupply_mod.f90 +++ b/src/psupply_mod.f90 @@ -1,309 +1,313 @@ module psupply use constants implicit none type power_supply logical :: active = .false. ! is the power supply active real(kind=db):: geomcapacitor = 1 ! capacitance of the metalic vessel normalised to the neutral density in vessel used in the neutcol module real(kind=db):: PSresistor = 1 ! internal resistance of the power supply normalised to the actual neutral density in vessel real(kind=db):: targetbias = 0 ! Set voltage on the power supply real(kind=db):: bias = 0 ! current voltage on the power supply integer :: nbbounds = 2 ! number of boundaries defined in the geometry integer :: lststp = 0 ! previous step on which the bias was updated real(kind=db):: current(3) = 0 ! current collected on the boundaries normalised to the simulated collision neutral density set in neutcol module ! 1 is at time i-2nbhdt, 2 is at time i-nbhdt and 3 is at time i integer, allocatable:: bdpos(:) ! sign of each boundary for collected charge to determine direction of current real(kind=db),allocatable:: charge(:) ! Charge collected on each boundary and per nbdt real(kind=db),allocatable:: biases(:) ! Actual potentials at each boundary integer :: nbhdt = 10 ! half of the number of time steps between each calls to RK4 real(kind=db):: expdens ! [m-3] experimental neutral density real(kind=db):: neutcoldens ! [m-3] neutral density used in neutcol module end type type(power_supply):: the_ps contains ! read the input parameters from the input file and setup the necesary variables for the module to work subroutine psupply_init(fileid,cstep,nbbounds,neutcoldens,rstbias) use splinebound use basic, only: phinorm, tnorm, mpirank, qnorm, potinn, potout use constants use mpi use geometry use weighttypes use fields integer:: fileid, cstep, nbbounds, istat, nbhdt, ierr,i real(kind=db),OPTIONAL, INTENT(IN):: rstbias real(kind=db):: neutcoldens ! [m-3] real(kind=db):: expneutdens = 1 ! [m-3] real(kind=db):: PsResistor = 1 ! [Ohm] real(kind=db):: geomcapacitor = 1 ! [F] real(kind=db):: targetbias = 0 ! [V] integer, allocatable:: bdpos(:) character(len=1000) :: line logical :: active = .false. NAMELIST /psupplyparams/ expneutdens, PsResistor, geomcapacitor, targetbias, nbhdt, active, bdpos the_ps%lststp=cstep the_ps%nbbounds=nbbounds allocate(the_ps%bdpos(nbbounds),bdpos(nbbounds)) allocate(the_ps%charge(nbbounds),the_ps%biases(nbbounds)) the_ps%bdpos=0 bdpos=0 the_ps%charge=0 the_ps%biases=0 the_ps%current=0 ! read the input parameters from file Rewind(fileid) READ(fileid, psupplyparams, iostat=istat) if (istat.gt.0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in pssupplyparams: '//trim(line) call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if ! save the parameters on output IF(mpirank .eq. 0) THEN WRITE(*, psupplyparams) END IF ! rescale the targetbias set on the power supply the_ps%targetbias=abs(targetbias)/phinorm the_ps%bdpos=bdpos ! save the experimental neutral density the_ps%expdens=expneutdens ! save the neutral collision density the_ps%neutcoldens=neutcoldens if(present(rstbias))then ! Initialize the current bias from the restart value the_ps%bias=rstbias else ! initialize with the file input parameters if (the_domain%nbsplines.gt.0) then do i=1,the_ps%nbbounds if(the_ps%bdpos(i) .lt. 0)then the_ps%bias=-the_domain%boundaries(i)%Dirichlet_val exit end if end do else the_ps%bias=(potout-potinn) end if end if ! set the initial bias where(the_ps%bdpos .lt. 0) the_ps%biases=the_ps%bias*the_ps%bdpos end where ! Normalise resistor and capacitor to adapt to experimental pressure the_ps%PSresistor = PSresistor*the_ps%expdens/the_ps%neutcoldens*qnorm/(tnorm*phinorm) the_ps%geomcapacitor = geomcapacitor*phinorm/qnorm the_ps%nbhdt = nbhdt the_ps%active = active if( .not. the_ps%active) return ! Initialize the biases if (the_domain%nbsplines.gt.0) then do i=1,the_ps%nbbounds the_domain%boundaries(i)%Dirichlet_val=the_ps%biases(i) end do else potinn=the_ps%biases(1) Potout=0 + Phidown=Potinn + Phiup=Potout end if ! recalculate gtilde to adapt for the new biases CALL total_gtilde(vec1, vec2, gtilde, gridwgeom) call comp_gradgtilde end subroutine ! save to the result file the parameters of this module read from the input Subroutine psupply_diag(File_handle, str) use mpi Use futils use basic, only: tnorm, phinorm, qnorm implicit none Integer:: File_handle Character(len=*):: str CHARACTER(len=256):: grpname Integer:: ierr, mpirank CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0 .and. the_ps%active) THEN Write(grpname,'(a,a)') trim(str),"/psupply" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "expdens", the_ps%expdens) Call attach(File_handle, trim(grpname), "targetbias", the_ps%targetbias*phinorm) Call attach(File_handle, trim(grpname), "PSresistor", the_ps%PSresistor/the_ps%expdens*the_ps%neutcoldens/qnorm*(tnorm*phinorm)) Call attach(File_handle, trim(grpname), "geomcapacitor", the_ps%geomcapacitor/phinorm*qnorm) Call attach(File_handle, trim(grpname), "nbhdt", the_ps%nbhdt) Call putarr(File_handle,trim(grpname)//"/bdpos", the_ps%bdpos) END IF End subroutine psupply_diag ! gneral routine called from stepon to update the psupply bias subroutine psupply_step(ps,p,cstep) use particletypes use geometry use weighttypes use fields use basic, only: Potinn, potout type(power_supply):: ps type(particles):: p(:) integer:: cstep, i if (.not. ps%active ) return ! calculate the charge collected on each boundary due to the contribution of each specie call add_charge(ps,p,cstep) ! calculate the current flowing between the electrodes due to the cloud call calc_current(ps,cstep) ! calculate the bias at the new time step call updt_bias(ps,cstep) if(mod(cstep-ps%lststp,2*ps%nbhdt) .ne. 0) return ! update the bias on the geometry for the Dirichlet b.c. if (the_domain%nbsplines.gt.0) then do i=1,ps%nbbounds the_domain%boundaries(i)%Dirichlet_val=ps%biases(i) end do else potinn=ps%biases(1) Potout=0 + Phidown=Potinn + Phiup=Potout end if ! recalculate gtilde to adapt for the new biases CALL total_gtilde(vec1, vec2, gtilde, gridwgeom) call comp_gradgtilde end subroutine ! calculates the current flowing between the electrodes due to the cloud subroutine calc_current(ps,cstep) use geometry use basic, only: phinorm, dt use fields type(power_supply):: ps integer:: cstep if(mod(cstep-ps%lststp,ps%nbhdt).eq.0) then ! communicate the charge accumulation in this timestep call reduce_charge(ps) if(mod(cstep-ps%lststp,ps%nbhdt*2).eq.0)then ! calculates the current by adding the contribution of each boundary if (mpirank .eq. 0)then ps%current(3)=sum(-ps%charge*ps%bdpos)/(ps%nbhdt*dt) end if ps%lststp=cstep else ! calculates the current by adding the contribution of each boundary if (mpirank .eq. 0)then ps%current(2)=sum(-ps%charge*ps%bdpos)/(ps%nbhdt*dt) end if end if ps%charge=0 end if end subroutine ! calculate the charge deposited by each specie on the electrodes (used to calculate the resulting current) subroutine add_charge(ps,p, cstep) use particletypes use basic, only: qnorm type(power_supply):: ps type(particles):: p(:) integer::cstep, i do i=1,size(p,1) if(.not. p(i)%is_field) cycle !Add the normalised contribution of each specie ps%charge=ps%charge+p(i)%nblost(5:)*p(i)%weight*p(i)%q/qnorm end do end subroutine ! Time integrate the ODE of the actual bias between the accelerating electrodes ! and broadcast it to all the workers subroutine updt_bias(ps,cstep) use basic, only: dt, phinorm implicit none type(power_supply):: ps integer:: cstep real(kind=db):: bias,k1,k2,k3,k4, hdeltat if(mod(cstep-ps%lststp,2*ps%nbhdt) .ne. 0) return ! half delta t hdeltat=dt*ps%nbhdt/(ps%PSresistor*ps%geomcapacitor) bias=ps%bias ! we update the bias using RK4 k1=-(bias+ps%current(1)*ps%PSresistor-ps%targetbias) k2=-(bias+hdeltat*k1+ps%current(2)*ps%PSresistor-ps%targetbias) k3=-(bias+hdeltat*k2+ps%current(2)*ps%PSresistor-ps%targetbias) k4=-(bias+2*hdeltat*k3+ps%current(3)*ps%PSresistor-ps%targetbias) ps%bias=bias+(k1+2*k2+2*k3+k4)*2*hdeltat/6 Write(*,*) " new bias ", ps%bias*phinorm where (ps%bdpos .lt. 0) ps%biases=-ps%bias end where ! broadcast the bias to all the mpi processes call bcast_bias(ps) ps%current(1)=ps%current(3) end subroutine updt_bias ! gather on node 0 the collected charge on each metallic boundary subroutine reduce_charge(ps) use mpi use mpihelper use basic, ONLY: mpirank type(power_supply):: ps integer:: ierr if(mpirank .eq. 0) then call MPI_REDUCE(MPI_IN_PLACE,ps%charge,ps%nbbounds,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) Write(*,*) "curr charge ", ps%charge else call MPI_REDUCE(ps%charge,ps%charge,ps%nbbounds,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) end if end subroutine ! broadcast to all the nodes the new bias imposed by the power supply on the electrodes subroutine bcast_bias(ps) use mpi use mpihelper type(power_supply):: ps integer:: ierr call MPI_BCAST(ps%biases,ps%nbbounds,db_type,0,MPI_COMM_WORLD,ierr) end subroutine end module psupply diff --git a/src/resume.f90 b/src/resume.f90 index b59b155..72554bd 100644 --- a/src/resume.f90 +++ b/src/resume.f90 @@ -1,76 +1,73 @@ SUBROUTINE resume ! ! Resume from previous run ! Use beam Use basic Use fields Use sort Use maxwsrce Use geometry Use neutcol IMPLICIT NONE ! ! Local vars and arrays INTEGER:: i !________________________________________________________________________________ WRITE(*,'(a)') ' Resume from previous run' !________________________________________________________________________________ ! ! Open and read initial conditions from restart file CALL chkrst(0) ! If we want to start a new run with a new file IF (newres) THEN ! we start time and step from 0 cstep=0 time=0 END IF qnorm=abs(partslist(1)%weight*partslist(1)%q) CALL fields_init call timera(0, "read_geom") call read_geom(lu_in, rnorm, splrz, Potinn, Potout) call timera(1, "read_geom") ! Initialize the magnetic field and the electric field solver call mag_init CALL fields_start call boundary_loss(partslist(1)) - call localisation(partslist(1)) ! Add the requested test particles and read them from input files IF (nbaddtestspecies .gt. 0) THEN Do i=1,nbaddtestspecies CALL load_part_file(partslist(nbspecies+i),addedtestspecfile(i)) call boundary_loss(partslist(nbspecies+i)) - call localisation(partslist(nbspecies+i)) END DO nbspecies=nbspecies+nbaddtestspecies END IF IF(mpisize .gt. 1) THEN CALL calc_Zbounds(partslist(1), Zbounds, femorder) DO i=1,nbspecies CALL keep_mpi_self_parts(partslist(i),Zbounds) call collectparts(partslist(i)) END DO END IF CALL bound(partslist(1)) call boundary_loss(partslist(1)) - CALL localisation(partslist(1)) ! Allocate the variables for MPI communications CALL fields_comm_init(Zbounds) ! Solve Poisson for the initial conditions CALL rhscon(partslist) CALL poisson(splrz) ! Compute the fields at the particles positions CALL EFieldscompatparts(partslist(1)) if(mpirank .eq. 0) WRITE(*,*) "Initial forces computed" ! END SUBROUTINE resume diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index 845c2fd..cd8fa4f 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,806 +1,771 @@ MODULE splinebound USE constants USE bsplines USE forSISL, only: newCurve, freeCurve, freeIntCurve, writeSISLcurve, writeSISLpoints Use forSISLdata IMPLICIT NONE INTEGER(SELECTED_INT_KIND(4)), PARAMETER :: bd=-1, bd_Dirichletconst=0, bd_Dirichletvar=1, bd_Neumann=2 type cellkind integer:: spldirkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate integer:: spltotkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate integer:: linkedboundaries(2)=0 !< stores the spline curve indices in the spline_domain of the spline boundaries that are at a distance lower than dist_extent integer:: leftknot(2)=0 !< for each linkedboundary saves a lower initial guess for calculating the distance real(kind=db):: lguess(2)=-1 !< Spline parameter left limit as start guess real(kind=db):: rguess(2)=-1 !< Spline parameter right limit as start guess end type cellkind TYPE spline_boundary ! all curves assume right handedness to set which side of the curve is inside or outside type(SISLCurve):: curve Real(kind=db):: Dirichlet_val !< Value for the dirichlet boundary condition created by this boundary Real(kind=db):: epsge=1.0e-5 !< geometric resolution used for calculating distances Real(kind=db):: epsce=1.0e-9 !< value of weight below which it is 0 INTEGER(kind(bd)):: type=bd_Dirichletconst !< type of boundary conditions END TYPE spline_boundary type spline_domain integer:: nbsplines = 0 !< number of spline boundaries in the domain type(spline_boundary), allocatable:: boundaries(:) !< List of boundaries in the domain Real(kind=db):: dist_extent=0.1 !< distance used for the merging with the plateau function for the weight type(cellkind), ALLOCATABLE:: cellk(:,:) !< Precomputed parameters at each cell for faster weight computation type(spline2d), pointer:: splrz => null() !< Pointer to the main spline grid used for the FEM solver Integer:: nb1 !< Number of grid points in the 1st dimension Integer:: nb2 !< Number of grid points in the 2nd dimension real(kind=db), ALLOCATABLE:: x1(:) !< Grid points in first direction for weight interpolation real(kind=db), ALLOCATABLE:: x2(:) !< Grid points in 2nd direction for weight interpolation real(kind=db), ALLOCATABLE:: dx1(:) !< inverse cell width in first direction for weight interpolation real(kind=db), ALLOCATABLE:: dx2(:) !< inverse cell width in 2nd direction for weight interpolation type(SISLsurf):: Dirdomweight !< structure storing precalculated geometric weight for faster evaluation type(SISLsurf):: totdomweight !< structure storing precalculated total weight for faster evaluation end type spline_domain CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Reads a spline domain from the namelist or from a h5 file. Needs to be called for the initialization of the module !> @param[in] Fileid Text file id of the input file containing namelists !> @param[out] spldom spline domain !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] rnorm distance normalization constant !> @param[in] phinorm electric potential normalization constant !--------------------------------------------------------------------------- subroutine read_splinebound(Fileid, spldom, splrz, rnorm, Phinorm) use mpi Integer:: Fileid type(spline_domain):: spldom type(spline2d):: splrz real(kind=db):: rnorm, phinorm Integer:: nbsplines, istat, mpirank, ierr - real(kind=db):: dist_extent, Potinn, Potout, epsge, epsce + real(kind=db):: dist_extent, Potinn, Potout Character(len=128):: h5fname="", line real(kind=db) :: Dvals(30)=0 integer:: nbpoints - real(kind=db),allocatable:: cpoints(:,:) - real(kind=db):: dz, dr, ctr, Lz, ra, rb integer:: degree=4, i - namelist /spldomain/ nbsplines, dist_extent, h5fname, dr, Lz, ra, rb, epsge, epsce, Dvals + namelist /spldomain/ nbsplines, dist_extent, h5fname, Dvals CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) REWIND(fileid) READ(fileid,spldomain, iostat=istat) if (istat.gt.0) then if(mpirank .eq. 0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in geomparams: '//trim(line) end if call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(mpirank .eq. 0) WRITE(*, spldomain) Dvals=Dvals/phinorm dist_extent=dist_extent/rnorm if (.not. trim(h5fname)=='' ) then call setspline_domain(spldom, splrz, dist_extent, 0) call splinebound_readh5domain(h5fname,spldom, rnorm, phinorm) call classifycells(spldom) do i=1,spldom%nbsplines spldom%boundaries(i)%Dirichlet_val=Dvals(i) end do return + else + WRITE(*,*) "Error the filename h5fname is not defined. No boundary has been set!" + call mpi_Abort(MPI_COMM_WORLD, -1, ierr) end if - - nbpoints=600 - allocate(cpoints(2,nbpoints)) - - - nbsplines=2 - call setspline_domain(spldom, splrz, dist_extent, nbsplines) - - dz=(splrz%sp1%knots(splrz%sp1%nints)-splrz%sp1%knots(0))/(nbpoints-1) - Lz=Lz/rnorm - ctr=(splrz%sp1%knots(splrz%sp1%nints)+splrz%sp1%knots(0))/2 - do i=1,nbpoints - cpoints(1,i)=splrz%sp1%knots(0)+(i-1)*dz - end do - - cpoints(2,:)=ra/rnorm - call setspline_boundary(spldom%boundaries(1), cpoints, degree, Dvals(1), epsge, epsce) - - !spldom%boundaries(1)%type=bd_Neumann - - do i=1,nbpoints - cpoints(1,i)=splrz%sp1%knots(0)+(nbpoints-i)*dz - end do - dr=dr/rnorm - cpoints(2,:)=rb/rnorm - do i=1,nbpoints - if(((cpoints(1,i)-ctr)/Lz)**2.le.1) then - cpoints(2,i)=rb/rnorm+dr*sqrt(1-(cpoints(1,i)-ctr)**2/Lz**2) - endif - end do - - call setspline_boundary(spldom%boundaries(2), cpoints, degree, Dvals(2), epsge, epsce) - do i=1,spldom%nbsplines - spldom%boundaries(i)%Dirichlet_val=Dvals(i) - end do - call classifycells(spldom) end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Saves the spline boundaries to the result file !> @param[in] File_handle futils h5 file id !> @param[in] curr_grp groupname under which the boundaries must be saved !> @param[in] spldom spline domain !--------------------------------------------------------------------------- Subroutine splinebound_diag(File_handle, curr_grp, spldom) use mpi Use futils Use basic, ONLY: rnorm, phinorm Integer:: File_handle type(spline_domain):: spldom Character(len=*):: curr_grp CHARACTER(len=128):: grpname Integer:: ierr, mpirank, i CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0) THEN Write(grpname,'(a,a)') trim(curr_grp),"/geometry_spl" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "dist_extent",spldom%dist_extent) Call attach(File_handle, trim(grpname), "nbsplines", spldom%nbsplines) do i=1,spldom%nbsplines Write(grpname,'(a,a,i2.2)') trim(curr_grp),"/geometry_spl/",i If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "Dirichlet_val", spldom%boundaries(i)%Dirichlet_val*phinorm) Call attach(File_handle, trim(grpname), "order", spldom%boundaries(i)%curve%ik) Call attach(File_handle, trim(grpname), "kind", spldom%boundaries(i)%curve%ikind) Call attach(File_handle, trim(grpname), "dim", spldom%boundaries(i)%curve%idim) CALL putarr(File_handle, TRIM(grpname)//"/pos", spldom%boundaries(i)%curve%ecoef*rnorm) CALL putarr(File_handle, TRIM(grpname)//"/knots", spldom%boundaries(i)%curve%et) end do END IF End subroutine splinebound_diag !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Read a spline boundary domain from an h5 file structure !> @param[out] spldom new spline domain !> @param[in] filename filename of the h5 file !> @param[in] rnorm distance normalization constant !> @param[in] phinorm electric potential normalization constant !--------------------------------------------------------------------------- subroutine splinebound_readh5domain(filename, spldom, rnorm, phinorm) use futils use forSISL implicit none Character(len=*),intent(in) :: filename type(spline_domain),intent(inout) :: spldom integer:: h5id, i real(kind=db):: rnorm, phinorm CHARACTER(len=128):: grpname integer:: periodic integer:: order, dim, bdtype INTEGER:: posrank, posdim(2), err real(kind=db):: Dval, epsge, epsce real(kind=db),allocatable:: points(:,:) call openf(filename, h5id,'r','d') call getatt(h5id, '/geometry_spl/','nbsplines', spldom%nbsplines) ! prepare memory if (allocated(spldom%boundaries)) then do i=1,size(spldom%boundaries,1) call free_bsplinecurve(spldom%boundaries(i)) end do DEALLOCATE(spldom%boundaries) end if allocate(spldom%boundaries(spldom%nbsplines)) ! Read each boundary curve individually do i=1,spldom%nbsplines Write(grpname,'(a,i2.2)') "/geometry_spl/",i If(.not. isgroup(h5id, trim(grpname))) THEN Write(*,*) "Error the geometry definition file is invalid" END IF periodic=0 Call getatt(h5id, trim(grpname), "Dirichlet_val", Dval) Call getatt(h5id, trim(grpname), "epsge", epsge) Call getatt(h5id, trim(grpname), "epsce", epsce) Call getatt(h5id, trim(grpname), "order", order) Call getatt(h5id, trim(grpname), "dim", dim) err=0 Call getatt(h5id, trim(grpname), "periodic", periodic,err) if(err .lt.0) periodic=0 CALL getdims(h5id, TRIM(grpname)//"/pos", posrank, posdim) allocate(points(posdim(1),posdim(2))) CALL getarr(h5id, TRIM(grpname)//"/pos", points) points=points/rnorm Call setspline_boundary(spldom%boundaries(i),transpose(points), order-1, Dval/phinorm, epsge,epsce, periodic) bdtype=bd err=0 Call getatt(h5id, trim(grpname), "type", bdtype,err) if(err.ge.0) spldom%boundaries(i)%type=bdtype deallocate(points) end do call closef(h5id) end subroutine splinebound_readh5domain !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> initialize a spline domain and allocate the necessary memory !> @param[out] spldom new spline domain !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] dist_extent normalized characteristic fall lenght of the weight !> @param[in] nb_splines number of boundary splines to allocate !--------------------------------------------------------------------------- subroutine setspline_domain(spldom,splrz,dist_extent, nb_splines) type(spline_domain):: spldom type(spline2d), TARGET:: splrz real(kind=db):: dist_extent integer:: nb_splines, nb1, nb2 ! Store the grid parameters to speed-up calculations nb1=splrz%sp1%nints nb2=splrz%sp2%nints spldom%nb1=nb1 spldom%nb2=nb2 spldom%splrz=>splrz allocate(spldom%cellk(0:nb1-1,0:nb2-1)) allocate(spldom%x1(0:nb1)) allocate(spldom%x2(0:nb2)) allocate(spldom%dx1(0:nb1-1)) allocate(spldom%dx2(0:nb2-1)) spldom%x1(0:)=splrz%sp1%knots(0:nb1) spldom%x2(0:)=splrz%sp2%knots(0:nb2) spldom%dx1(0:)=1/(spldom%x1(1:nb1)-spldom%x1(0:nb1-1)) spldom%dx2(0:)=1/(spldom%x2(1:nb2)-spldom%x2(0:nb2-1)) !Prepare structures to host singular spline boundaries spldom%nbsplines=nb_splines if(spldom%nbsplines.gt. 0) allocate(spldom%boundaries(nb_splines)) spldom%dist_extent=dist_extent end subroutine setspline_domain subroutine setspline_boundary(b_curve, cpoints, degree, D_val, epsge, epsce, periodic) Use bsplines use forSISL,ONLY: newcurve, s1630 use mpi type(spline_boundary):: b_curve Real(kind=db):: cpoints(:,:) Real(REAL64),ALLOCATABLE:: points(:) Real(REAL64):: astpar integer:: degree integer, optional:: periodic Integer:: order, ierr Real(kind=db):: D_val Real(kind=db),OPTIONAL :: epsge, epsce Integer:: nbpoints,i, dim, jstat, bsptype integer:: period period=0 if(present(periodic))period=periodic nbpoints= size(cpoints,2) dim=size(cpoints,1) order=degree+1 if(nbpoints .lt. order) then WRITE(*,'(a,i3,a,i5)') "Error: the number of points", nbpoints, " is insuficient for the required order ", order CALL mpi_finalize(ierr) call EXIT(-1) end if allocate(points(dim*nbpoints)) points=reshape(cpoints,(/dim*nbpoints/)) !CALL newCurve(nbpoints, order, knots, points, kind, dim, copy, b_curve%curve) bsptype=1 ! open boundaries b-spline if(period.gt.0) bsptype=-1 ! closed periodic curve astpar=0.0 ! starting parameter for the knots vector CALL s1630(points, nbpoints, astpar, bsptype, dim, order, b_curve%curve, jstat) if (jstat > 0 ) WRITE(*,*) "Warning ", jstat," in curve initialisation s1630 for splineweight" if (jstat < 0 ) WRITE(*,*) "Error ", jstat," in curve initialisation s1630 for splineweight" b_curve%Dirichlet_val=D_val if(present(epsge)) b_curve%epsge=epsge if(present(epsce)) b_curve%epsce=epsce end subroutine setspline_boundary !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the Dirichlet boundary weight from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] x1(:) array of axial positions where the weights are evaluated !> @param[in] x2(:) array of radial positions where the weights are evaluated !> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_w(spldom,x1,x2,w) use forSISL,ONLY: s1424 type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(:,0:) Integer:: k,sstat,l type(cellkind):: cellk REAL(real64):: wtmp(4), point(2) Integer,allocatable::i(:),j(:) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) Do k=1,size(x2,1) cellk=spldom%cellk(i(k),j(k)) SELECT CASE (cellk%spldirkind) CASE(-1) w(k,:)=0 w(k,0)=-1 CYCLE CASE(1) w(k,:)=0 w(k,0)=1 CYCLE CASE DEFAULT !call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,rguess=cellk%rguess(1),lguess=cellk%lguess(1)) point=(/x1(k),x2(k)/) call s1424(spldom%Dirdomweight,1,1,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat) if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splineweight" if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splineweight" w(k,:)=wtmp(1:size(w,2)) END SELECT end DO End SUBROUTINE spline_w !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the total geometric weight from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] x1(:) array of axial positions where the weights are evaluated !> @param[in] x2(:) array of radial positions where the weights are evaluated !> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_wtot(spldom,x1,x2,w,idwall) use forSISL,ONLY: s1424 type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(:,0:) INTEGER, optional, INTENT(OUT):: idwall(:) Integer:: k,sstat,l type(cellkind):: cellk REAL(real64):: wtmp(4), point(2) Integer,allocatable::i(:),j(:) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) Do k=1,size(x2,1) cellk=spldom%cellk(i(k),j(k)) if(present(idwall)) idwall(k)=cellk%linkedboundaries(2) SELECT CASE (cellk%spltotkind) CASE(-1) w(k,:)=0 w(k,0)=-1 CYCLE CASE(1) w(k,:)=0 w(k,0)=1 CYCLE CASE DEFAULT !call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,rguess=cellk%rguess(1),lguess=cellk%lguess(1)) point=(/x1(k),x2(k)/) if(size(w,2).gt.1) then call s1424(spldom%totdomweight,1,1,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat) else call s1424(spldom%totdomweight,0,0,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat) end if if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splinedomweight" if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splinedomweight" w(k,:)=wtmp(1:size(w,2)) END SELECT end DO End SUBROUTINE spline_wtot !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the interpolation in the domain of the Dirichlet boundary conditions from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] z(:) array of axial positions where the weights are evaluated !> @param[in] r(:) array of radial positions where the weights are evaluated !> @param[out] g(:,0:) matrix of boundary interpolations g with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_g(spldom,x1,x2,g,w) use forSISL,ONLY: s1424 type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: g(:,0:) Real(kind=db), INTENT(IN),OPTIONAL::w(:,0:) REAL(real64),allocatable:: gtmp(:,:) Integer:: k,sstat,sizew Integer,allocatable::i(:),j(:) type(cellkind):: cellk allocate(gtmp(size(x2,1),0:size(g,2)-1)) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) if(present(w)) then gtmp=w else Do k=1,size(x2,1) cellk=spldom%cellk(i(k),j(k)) call s1424(spldom%Dirdomweight,1,1,(/x1(k),x2(k)/),cellk%leftknot(1),cellk%leftknot(2),gtmp(k,:),sstat) if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splineweight" if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splineweight" end DO end if Do k=1,size(x2,1) cellk=spldom%cellk(i(k),j(k)) if(cellk%spldirkind.eq.0)then if(gtmp(k,0) .ge. 0) then if(size(g,2) .gt. 1) then g(k,1:2)=-gtmp(k,1:2)*spldom%boundaries(cellk%linkedboundaries(1))%Dirichlet_val end if g(k,0)=(1-gtmp(k,0))*spldom%boundaries(cellk%linkedboundaries(1))%Dirichlet_val else g(k,0)=spldom%boundaries(cellk%linkedboundaries(1))%Dirichlet_val if(size(g,2).gt. 1) then g(k,1:2)=0 end if end if else g(k,:)=0 end if end DO End SUBROUTINE spline_g !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Evaluates the geometric weight induced by the spline curve defined by b_curve at position (z,r) !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] z axial position where the weight is evaluated !> @param[in] r radial position where the weight is evaluated !> @param[out] weight(:) weight index defines the order of derivation by r or z !> @param[in] h distance from the spline at which the weight is 1 !> @param[out] distance unscaled distance between evaluation point and spline b_curve !> @param[inout] leftknot initial guess for the closest spline knot of the points (r,z) !--------------------------------------------------------------------------- subroutine splineweight(b_curve, z, r, weight, h, distance, guess, lguess, rguess) Use forSISL, ONLY: s1227,s1221, s1774 type(spline_boundary):: b_curve Real(kind=db)::r,z Real(kind=db):: weight(0:) Real(kind=db),OPTIONAL:: distance real(kind=db),OPTIONAL:: guess real(kind=db),OPTIONAL:: lguess real(kind=db),OPTIONAL:: rguess integer:: sstatus, der, left,siz real(kind=db):: h, d, tpos, proj real(kind=real64):: curvepos(2*b_curve%curve%idim) real(kind=db):: leftpar, rightpar,guesspar weight=0 der=1 sstatus=-1 guesspar=-1.0_db if(present(lguess) .and. present(rguess)) then leftpar=lguess rightpar=rguess guesspar=(lguess+rguess)/2 call s1774(b_curve%curve,(/z,r/),b_curve%curve%idim,b_curve%epsge,leftpar,rightpar,guesspar,tpos,sstatus) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1774 for splineweight at ", z, r else call dist(b_curve,(/z,r/),d,tpos) end if ! position and derivative wrt r,z call s1221(b_curve%curve,der,tpos,left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1227 for splineweight at ", z, r if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1227 for splineweight at ", z, r d=sqrt((curvepos(1)-z)**2+(curvepos(2)-r)**2) weight(0)=1-max((h-d)/h,0.0_db)**3 curvepos(3:4)=curvepos(3:4)/sqrt(curvepos(3)**2+curvepos(4)**2) ! if the projection of the distance vector on the normal is negative, the weight is negative proj=(-(z-curvepos(1))*curvepos(4)+(r-curvepos(2))*curvepos(3)) if (proj .lt. 0 .or. abs(abs(proj) -sqrt((z-curvepos(1))**2+(r-curvepos(2))**2)).gt.1e-10) weight(0)=-weight(0) !if (proj .lt. 0 ) weight(0)=-weight(0) siz=size(weight,1) if (size(weight,1).gt.1 .and. abs(weight(0)) .lt. 1) then weight(1)=-3*curvepos(4)*(h-d)**2/h**3 weight(2)=+3*curvepos(3)*(h-d)**2/h**3 end if if(present(distance)) distance=d if(present(guess)) guess=tpos end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the closest distance between the point and the selected spline b_curve !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] point(:) array containing the position from which to calculate the distance !> @param[out] distance distance from the point to the spline !> @param[in] pos parameter value of the closest point on the spline !--------------------------------------------------------------------------- subroutine dist(b_curve, point, distance, pos) Use forSISL, ONLY: s1957,s1953, s1221,s1227 type(spline_boundary):: b_curve Real(kind=db):: point(:) real(kind=db):: distance Real(kind=db),optional::pos REAL(real64):: posres, epsco, epsge,curvepos(2),d,distmin REAL(real64),allocatable::intpar(:) integer:: numintpt, sstat, numintcu,i,left,sstatus type(SISLIntCurve),ALLOCATABLE:: intcurve(:) epsco=1.0e-15 epsge=1.0e-15 !epsco=0 !epsge=b_curve%epsge numintpt=0 sstatus=0 distmin=HUGE(d) call s1953(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,numintpt,intpar,numintcu,intcurve,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1953 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1953 for splineweight at ",point(1), point(2) if(numintpt .gt. 1) then Do i=1,numintpt call s1227(b_curve%curve,0,intpar(i),left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1221 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1221 for splineweight at ",point(1), point(2) d=(curvepos(1)-point(1))**2+(curvepos(2)-point(2))**2 if(d .lt. distmin) then distmin=d posres=intpar(i) end if end do else if(numintpt .gt. 0) then posres=intpar(1) end if if(numintcu.ge.1) then posres=0.5*(intcurve(1)%epar1(1)+intcurve(1)%epar1(2)) end if if (present(pos)) pos=posres END subroutine SUBROUTINE classify(x1, x2, cellk, spldom, wpredir, wpretot) real(kind=db), INTENT(IN):: x2(2), x1(2) type(cellkind), intent(INOUT):: cellk type(spline_domain)::spldom Real(kind=db):: zeval(4),reval(4), wpretot, wpredir real(kind=db), allocatable:: guess(:,:), w(:,:,:) Real(kind=db):: dmin, insidedir, insidetot, distance integer:: i,k allocate(guess(spldom%nbsplines,4)) allocate(w(0:2,spldom%nbsplines,4)) w=0 cellk%spldirkind=0 guess=-1.0_db dmin=HUGE(spldom%dist_extent) cellk%linkedboundaries=0 zeval=(/ x1(1),x1(2),x1(1),x1(2) /) reval=(/ x2(1),x2(1),x2(2),x2(2) /) insidedir=1 do i=1,spldom%nbsplines do k=1,4 ! calculate the weight for each spline boundaries at each cell corner call splineweight(spldom%boundaries(i),zeval(k),reval(k),w(:,i,k),spldom%dist_extent,distance,guess(i,k)) ! We find the closest boundary to this point if(distance .lt. dmin) then ! If we are close enough we check if we are below dist_extent and need to calculate the distance each time if(distance .lt. spldom%dist_extent) then if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then cellk%linkedboundaries(1)=i cellk%spldirkind=0 end if cellk%linkedboundaries(2)=i cellk%spltotkind=0 end if dmin=distance ! Otherwise we define the interior by the closest spline if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then insidedir=w(0,i,k) end if insidetot=w(0,i,k) end if end do end do if(cellk%linkedboundaries(1) .gt. 0) then i=cellk%linkedboundaries(1) cellk%lguess(1)=minval(guess(i,:),1,guess(i,:).ge.0) cellk%rguess(1)=maxval(guess(i,:),1) wpredir=w(0,i,1) else cellk%spldirkind=sign(1,int(insidedir)) wpredir=insidedir end if if(cellk%linkedboundaries(2) .gt. 0) then i=cellk%linkedboundaries(2) wpretot=w(0,i,1) else cellk%spltotkind=sign(1,int(insidetot)) wpretot=insidetot end if end subroutine subroutine classifycells(spldom) use forSISL, ONLY: s1537, s1424 type(spline_domain):: spldom integer:: i,j,sstat type(cellkind):: cellk real(kind=db), allocatable:: wpretot(:,:,:), wpredir(:,:,:) allocate(wpretot(1:1,0:spldom%nb1,0:spldom%nb2)) allocate(wpredir(1:1,0:spldom%nb1,0:spldom%nb2)) wpretot=0 wpredir=0 !$OMP PARALLEL DO private(i,j) do i=0,spldom%nb1-1 !DIR$ UNROLL do j=0,spldom%nb2-1 call classify(spldom%x1(i:i+1),spldom%x2(j:j+1),spldom%cellk(i,j),spldom, wpredir(1,i,j),wpretot(1,i,j)) end do end do !$OMP END PARALLEL DO i=spldom%nb1 !$OMP PARALLEL DO private(j,cellk) do j=0,spldom%nb2 cellk=spldom%cellk(min(i,spldom%nb1-1),min(j,spldom%nb2-1)) If(abs(cellk%spldirkind) .eq. 1) Then wpredir(1,i,j)=cellk%spldirkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), spldom%x1(i),spldom%x2(j), wpredir(:,i,j),spldom%dist_extent) end IF If(abs(cellk%spltotkind) .eq. 1) Then wpretot(1,i,j)=cellk%spltotkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(2)), spldom%x1(i),spldom%x2(j), wpretot(:,i,j),spldom%dist_extent) end IF end do !$OMP END PARALLEL DO !$OMP PARALLEL DO private(i,cellk) do i=0,spldom%nb1 j=spldom%nb2 cellk=spldom%cellk(min(i,spldom%nb1-1),min(j,spldom%nb2-1)) If(abs(cellk%spldirkind) .eq. 1) Then wpredir(1,i,j)=cellk%spldirkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), spldom%x1(i),spldom%x2(j), wpredir(:,i,j),spldom%dist_extent) end IF If(abs(cellk%spltotkind) .eq. 1) Then wpretot(1,i,j)=cellk%spltotkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(2)), spldom%x1(i),spldom%x2(j), wpretot(:,i,j),spldom%dist_extent) end IF end do !$OMP END PARALLEL DO call s1537(reshape(wpredir(1,:,:),(/(spldom%nb1+1)*(spldom%nb2+1)/)),spldom%nb1+1,spldom%nb2+1,1,spldom%x1,spldom%x2,0,0,0,0,4,4,1,1,spldom%Dirdomweight,sstat) if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1537 for splineweight" if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1537 for splineweight" call s1537(reshape(wpretot(1,:,:),(/(spldom%nb1+1)*(spldom%nb2+1)/)),spldom%nb1+1,spldom%nb2+1,1,spldom%x1,spldom%x2,0,0,0,0,4,4,1,1,spldom%totdomweight,sstat) if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1537 for splineweight" if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1537 for splineweight" end subroutine subroutine getindex(x1,x2,spldom, i, j) use distrib, ONLY: closest type(spline_domain):: spldom real(kind=db):: x1(:), x2(:) integer:: i(:),j(:) call locintv(spldom%splrz%sp1,x1, i) call locintv(spldom%splrz%sp2,x2, j) end subroutine subroutine free_bsplinecurve(b_curve) type(spline_boundary):: b_curve call freeCurve(b_curve%curve) !call freeIntCurve(b_curve%intcurve) end subroutine END MODULE splinebound diff --git a/src/stepon.f90 b/src/stepon.f90 index 03951ad..c3ccd90 100644 --- a/src/stepon.f90 +++ b/src/stepon.f90 @@ -1,101 +1,99 @@ SUBROUTINE stepon ! ! Advance one time step ! USE basic USE constants USE fields USE beam USE maxwsrce USE celldiag USE neutcol USE sort Use psupply INTEGER:: i DO i=1,nbspecies ! Boundary conditions for plasma particles outside the plasma region CALL bound(partslist(i)) - call boundary_loss(partslist(i)) - ! Localisation of particles in cells (calculation of the r and z indices) - CALL localisation(partslist(i)) + call boundary_loss(partslist(i)) END DO ! Cell diag quantities IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN CALL celldiag_save(time, fidres) END IF ! We compute collisions on the main particles IF(modulo(step,itcol).eq. 0) THEN CALL neutcol_step(partslist) END IF ! The particles are injected by the source CALL maxwsrce_inject(time) ! Sort particles for faster rhscon run time ! DO i=1,nbspecies ! IF(modulo(step,it2d) .eq. 0) THEN ! CALL gridsort(partslist(i),1,partslist(i)%Nploc) ! END IF ! END DO ! Assemble right hand side of Poisson equation CALL rhscon(partslist) if (.not. nlfreezephi) THEN ! Solve Poisson equation CALL poisson(splrz) end if DO i=1,nbspecies ! Compute the electric field at the particle position CALL EFieldscompatparts(partslist(i)) ! Compute the magnetic field at the particle position call comp_mag_p(partslist(i)) ! Solve Newton eq. and advance velocity by delta t CALL comp_velocity(partslist(i)) ! Compute the energy of added particles CALL calc_newparts_energy(partslist(i)) END DO ! Calculate main physical quantities CALL partdiagnostics IF (modulo(step,it2d).eq. 0 .or. nlend) THEN Do i=1,nbspecies if(partslist(i)%calc_moments) CALL momentsdiag(partslist(i)) End do END IF ! update the power supply voltage if necessary call psupply_step(the_ps,partslist,cstep) ! Save variables to file CALL diagnose(step) Do i=1,nbspecies ! Calculate new positions of particles at time t+delta t CALL push(partslist(i)) END DO ! We recalculate the mpi axial boundaries and we adapt them if necessary IF(modulo(step,50) .eq. 0) THEN CALL calc_Zbounds(partslist(1),Zbounds, femorder) CALL fields_comm_init(Zbounds) CALL maxwsrce_calcfreq(Zbounds) END IF END SUBROUTINE stepon