diff --git a/cmake/FindforSISL.cmake b/cmake/FindforSISL.cmake index da051cc..6fa409a 100644 --- a/cmake/FindforSISL.cmake +++ b/cmake/FindforSISL.cmake @@ -1,6 +1,6 @@ find_library(forSISL_LIBRARY NAMES forsisl.a PATHS ${forSISL}/lib) -find_path(forSISL_INCLUDE_DIR forsisl.mod ${forSISL}/include) +find_path(forSISL_INCLUDE_DIR forsisl.mod ${forSISL}/include ${forSISL}/modules) mark_as_advanced(forSISL_LIBRARY forSISL_INCLUDE_DIR) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(forSISL DEFAULT_MSG forSISL_LIBRARY) diff --git a/cmake/Findsisl.cmake b/cmake/Findsisl.cmake index a81903b..0e6a502 100644 --- a/cmake/Findsisl.cmake +++ b/cmake/Findsisl.cmake @@ -1,7 +1,8 @@ +message(STATUS ${SISL}) find_library(sisl_LIBRARY NAMES sisl PATHS ${SISL}/lib) find_path(sisl_INCLUDES sisl.h sislP.h ${SISL}/include) mark_as_advanced(sisl_LIBRARY sisl_INCLUDES) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(sisl DEFAULT_MSG sisl_LIBRARY) diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index 77f1bb3..cb6c465 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,2299 +1,2300 @@ MODULE beam !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Guillaume Le Bars EPFL/SPC !> Patryk Kaminski EPFL/SPC !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> Module responsible for loading, advancing and computing the necessary diagnostics for the simulated particles. !------------------------------------------------------------------------------ ! USE constants use mpi USE mpihelper USE basic, ONLY: mpirank, mpisize USE distrib USE particletypes USE weighttypes IMPLICIT NONE ! !TYPE(particles) :: parts !< Storage for all the particles !SAVE :: parts TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist_towrite ! Diagnostics (scalars) REAL(kind=db) :: ekin=0 !< Total kinetic energy (J) REAL(kind=db) :: epot=0 !< Total potential energy (J) REAL(kind=db) :: etot=0 !< Current total energy (J) REAL(kind=db) :: etot0=0 !< Initial total energy (J) REAL(kind=db) :: loc_etot0=0 !< theoretical local total energy (J) REAL(kind=db) :: Energies(4) !< (1) kinetic energy, (2) potential energy, (3) total energy and (4) gained/lossed energy due to gain or loss of particles (J) ! INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process INTERFACE add_created_part MODULE PROCEDURE add_linked_created_part, add_list_created_part END INTERFACE add_created_part ! abstract interface subroutine rloader(nbase,y,rminus,rplus) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rplus, rminus end subroutine REAL(kind=db) FUNCTION gamma(UZ, UR, UTHET) USE constants REAL(kind=db), INTENT(IN):: UR,UZ,UTHET end FUNCTION end interface CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Loads the particles at the beginning of the simulation and create the parts variable if necessary !--------------------------------------------------------------------------- SUBROUTINE load_parts USE basic, ONLY: nplasma, mpirank, ierr, distribtype, nlclassical, nbspecies, partfile use mpi INTEGER:: i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V, uniform in Z and 1/R in R CALL loaduniformRZ(partslist(1), VR, VZ, VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(3) !Stable distribution from Davidson 4.95 p.119 but with constant distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodinvr) CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodgausr) CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz @brief Checks for each particle if the z position is outside of the local/global simulation space. !> Depending on the boundary conditions, the leaving particles are sent to the correct neighbouring MPI process !> or deleted. ! !> @param[in] p particles structure ! !> @author Guillaume Le Bars EPFL/SPC !--------------------------------------------------------------------------- SUBROUTINE bound(p) USE basic, ONLY: zgrid, nz, rgrid, nr, Zbounds, mpirank, step, leftproc, rightproc, partperiodic use omp_lib IMPLICIT NONE type(particles), INTENT(INOUT):: p INTEGER :: i,j, rsendnbparts, lsendnbparts, nblostparts INTEGER :: receivednbparts, partdiff LOGICAL:: leftcomm, rightcomm INTEGER, ALLOCATABLE:: partstoremove(:) INTEGER,allocatable :: nblost(:) allocate(nblost(size(p%nblost,1))) nblost=0 IF (p%Nploc .gt. 0) THEN ! We communicate with the left processus leftcomm = leftproc .ne. -1 ! We communicate with the right processus rightcomm = rightproc .ne. -1 ! Boundary condition at z direction !$OMP DO SIMD DO i=1,p%Nploc p%losthole(i)=0 p%sendhole(i)=0 ! If the particle is above or below the simulation domain IF(p%pos(1,i) .gt. rgrid(nr)) THEN p%losthole(i)=4 cycle else if(p%pos(1,i) .lt. rgrid(0))then p%losthole(i)=3 cycle end if ! If the particle is to the right of the local simulation space, it is sent to the right MPI process IF (p%pos(3,i) .ge. zgrid(Zbounds(mpirank+1))) THEN IF(partperiodic) THEN DO WHILE (p%pos(3,i) .GT. zgrid(nz)) p%pos(3,i) = p%pos(3,i) - zgrid(nz) + zgrid(0) END DO END IF !!$OMP CRITICAL (nbparts) IF(rightcomm) THEN p%sendhole(i)=i ELSE IF(.not. partperiodic) THEN p%losthole(i)=2 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%pos(3,i) .lt. zgrid(Zbounds(mpirank))) THEN IF(partperiodic) THEN DO WHILE (p%pos(3,i) .LT. zgrid(0)) p%pos(3,i) = p%pos(3,i) + zgrid(nz) - zgrid(0) END DO END IF !!$OMP CRITICAL (nbparts) IF(leftcomm) THEN ! We send the particle to the left process p%sendhole(i)=-i ELSE IF(.not. partperiodic) THEN ! we destroy the particle p%losthole(i)=1 END IF !!$OMP END CRITICAL (nbparts) END IF END DO !$OMP END DO SIMD END IF !$OMP MASTER receivednbparts=0 nblost=0 j=1 rsendnbparts=0 lsendnbparts=0 Do i=1,p%Nploc if(p%sendhole(i) .eq. 0) cycle p%sendhole(j)=p%sendhole(i) if(p%sendhole(i).gt.0)then rsendnbparts=rsendnbparts+1 else lsendnbparts=lsendnbparts+1 end if j=j+1 end do j=1 nblostparts=0 Do i=1,p%Nploc if(p%losthole(i) .eq. 0) cycle nblost(p%losthole(i))=nblost(p%losthole(i))+1 p%losthole(j)=i j=j+1 nblostparts=nblostparts+1 end do p%nblost=nblost+p%nblost IF(mpisize .gt. 1) THEN ! We send the particles leaving the local simulation space to the closest neighbour CALL particlescommunication(p, lsendnbparts, rsendnbparts, receivednbparts, (/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, p%losthole(i), .false. ) END DO !WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank END IF ! computes if we received less particles than we sent partdiff=max(lsendnbparts+rsendnbparts-receivednbparts,0) IF(nblostparts + partdiff .gt. 0) THEN ALLOCATE(partstoremove(nblostparts+partdiff)) partstoremove(1:partdiff)=abs(p%sendhole(receivednbparts+1:receivednbparts+partdiff)) partstoremove(partdiff+1:partdiff+nblostparts)=abs(p%losthole(1:nblostparts)) call LSDRADIXSORT(partstoremove,nblostparts + partdiff) ! If we received less particles than we sent, or lost particles we fill the remaining holes with the particles from the end of the parts arrays DO i=nblostparts+partdiff,1,-1 CALL move_part(p, p%Nploc, partstoremove(i)) p%partindex(p%Nploc)=-1 p%Nploc = p%Nploc-1 END DO deallocate(partstoremove) END IF !$OMP END MASTER deallocate(nblost) END subroutine bound !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Check if a particle is outside the simulation domain and remove it if needed !> @param[in] p particles structure !--------------------------------------------------------------------------- SUBROUTINE boundary_loss(p) USE basic, ONLY: rgrid, nr Use geometry, ONLY: geom_weight, dom_weight,is_insidegeom Use omp_lib !--------------------------------------------------------------------------- ! add below the usage of module iiee USE iiee IMPLICIT NONE type(particles), INTENT(INOUT):: p INTEGER :: i,j,isup, nblostparts, iend,nbunch INTEGER, DIMENSION(16)::idwall INTEGER :: nblost(size(p%nblost,1)), ii, Nploc_init, Nploc_new logical::inside nblost=0 nbunch=16 IF (p%Nploc .le. 0) return !$OMP DO DO i=1,p%Nploc,nbunch ! Avoid segmentation fault caused by accessing non relevant data iend=min(i+nbunch-1,p%Nploc) p%losthole(i:iend)=0 ! calculate the weight do determine if a particle is inside the simulation domain. do j=i,iend call p_calc_cellindex(p,j) call is_insidegeom(p%pos(1,j), p%pos(3,j),p%cellindex(1,j),p%cellindex(3,j),idwall(j-i+1),inside) if(.not. inside) then ! If the particle is outside of the vacuum region it is deleted. p%losthole(j)=max(idwall(j-i+1),1) end if end do call geom_weight(p%pos(1,i:iend), p%pos(3,i:iend), p%geomweight(:,i:iend)) END DO !$OMP END DO NOWAIT !!$OMP critical (lostparts_red) ! p%nblost=nblost+p%nblost !!$OMP END CRITICAL (lostparts_red) !$OMP BARRIER !$OMP MASTER nblostparts=0 nblost=0 j=1 Do i=1,p%Nploc if(p%losthole(i) .eq. 0) cycle nblost(p%losthole(i)+4)=nblost(p%losthole(i)+4)+1 p%losthole(j)=i j=j+1 nblostparts=nblostparts+1 end do p%nblost=nblost+p%nblost IF(nblostparts.gt.0) THEN IF(p%iiee_id.gt.0) THEN Nploc_init = partslist(p%iiee_id)%Nploc CALL ion_induced(p, p%losthole, partslist(p%iiee_id), nblostparts) Nploc_new = partslist(p%iiee_id)%Nploc if (Nploc_new-Nploc_init .ge. 1) then DO ii =Nploc_init+1,Nploc_new Call p_calc_cellindex(partslist(p%iiee_id),ii) call geom_weight(partslist(p%iiee_id)%pos(1,ii), partslist(p%iiee_id)%pos(3,ii), partslist(p%iiee_id)%geomweight(:,ii)) END DO end if !---------------------------------------------------------- ! CALL ion_induced(p,losthole,partslist(indpelec)) ! here we call our routine to create electrons out of ! eliminated ions. ! need to define in this file: indpelec (need not to since) ! we have the index p%iiee_id !---------------------------------------------------------- END IF DO i=nblostparts,1,-1 CALL delete_part(p,p%losthole(i)) END DO END IF !$OMP END MASTER !$OMP BARRIER END SUBROUTINE boundary_loss !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Computes the radial and axial cell index of the particle i !> @param[in] p particles structure !> @param[in] i index in p of the particle !--------------------------------------------------------------------------- subroutine p_calc_cellindex(p,i) - use basic, only: rgrid,zgrid,invdz,invdr, nnr, nsubr,nsubz, nnz, invdthet,nr + use basic, only: rgrid,zgrid,invdz,invdr, nnr, nsubr,nsubz, nnz, invdthet, nr, nz integer::i,j,k type(particles)::p k=0 p%cellindex(1,i)=0 do j=1,nsubr - IF (p%pos(1,i) .GE. rgrid(k) .AND. p%pos(1,i) .LE. rgrid(k+nnr(j))) THEN + IF (p%pos(1,i) .GE. rgrid(k) .AND. p%pos(1,i) .LT. rgrid(k+nnr(j))) THEN p%cellindex(1,i)=floor((p%pos(1,i)-rgrid(k))*invdr(j))+k exit end if k=k+nnr(j) end do - IF (p%pos(1,i) .GE. rgrid(nr))p%pos(1,i)=nr-1 + if (p%pos(1,i) .GE. rgrid(nr)) p%pos(1,i)=nr-1 k=0 p%cellindex(3,i)=0 do j=1,nsubz - IF (p%pos(3,i) .GE. zgrid(k) .AND. p%pos(3,i) .LE. zgrid(k+nnz(j))) THEN + IF (p%pos(3,i) .GE. zgrid(k) .AND. p%pos(3,i) .LT. zgrid(k+nnz(j))) THEN p%cellindex(3,i)=floor((p%pos(3,i)-zgrid(k))*invdz(j))+k exit end if k=k+nnz(j) end do + if (p%pos(3,i) .GE. zgrid(nz)) p%pos(3,i)=nz-1 p%cellindex(2,i)=floor(p%pos(2,i)*invdthet) !p%zindex(i)=floor((p%Z(i)-zgrid(0))*invdz) end subroutine p_calc_cellindex !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Computes the magnetic field amplitude at each particle position interpolated from the magnetic field at the closeset grid point !> @param[in] p particles structure !--------------------------------------------------------------------------- SUBROUTINE comp_mag_p(p) USE basic, ONLY: zgrid, rgrid, BZ, BR, nz, nr, invdz type(particles), INTENT(INOUT):: p INTEGER :: i Real(kind=db):: WZ,WR INTEGER:: j1,j2,j3,j4 !$OMP DO SIMD DO i=1,p%Nploc WZ=(p%pos(3,i)-zgrid(p%cellindex(3,i)))/(zgrid(p%cellindex(3,i)+1)-zgrid(p%cellindex(3,i))); WR=(p%pos(1,i)-rgrid(p%cellindex(1,i)))/(rgrid(p%cellindex(1,i)+1)-rgrid(p%cellindex(1,i))); J1=(p%cellindex(3,i))*(nr+1) + p%cellindex(1,i)+1 J2=(p%cellindex(3,i))*(nr+1) + p%cellindex(1,i)+2 J3=(p%cellindex(3,i)+1)*(nr+1)+p%cellindex(1,i)+1 J4=(p%cellindex(3,i)+1)*(nr+1)+p%cellindex(1,i)+2 ! Interpolation for magnetic field p%B(2,i)=(1-WZ)*(1-WR)*Bz(J4) & & +WZ*(1-WR)*Bz(J3) & & +(1-WZ)*WR*Bz(J2) & & +WZ*WR*Bz(J1) p%B(1,i)=(1-WZ)*(1-WR)*Br(J4) & & +WZ*(1-WR)*Br(J3) & & +(1-WZ)*WR*Br(J2) & & +WZ*WR*Br(J1) END DO !$OMP END DO SIMD NOWAIT end subroutine comp_mag_p !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the classical simulations. !> This routine systematically returns 1.0 to treat the system according to classical dynamic. ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- REAL(kind=db) FUNCTION gamma_classical(UZ, UR, UTHET) !!#if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_classical) !!#endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_classical=1.0 END FUNCTION gamma_classical !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the relativistic simulations. !> This routine computes the Lorentz factor \f$\gamma=\sqrt{1+\mathbf{\gamma\beta}^2}\f$ ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- REAL(kind=db) FUNCTION gamma_relativistic(UZ, UR, UTHET) !!#if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_relativistic) !!#endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_relativistic=sqrt(1+UZ**2+UR**2+UTHET**2) END FUNCTION gamma_relativistic !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief General routine to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity(p) ! ! 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) !$OMP master CALL swappointer2(p%Uold, p%U) !CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) !$OMP end master !$OMP BARRIER IF (nlclassical) THEN CALL comp_velocity_fun(p, gamma_classical) ELSE CALL comp_velocity_fun(p, gamma_relativistic) END IF END SUBROUTINE comp_velocity !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine called by comp_velocity to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using the routine computing gamma as an input. This avoid the nlclassical flag check on each particle. ! !> @param[in] gamma the function used to compute the value of the lorentz factor \f$\gamma\f$ !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity_fun(p, gammafun) ! ! 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):: ZBR(16), ZBZ(16), ZPR(16), ZPZ(16), ZPTHET(16), SQR(16), ZBZ2(16), ZBR2(16) INTEGER:: i, iend, nbunch, j,k nbunch=16 ! Normalized time increment tau=p%qmRatio*bnorm*0.5*dt*tnorm IF (p%Nploc .NE. 0) THEN !$OMP DO DO i=1,p%Nploc,nbunch iend = min(i + nbunch - 1, p%Nploc) !$DIR SIMD DO j=i,iend k=j-i+1 ! First half of electric pulse p%U(3,j)=p%Uold(3,j)+p%E(2,j)*tau p%U(1,j)=p%Uold(1,j)+p%E(1,j)*tau p%Gamma(j)=gammafun(p%U(3,j), p%U(1,j), p%Uold(2,j)) ! Rotation along magnetic field ZBZ(k)=tau*p%B(2,j)/p%Gamma(j) ZBR(k)=tau*p%B(1,j)/p%Gamma(j) ZPZ(k)=p%U(3,j)-ZBR(k)*p%Uold(2,j) !u'_{z} ZPR(k)=p%U(1,j)+ZBZ(k)*p%Uold(2,j) !u'_{r} ZPTHET(k)=p%Uold(2,j)+(ZBR(k)*p%U(3,j)-ZBZ(k)*p%U(1,j)) !u'_{theta} SQR(k)=1+ZBZ(k)*ZBZ(k)+ZBR(k)*ZBR(k) ZBZ2(k)=2*ZBZ(k)/SQR(k) ZBR2(k)=2*ZBR(k)/SQR(k) p%U(3,j)=p%U(3,j)-ZBR2(k)*ZPTHET(k) !u+_{z} p%U(1,j)=p%U(1,j)+ZBZ2(k)*ZPTHET(k) !u+_{r} p%U(2,j)=p%Uold(2,j)+(ZBR2(k)*ZPZ(k)-ZBZ2(k)*ZPR(k)) !u+_{theta} ! Second half of acceleration p%U(3,j)=p%U(3,j)+p%E(2,j)*tau p%U(1,j)=p%U(1,j)+p%E(1,j)*tau ! Final computation of the Lorentz factor p%Gamma(j)=gammafun(p%U(3,j), p%U(1,j), p%U(2,j)) END DO END DO !$OMP END DO NOWAIT END IF p%collected=.false. END SUBROUTINE comp_velocity_fun !--------------------------------------------------------------------------- !> @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_fun3d(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*0.5*dt*tnorm IF (p%Nploc .NE. 0) THEN !$OMP DO SIMD DO i=1,p%Nploc ! First half of electric pulse p%U(3,i)=p%Uold(3,i)+p%E(3,i)*tau p%U(2,i)=p%Uold(2,i)+p%E(2,i)*tau p%U(1,i)=p%Uold(1,i)+p%E(1,i)*tau p%Gamma(i)=gammafun(p%U(3,i), p%U(1,i), p%Uold(2,i)) ! Rotation along magnetic field ZBZ=tau*p%B(2,i)/p%Gamma(i) ZBR=tau*p%B(1,i)/p%Gamma(i) ZPZ=p%U(3,i)-ZBR*p%Uold(2,i) !u'_{z} ZPR=p%U(1,i)+ZBZ*p%Uold(2,i) !u'_{r} ZPTHET=p%Uold(2,i)+(ZBR*p%U(3,i)-ZBZ*p%U(1,i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR p%U(3,i)=p%U(3,i)-ZBR2*ZPTHET !u+_{z} p%U(1,i)=p%U(1,i)+ZBZ2*ZPTHET !u+_{r} p%U(2,i)=p%Uold(2,i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration p%U(3,i)=p%U(3,i)+p%E(3,i)*tau p%U(2,i)=p%U(2,i)+p%E(2,i)*tau p%U(1,i)=p%U(1,i)+p%E(1,i)*tau ! Final computation of the Lorentz factor p%Gamma(i)=gammafun(p%U(3,i), p%U(1,i), p%U(2,i)) END DO !$OMP END DO SIMD NOWAIT END IF p%collected=.false. END SUBROUTINE comp_velocity_fun3d !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles position at time t+1 !> This routine computes the particles position at time t+1 according to the Bunemann algorithm. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE push(p) Use basic, ONLY: dt type(particles), INTENT(INOUT):: p REAL(kind=db):: XP, YP, COSA, SINA, 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) !$OMP DO SIMD DO i=1,p%Nploc ! Local Cartesian coordinates XP=p%pos(1,i)+dt*p%U(1,i)/p%Gamma(i) IF(XP.le.0)THEN XP=abs(XP) p%U(1,i)=-p%U(1,i) end IF YP=dt*p%U(2,i)/p%Gamma(i) ! Conversion to cylindrical coordiantes p%pos(3,i)=p%pos(3,i)+dt*p%U(3,i)/p%Gamma(i) p%pos(1,i)=sqrt(XP**2+YP**2) ! Computation of the rotation angle IF (p%pos(1,i) .EQ. 0) THEN COSA=1 SINA=0 ALPHA=0 ELSE COSA=XP/p%pos(1,i) SINA=YP/p%pos(1,i) ALPHA=asin(SINA) END IF ! New azimuthal position p%pos(2,i)=MOD(p%pos(2,i)+ALPHA,2*pi) ! Velocity in rotated reference frame U1=COSA*p%U(1,i)+SINA*p%U(2,i) U2=-SINA*p%U(1,i)+COSA*p%U(2,i) p%U(1,i)=U1 p%U(2,i)=U2 END DO !$OMP END DO SIMD NOWAIT END IF !$OMP SINGLE p%collected=.false. !$OMP END SINGLE NOWAIT END SUBROUTINE push !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes several diagnostic quantities !> This routine computes the total kinetic and electric potential energy. !> It keeps track of the reference energy and the number of particle per mpi node. ! !--------------------------------------------------------------------------- SUBROUTINE partdiagnostics ! ! Compute energies ! USE constants, ONLY: vlight USE basic, ONLY: phinorm, cstep, nlclassical, ierr, nbspecies INTEGER:: i,j ! Reset the quantities !$OMP SINGLE ekin=0 epot=0 etot=0 !$OMP END SINGLE NOWAIT ! Computation of the kinetic and potential energy as well as fluid velocities and density !!$OMP PARALLEL DO REDUCTION(+:epot, ekin) DEFAULT(SHARED), PRIVATE(i,j) Do j=1,nbspecies if(.not. partslist(j)%is_field) CYCLE !$OMP DO reduction(+:epot,ekin) DO i=1,partslist(j)%Nploc ! Potential energy epot=epot+(partslist(j)%pot(i)+partslist(j)%potxt(i))*partslist(j)%q*partslist(j)%weight ! Kinetic energy IF(.not. nlclassical) THEN ekin=ekin+(0.5*(partslist(j)%Gammaold(i)+partslist(j)%Gamma(i))-1)*partslist(j)%m*partslist(j)%weight ELSE ekin=ekin+0.5*( partslist(j)%U(1,i)*partslist(j)%Uold(1,i) & & + partslist(j)%U(3,i)*partslist(j)%Uold(3,i) & & + partslist(j)%U(2,i)*partslist(j)%Uold(2,i) )*partslist(j)%m*partslist(j)%weight END IF END DO !$OMP END DO NOWAIT END DO !$OMP BARRIER !$OMP MASTER !!$OMP END PARALLEL DO epot=epot*phinorm*0.5 ekin=ekin*vlight**2 ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ. 1) THEN ! Compute the local total energy loc_etot0 = epot+ekin etot0=0 END IF ! Compute the total energy etot=epot+ekin Energies=(/ekin,epot,etot,loc_etot0/) ! The computed energy is sent to the root process IF(mpisize .gt.1) THEN IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) etot0=etot0+Energies(4) ekin=Energies(1) epot=Energies(2) etot=Energies(3) ELSE CALL MPI_REDUCE(Energies, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) END IF ELSE etot0=etot0+loc_etot0 END IF loc_etot0=0 ! Send the local number of particles on each node to the root process IF(mpisize .gt. 1) THEN Nplocs_all(mpirank)=partslist(1)%Nploc IF(mpirank .eq.0 ) THEN CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) !CALL MPI_REDUCE(MPI_IN_PLACE,partslist(1)%nudcol,3,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) partslist(1)%Nptot=sum(Nplocs_all) !partslist(1)%nudcol=partslist(1)%nudcol/partslist(1)%Nptot ELSE CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) !CALL MPI_REDUCE(partslist(1)%nudcol,partslist(1)%nudcol,3,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) END IF ELSE partslist(1)%Nptot=partslist(1)%Nploc END IF !$OMP END MASTER !$OMP BARRIER end subroutine partdiagnostics !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Collect the particles positions and velocities on the root process. !> If the collection has already been performed at this time step, the routine does nothing. ! !--------------------------------------------------------------------------- SUBROUTINE collectparts(p,p_towrite) USE basic, ONLY: mpirank, mpisize, ierr type(particles), INTENT(INOUT):: p type(particles), INTENT(INOUT), OPTIONAL:: p_towrite INTEGER, DIMENSION(:), ALLOCATABLE :: displs, Nploc INTEGER:: i INTEGER:: particles_type(0:mpisize-1) !< Stores the MPI data type used for particles gathering on node 0 and broadcast from node 0 INTEGER :: part_requests(0:mpisize-1) INTEGER:: stats(MPI_STATUS_SIZE,mpisize) INTEGER:: root_send_particles_type INTEGER:: root_send_part_request part_requests=MPI_REQUEST_NULL root_send_part_request=MPI_REQUEST_NULL particles_type=MPI_DATATYPE_NULL root_send_particles_type=MPI_DATATYPE_NULL IF(p%collected) RETURN ! exit subroutine if particles have already been collected during this time step ALLOCATE(Nploc(0:mpisize-1)) ALLOCATE(displs(0:mpisize-1)) displs=0 Nploc(mpirank)=p%Nploc CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, Nploc, 1, MPI_INTEGER,& & MPI_COMM_WORLD, ierr) p%Nptot=sum(Nploc) IF(p%Nptot .eq. 0 ) THEN p%partindex=-1 p%collected=.true. RETURN END IF Do i=1,mpisize-1 displs(i)=displs(i-1)+Nploc(i-1) END DO IF (present(p_towrite))then IF(mpirank.eq.0 .and. p%Nptot .gt. size(p_towrite%pos,2)) THEN CALL change_parts_allocation(p_towrite,max(p%Nptot-size(p_towrite%pos,2),floor(0.5*size(p_towrite%pos,2)))) END IF else IF(mpirank.eq.0 .and. p%Nptot .gt. size(p%pos,2)) THEN CALL change_parts_allocation(p,max(p%Nptot-size(P%pos,2),floor(0.5*size(P%pos,2)))) END IF end if IF(mpirank .ne. 0) THEN if(Nploc(mpirank) .gt. 0) THEN Call init_particles_gather_mpi(p,1,Nploc(mpirank),particles_type(mpirank)) ! Send Particles informations to root process CALL MPI_SEND(p, 1, particles_type(mpirank), 0, partsgather_tag, MPI_COMM_WORLD, ierr) CALL MPI_TYPE_FREE(particles_type(mpirank),ierr) END IF ELSE ! Receive particle information from all processes if(present(p_towrite))then if(Nploc(mpirank) .gt. 0) THEN Call init_particles_gather_mpi(p,1,Nploc(mpirank),root_send_particles_type) ! Send Particles informations to root process CALL MPI_ISEND(p, 1, root_send_particles_type, 0, partsgather_tag, MPI_COMM_WORLD, root_send_part_request, ierr) end if DO i=0,mpisize-1 if(Nploc(i) .lt. 1) cycle Call init_particles_gather_mpi(p_towrite,displs(i)+1,Nploc(i),particles_type(i)) CALL MPI_IRECV(p_towrite,1,particles_type(i),i,partsgather_tag,MPI_COMM_WORLD, part_requests(i), ierr) END DO CALL MPI_WAITALL(mpisize,part_requests, stats, ierr) if(Nploc(mpirank).gt. 0) CALL MPI_TYPE_FREE(root_send_particles_type,ierr) p_towrite%partindex(sum(Nploc)+1:)=-1 else DO i=1,mpisize-1 if(Nploc(i) .lt. 1) cycle Call init_particles_gather_mpi(p,displs(i)+1,Nploc(i),particles_type(i)) CALL MPI_IRECV(p,1,particles_type(i),i,partsgather_tag,MPI_COMM_WORLD, part_requests(i), ierr) END DO CALL MPI_WAITALL(mpisize,part_requests, stats, ierr) p%partindex(sum(Nploc)+1:)=-1 end if Do i=1,mpisize-1 if(Nploc(i) .lt. 1) cycle CALL MPI_TYPE_FREE(particles_type(i),ierr) END DO END IF p%collected=.TRUE. END SUBROUTINE collectparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Computes the velocities at time t-1/2 delta t to keep the second order precision in time on the velocity. !> This should only be used at particle initialisation time, ot in the case of a restart. ! !--------------------------------------------------------------------------- SUBROUTINE adapt_vinit(p) !! Computes the velocity at time -dt/2 from velocities computed at time 0 ! USE basic, ONLY : bnorm, dt, tnorm, nlclassical, phinorm, distribtype, vnorm type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, Vperp, v2 INTEGER :: i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ! In case Davidson distribution is used the longitudinal and radial velocities are adapted to take into account the ! electric potential. IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4 .or. p%Davidson) THEN ALLOCATE(VR(p%Nploc),VZ(p%Nploc),VTHET(p%Nploc)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=p%U(2,:)/p%Gamma*vnorm DO i=1,p%Nploc Vperp=sqrt(MAX(2*p%H0/p%m-2*p%qmRatio*p%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN p%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 p%Gamma(i)=sqrt(1/(1-v2/vnorm**2)) END IF p%U(1,i)=p%Gamma(i)*VR(i)/vnorm p%U(3,i)=p%Gamma(i)*VZ(i)/vnorm p%U(2,i)=p%Gamma(i)*VTHET(i)/vnorm END DO DEALLOCATE(VR,VZ,VTHET) END IF ! Normalised time increment !tau=-omegac/2/omegap*dt/tnorm tau=-p%qmRatio*bnorm*0.5*dt*tnorm ! Store old Velocities CALL swappointer2(p%Uold, p%U) !CALL swappointer(p%URold, p%UR) !CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR) DO i=1,p%Nploc ! Half inverse Rotation along magnetic field ZBZ=tau*p%B(2,i)/p%Gammaold(i) ZBR=tau*p%B(1,i)/p%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(p%Uold(3,i)-ZBR*p%Uold(2,i))/SQR !u-_{z} ZPR=(p%Uold(1,i)+ZBZ*p%Uold(2,i))/SQR !u-_{r} ZPTHET=p%Uold(2,i)+(ZBR*p%Uold(3,i)-ZBZ*p%Uold(1,i))/SQR !u-_{theta} p%U(3,i)=ZPZ p%U(1,i)=ZPR p%U(2,i)=ZPTHET ! half of decceleration p%U(3,i)=p%U(3,i)+p%E(2,i)*tau p%U(1,i)=p%U(1,i)+p%E(1,i)*tau IF(.not. nlclassical) THEN p%Gamma(i)=sqrt(1+p%U(3,i)**2+p%U(1,i)**2+p%U(2,i)**2) END IF END DO !$OMP END PARALLEL DO SIMD END IF END SUBROUTINE adapt_vinit !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Calculates the number of particles per column of the spatial grid ( at fixed axial cell position) !> This facilitate the computation of the axial grid limits for each MPI worker ! !--------------------------------------------------------------------------- SUBROUTINE calcnbperz(p,nbperz) USE basic, only: nz IMPLICIT NONE type(particles):: p INTEGER, INTENT(INOUT):: nbperz(0:) Integer::i, zindex !$OMP PARALLEL DO DEFAULT(SHARED) reduction(+:nbperz), private(zindex,i) Do i=1,p%Nploc ! we make sure zindex is in [0, nz-1] to avoid segmentation faults zindex=min(nz-1,max(p%cellindex(3,i),0)) nbperz(zindex)=nbperz(zindex)+1 END DO !$OMP END PARALLEL DO END SUBROUTINE calcnbperz !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief In the case of MPI parallelism, computes the axial limits of the local domain. !--------------------------------------------------------------------------- SUBROUTINE calc_Zbounds(p, Zbounds, norder) ! Computes the start and end indices for the Z boundaries on local processus ! Computes the particle indices from initial particle loading vector, that stay in current process USE basic, ONLY: nz, cstep, mpirank, mpisize, step, nbspecies USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER:: Zbounds(0:) INTEGER:: norder(2) INTEGER:: old_Zbounds(0:size(Zbounds,1)-1) INTEGER:: k, i, nbparts REAL(kind=db):: idealnbpartsperproc INTEGER, DIMENSION(0:nz-1):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1) INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction INTEGER:: Zperproc, ierr, remparts CHARACTER(12)::fmt partspercol=0 ! calculatese the axial disstibution integrated along the radial direction do i=1,nbspecies call calcnbperz(partslist(i),partspercol) end do ! gather this data on all nodes if(step .gt. 0 .and. mpisize .gt. 1) THEN old_Zbounds=Zbounds CALL MPI_ALLREDUCE(MPI_IN_PLACE, partspercol, nz, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr) END IF ! estimate the ideal number of particle per MPI worker idealnbpartsperproc = 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 ! Find naive axial limits assuming uniform axial distribution IF(Zmax .le. 0) Zmax=nz-1 IF(Zmin .gt. nz) Zmin=0 Zperproc=(Zmax-Zmin)/mpisize IF (Zperproc .lt. 1) THEN !! No particles are present initially Zperproc=nz/mpisize Zmin=0 ! Define boundaries using naive guess on start or restart (allow to start with 0 parts) DO k=1,mpisize-1 IF(k .lt. mpisize-1-MODULO(Zmax-Zmin,mpisize)) THEN Zbounds(k)=Zmin+k*Zperproc-1 ELSE Zbounds(k)=Zmin+k*Zperproc-1+k-mpisize+2+MODULO(Zmax-Zmin,mpisize) END IF END DO ELSE i=0 ! Define axial boundaries using the axial distribution information. ! the subdomains are not equal remparts=p%Nptot DO k=1,mpisize-1 nbparts=0 DO WHILE(nbparts<0.98*idealnbpartsperproc .and. i .lt. Zmax .and. (nbparts+partspercol(i)).lt.1.05*idealnbpartsperproc) nbparts=nbparts+partspercol(i) i=i+1 END DO remparts=remparts-nbparts Zbounds(k)=i END DO END IF IF(step .gt. 0 .and. mpirank .eq. 0) THEN Do i=1,mpisize-1 !We check that the new limits will not exceed the old limits of the left and right process ! This avoids particle communications with process >mpirank+2 and < mpirank-2 ! However this should converge over time IF(Zbounds(i) .lt. old_Zbounds(i-1)) Zbounds(i)=old_Zbounds(i-1) if(Zbounds(i) .gt. old_Zbounds(i+1))Zbounds(i)=old_Zbounds(i+1) ! If a process would have an axial domain shoter than axial norder, we revert to the old boundaries. IF((Zbounds(i)-Zbounds(i-1)).lt. norder(1) .or. (Zbounds(i+1)-Zbounds(i)).lt. norder(1)) THEN Zbounds=old_Zbounds EXIT END IF END DO END IF ! send the new boundaries to all the workers CALL MPI_Bcast(Zbounds,mpisize+1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) DO k=0,mpisize-1 Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1)) END DO if(mpirank .eq. 0) THEN WRITE(fmt,'(a,i3,a)')"(a,",mpisize+1, "i5)" WRITE(*,fmt) "Zbounds: ", Zbounds WRITE(fmt,'(a,i3,a)')"(a,",mpisize, "i8)" WRITE(*,fmt) "Nplocs: ", Nplocs_all END IF END SUBROUTINE calc_Zbounds !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief After a restart keep only the particles in the local domain of the current MPI worker !--------------------------------------------------------------------------- SUBROUTINE keep_mpi_self_parts(p,Zbounds) TYPE(particles),INTENT(INOUT):: p INTEGER,INTENT(in)::Zbounds(0:) INTEGER :: i, partstart, old_sum,ierr partstart=1 p%Nploc=0 Do i=1,p%Nptot IF(p%cellindex(3,i) .ge.Zbounds(mpirank).and.p%cellindex(3,i).lt.Zbounds(mpirank+1)) THEN p%Nploc=p%Nploc+1 CALL move_part(p,i,p%Nploc) END IF END DO old_sum=p%Nptot CALL MPI_REDUCE(p%Nploc, p%Nptot,1,MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr) IF(p%Nptot .ne. old_sum) THEN WRITE(*,*) "Error in particle distribution kept: ", p%Nptot, "/",old_sum !call MPI_Abort(MPI_COMM_WORLD, -1, ierr) !stop END IF END SUBROUTINE keep_mpi_self_parts !_______________________________________________________________________________ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Manage the particle communication between neighbours. !> This routine is responsible to receive the incoming particles from the MPI neighbours and to send its outgoing !> particles to these neighbours ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !--------------------------------------------------------------------------- SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, receivednbparts, procs) USE mpihelper, ONLY: particle_type #ifdef _DEBUG USE basic, ONLY: step #endif type(particles), INTENT(INOUT):: p INTEGER, INTENT(inout) :: lsendnbparts, rsendnbparts INTEGER, INTENT(out) :: receivednbparts INTEGER, INTENT(in) :: procs(2) INTEGER :: rrecvnbparts=0, lrecvnbparts=0 INTEGER :: sendrequest(2), recvrequest(2) INTEGER :: sendstatus(MPI_STATUS_SIZE,2), recvstatus(MPI_STATUS_SIZE,2) TYPE(particle), ALLOCATABLE :: rrecvpartbuff(:), lrecvpartbuff(:), rsendpartbuff(:), lsendpartbuff(:) ! buffers to send and receive particle from left and right processes INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts, ierr lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL lrecvnbparts=0 rrecvnbparts=0 ! Send and receive the number of particles to exchange 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, lsendpartbuff, rsendpartbuff) END IF CALL MPI_Waitall(2,sendrequest(1:2), sendstatus(:,1:2), ierr) ! Send the particles to the left and right neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, procs(1), partsexchange_tag, MPI_COMM_WORLD, sendrequest(1), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", lsendnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, procs(2), partsexchange_tag, MPI_COMM_WORLD, sendrequest(2), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", rsendnbparts , " to right at step: ",step #endif END IF receivednbparts=rreceivednbparts+lreceivednbparts IF(p%Nploc+receivednbparts-lsendnbparts-rsendnbparts .gt. size(p%pos,2)) THEN CALL change_parts_allocation(p,receivednbparts) END IF ! Receive the incoming parts in the receive buffers IF ( lreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(1), recvstatus(:,1), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,1) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", lreceivednbparts , " from left at step: ",step #endif END IF IF ( rreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(2), recvstatus(:,2), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,2) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", rreceivednbparts , " from right at step: ",step #endif END IF ! Copy the incoming particles from the receive buffers to the simulation parts variable CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, & & lrecvpartbuff, rrecvpartbuff) ! Wait for the outgoing particles to be fully received by the neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(1), sendstatus(:,1), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", lsentnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(2), sendstatus(:,2), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", rsentnbparts , " to right at step: ",step #endif END IF ! ! END SUBROUTINE particlescommunication !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the receive buffers to the local simulation variable parts. !> The incoming particles will first be stored in the holes left by the outgoing particles, then they !> will be added at the end of the parts variable ! !> @param [in] rrecvnbparts number of particles received from the right neighbour (mpirank+1) !> @param [in] lrecvnbparts number of particles received from the left neighbour (mpirank-1) !> @param [in] sendnbparts total number of particles having left the local domain !--------------------------------------------------------------------------- SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts,lrecvpartbuff, rrecvpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts TYPE(particle), INTENT(IN) :: rrecvpartbuff(:), lrecvpartbuff(:) INTEGER k,partpos ! First import the particles coming from the right IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(p%sendhole(k)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, rrecvpartbuff(k), partpos) END DO END IF ! Then import the particles coming from the left IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(p%sendhole(k+rrecvnbparts)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, lrecvpartbuff(k), partpos) END DO END IF ! END SUBROUTINE Addincomingparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the local parts variable to the left and right send buffers. ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !--------------------------------------------------------------------------- SUBROUTINE AddPartSendBuffers(p, lsendnbparts, rsendnbparts, lsendpartbuff, rsendpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts TYPE(particle), INTENT(OUT) :: rsendpartbuff(:), lsendpartbuff(:) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 ! Loop over the outgoing particles and fill the correct send buffer Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(p%sendhole(k)) IF(p%sendhole(k) .GT. 0) THEN rsendpos=rsendpos+1 CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos) ELSE IF(p%sendhole(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(p, lsendpartbuff, lsendpos, partpos) END IF END DO ! ! END SUBROUTINE AddPartSendBuffers !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Add the particles stored in the buffer to the main particle storage p in particles form !> @param[in] p particles structure to add particles to !> @param[in] buffer memory containing the particles to be added !> @param[in] nb_ins number of particles stored in buffer !--------------------------------------------------------------------------- SUBROUTINE add_list_created_part(p, buffer,nb_ins) IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:) INTEGER, OPTIONAL:: nb_ins INTEGER:: i, nptotinit, parts_size_increase, nb_added nptotinit=p%Nploc+1 if(present(nb_ins)) THEN nb_added=nb_ins ELSE nb_added=size(buffer,1) end if IF(nb_added .le. 0) RETURN ! No particles to add ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nploc + nb_added .gt. size(p%pos,2)) THEN parts_size_increase=Max(floor(0.1*size(p%pos,2)),nb_added) CALL change_parts_allocation(p, parts_size_increase) END IF DO i=1,nb_added CALL add_created_particle(p,buffer(i)) END DO nb_added=p%Nploc-nptotinit+1 if(p%is_field) then IF(allocated(p%addedlist)) then call change_array_size_int(p%addedlist,2) else allocate(p%addedlist(2)) end if p%addedlist(size(p%addedlist)-1)=nptotinit p%addedlist(size(p%addedlist))=nb_added end if END SUBROUTINE add_list_created_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Add the particles stored in the linked buffer to the main particle storage p in particles form !> @param[in] p particles structure to add particles to !> @param[in] linked_buffer memory containing the particles to be added in linked list format !> @param[in] destroy Indicates if the memory of the linked buffer must be freed after copy to p !> @param[in] zerovelocity Define if the velocity of the particles in p is set to 0 or copied from the buffer !--------------------------------------------------------------------------- SUBROUTINE add_linked_created_part(p, linked_buffer, destroy, zerovelocity) IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(linked_part_row), INTENT(in) :: linked_buffer LOGICAL:: destroy, zerovelocity TYPE(linked_part), POINTER:: part INTEGER:: i, nptotinit, parts_size_increase, nb_added nptotinit=p%Nploc+1 nb_added=linked_buffer%n IF(nb_added .le. 0) RETURN ! No particles to add ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nploc + nb_added .gt. size(p%pos,2)) THEN parts_size_increase=Max(floor(0.2*size(p%pos,2)),nb_added) CALL change_parts_allocation(p, parts_size_increase) END IF part=>linked_buffer%start DO i=1,nb_added CALL add_created_particle(p,part%p) part=>part%next END DO nb_added=p%Nploc-nptotinit+1 if(p%is_field) then IF(allocated(p%addedlist)) then call change_array_size_int(p%addedlist,2) else allocate(p%addedlist(2)) end if p%addedlist(size(p%addedlist)-1)=nptotinit p%addedlist(size(p%addedlist))=nb_added end if if(zerovelocity)then p%U(:,nptotinit:p%Nploc)=0 !p%UTHET(nptotinit:p%Nploc)=0 !p%UZ(nptotinit:p%Nploc)=0 end if if (destroy) call destroy_linked_parts(linked_buffer%start) if (p%is_field) then ! we keep track of energy by removing the ionisation energy ! with conversion from electronvolt to joules loc_etot0=loc_etot0-sum(p%pot(nptotinit:p%Nploc)*elchar) end if END SUBROUTINE add_linked_created_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Add created particles from a buffer of type particle to the main species storages. ! !> @param [in] p specie memory where we want to add particles !> @param [in] part particle buffer storing the data we want to add to p !--------------------------------------------------------------------------- SUBROUTINE add_created_particle(p,part) USE geometry TYPE(particles):: p TYPE(particle):: part p%Nploc=p%Nploc+1 p%newindex=p%newindex+1 ! add the data to the p structure CALL Insertincomingpart(p, part, p%Nploc) p%partindex(p%Nploc)=p%newindex ! calculate the new domain weight CALL dom_weight(p%pos(1,p%Nploc),p%pos(3,p%Nploc),p%geomweight(0,p%Nploc)) ! delete the particle if it is outside of the computational domain if( .not. is_inside(p,p%Nploc) ) then p%Nploc=p%Nploc-1 p%newindex=p%newindex-1 RETURN end if ! Calculate the geometric weight for the Poisson solver and the grid indices CALL geom_weight(p%pos(1,p%Nploc),p%pos(3,p%Nploc),p%geomweight(:,p%Nploc)) call p_calc_cellindex(p,p%Nploc) END SUBROUTINE add_created_particle !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Checks if the particle id in p is inside of the simulation domain ! !> @param [in] p specie memory !> @param [in] id index of the particle we want to test !--------------------------------------------------------------------------- function is_inside(p,id) Use basic, ONLY: rgrid,zgrid, nr, nz IMPLICIT NONE logical :: is_inside type(particles) :: p integer :: id is_inside=.true. ! Check if the particle is in the simulation domain if(p%geomweight(0,id).le.0)then is_inside=.false. return end if ! check if the particle is in the simulation grid if(p%pos(1,id).ge.rgrid(nr) .or. p%pos(1,id) .le. rgrid(0))then is_inside=.false. return end if if(p%pos(3,id).ge.zgrid(nz) .or. p%pos(3,id) .le. zgrid(0))then is_inside=.false. return end if end function is_inside !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Calculate the energy added by new particles to the system for diagnostic purposes ! !> @param [in] p specie memory !--------------------------------------------------------------------------- SUBROUTINE calc_newparts_energy(p) USE basic, ONLY: phinorm, nlclassical USE omp_lib type(particles)::p integer::i,n,nptotinit,nbadded, nptotend, j Real(kind=db):: addedE ! exit if these particles dont participate in the Poisson equation if(.not. p%is_field) return if( allocated(p%addedlist)) then addedE=0.0_db n=size(p%addedlist) !$OMP SINGLE nbadded=0 Do i=1,n,2 nbadded=nbadded+p%addedlist(i+1) end do p%nbadded=p%nbadded+nbadded !$OMP END SINGLE !$OMP BARRIER ! For each set of added particles Do i=1,n,2 nptotinit=p%addedlist(i) nbadded=p%addedlist(i+1) nptotend=nptotinit+nbadded-1 ! Potential energy\ !$OMP DO DO j=nptotinit,nptotend,1 addedE=addedE+p%q*p%weight*p%pot(j)*phinorm ! Kinetic energy IF(.not. nlclassical) THEN addedE=addedE+p%m*p%weight*vlight**2*(0.5*(p%Gamma(j)+p%Gammaold(j))-1) ELSE addedE=addedE+0.5*p%m*p%weight*vlight**2*(p%U(1,j)*p%Uold(1,j) & & +p%U(3,j)*p%Uold(3,j) & & +p%U(2,j)*p%Uold(2,j)) END IF END DO !$OMP END DO end do !$OMP ATOMIC UPDATE loc_etot0=loc_etot0+addedE !$OMP END ATOMIC !$OMP BARRIER !$OMP MASTER deallocate(p%addedlist) !$OMP END MASTER end if end subroutine calc_newparts_energy !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Delete particle at given index removing its energy from the diagnosed quantities ! !> @param [in] index index of particle to be deleted !--------------------------------------------------------------------------- SUBROUTINE delete_part(p, index, replace) !! This will destroy particle at the given index USE constants, ONLY: vlight USE bsplines USE geometry USE basic, ONLY: phinorm, nlclassical TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index LOGICAL, OPTIONAL :: replace LOGICAL:: repl IF(present(replace)) THEN repl=replace ELSE repl=.true. END IF !Computes the potential at the particle position with phi_ext+phi_s IF(index .le. p%Nploc) THEN IF(p%is_field) THEN loc_etot0=loc_etot0-p%q*p%weight*(p%pot(index))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0-p%m*p%weight*vlight**2*(p%Gamma(index)-1) ELSE loc_etot0=loc_etot0-0.5*p%m*p%weight*vlight**2*(p%U(1,index)**2+p%U(3,index)**2+p%U(2,index)**2) END IF END IF IF(repl) THEN ! We fill the gap CALL move_part(p, p%Nploc, index) p%partindex(p%Nploc)=-1 ! Reduce the total number of simulated parts p%Nploc=p%Nploc-1 END IF END IF END SUBROUTINE delete_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Loads a uniform density of particles on a rectangular annulus qith maxwellian velocities ! !> @param [inout] p particle memory to load into !> @param [inout] VR array of radial velocity for the particles !> @param [inout] VTHET array of azimuthal velocity for the particles !> @param [inout] VZ array of axial velocity for the particles !--------------------------------------------------------------------------- SUBROUTINE loaduniformRZ(p, VR,VZ,VTHET) USE basic, ONLY: plasmadim, rnorm, temp, qsim, msim USE constants, ONLY: me, kb, elchar REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=size(VR,1) p%q=sign(elchar,qsim) p%weight=msim/me p%m=me p%qmRatio=qsim/msim ! Initial distribution in z with normalisation CALL loduni(1,p%pos(3,1:p%Nploc)) p%pos(3,1:p%Nploc)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*p%pos(3,1:p%Nploc))/rnorm ! Initial distribution in r with normalisation CALL lodlinr(2,p%pos(1,1:p%Nploc),plasmadim(3),plasmadim(4)) p%pos(1,1:p%Nploc)=p%pos(1,1:p%Nploc)/rnorm ! Initial velocities distribution CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) END SUBROUTINE loaduniformRZ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Loads a cloud of electrons trapped in a magnetic mirror according to Davidsons equilibrium !> p117 of physics of non-neutral plasma book ! !> @param [inout] p particle memory to load into !> @param [inout] VR array of radial velocity for the particles !> @param [inout] VTHET array of azimuthal velocity for the particles !> @param [inout] VZ array of axial velocity for the particles !--------------------------------------------------------------------------- SUBROUTINE loadDavidson(p, VR,VZ,VTHET, lodr) USE constants, ONLY: me, kb, elchar USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, width, qsim, msim, & & omegac, zgrid, nz, rnorm, n0, nblock, temp procedure(rloader)::lodr TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb, z REAL(kind=db) :: r0, deltar2, halfLz, Mirrorratio, Le, VOL INTEGER :: j, n, blockstart, blockend, addedpart, remainparts INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=p%Nploc Allocate(ra(nblock),rb(nblock), z(0:nblock)) !r0=(plasmadim(4)+plasmadim(3))/2 r0=sqrt(4*H0/(me*omegac**2)) halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) z(0)=plasmadim(1) DO n=1,nblock ! Compute limits in radius and load radii for each part Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1) z(n)=z(0)+n*(plasmadim(2)-plasmadim(1))/nblock deltar2=1-MirrorRatio*cos(2*pi*Le/width) rb(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2)) ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/H0*deltar2)) END DO VOL=SUM(2*pi*MINVAL(ra)*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) qsim=VOL*n0*elchar/nplasma msim=abs(qsim)/elchar*me p%weight=abs(qsim)/elchar p%m=me p%q=sign(elchar,qsim) p%qmRatio=p%q/p%m blockstart=1 blockend=0 ALLOCATE(blocksize(nblock)) WRITE(*,*) "blocksize: ", size(blocksize), nblock DO n=1,nblock blocksize(n)=nplasma/VOL*2*pi*MINVAL(ra)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock END DO remainparts=p%Nploc-SUM(blocksize) addedpart=1 n=nblock/2 j=1 DO WHILE(remainparts .GT. 0) blocksize(n)=blocksize(n)+addedpart remainparts=remainparts-addedpart n=n+j j=-1*(j+SIGN(1,j)) END DO CALL loadPartSlices(p, lodr, ra, rb, z, blocksize) IF(distribtype .eq. 5) THEN CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) VZ=VZ/4 VR=VR*8 VTHET=VTHET*8 ELSE Call loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) END IF END SUBROUTINE loadDavidson !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the velocities for a cloud of electrons trapped in a magnetic mirror according to Davidsons equilibrium !> p117 of physics of non-neutral plasma book. This equilibrium assume mono energy and mono canonical angular momentum ! !> @param [inout] p particle memory to load into !> @param [inout] VR array of radial velocity for the particles !> @param [inout] VTHET array of azimuthal velocity for the particles !> @param [inout] VZ array of axial velocity for the particles !> @param [in] H0 Total energy of each particle !> @param [in] P0 Initial canonical angular momentum of each particle !--------------------------------------------------------------------------- SUBROUTINE loadDavidsonVelocities(p, VR,VZ,VTHET, H0, P0) USE constants, ONLY: me, kb, elchar USE basic, ONLY: rnorm, Rcurv, B0, width, vnorm, zgrid, nz TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), INTENT(IN):: H0, P0 REAL(kind=db) :: athetpos, rg, zg, halfLz, Mirrorratio, Pcomp, Acomp INTEGER :: i MirrorRatio=(Rcurv-1)/(Rcurv+1) halfLz=(zgrid(nz)+zgrid(0))/2 ! Load velocities theta velocity ! Loading of r and z velocity is done in adapt_vinit to have ! access to parts%pot DO i=1,p%Nploc ! Interpolation for Magnetic potential rg=p%pos(1,i)*rnorm zg=(p%pos(3,i)-halfLz)*rnorm Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width)) Pcomp=P0/rg/p%m Acomp=-p%qmRatio*Athetpos VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/p%m)),Pcomp+Acomp) !VTHET(i)=Pcomp+Acomp END DO VTHET=VTHET/vnorm VZ=0._db VR=0._db p%Davidson=.true. p%H0=H0 p%P0=P0 END SUBROUTINE loadDavidsonvelocities !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles velocities according to a maxwellian distribution of temperature temperature [K] ! !> @param [inout] p particle memory to load into !> @param [inout] VR array of radial velocity for the particles !> @param [inout] VTHET array of azimuthal velocity for the particles !> @param [inout] VZ array of axial velocity for the particles !> @param [in] temperature temperature in [k] of the distribution function !--------------------------------------------------------------------------- SUBROUTINE loadGaussianVelocities(p, VR,VZ,VTHET, temperature) USE basic, ONLY: vnorm USE constants, ONLY: kb REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN):: temperature REAL(kind=db):: vth ! Initial velocities distribution vth=sqrt(2.0/3.0*kb*temperature/p%m)/vnorm !thermal velocity CALL lodgaus(3,VZ) CALL lodgaus(5,VR) CALL lodgaus(7,VTHET) VZ=VZ*vth VR=VR*vth VTHET=VTHET*vth p%temperature=temperature p%Davidson=.false. END SUBROUTINE loadGaussianVelocities !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles velocities with a uniform distribution centered in meanv and limited by meanv+spanv and meanv-spanv ! !> @param [inout] p particle memory to load into !> @param [inout] VR array of radial velocity for the particles !> @param [inout] VTHET array of azimuthal velocity for the particles !> @param [inout] VZ array of axial velocity for the particles !> @param [in] meanv mean velocity in each direction [m/s] !> @param [in] spanv extent of the velocity in each direction above and below the mean velocity [m/s] !--------------------------------------------------------------------------- SUBROUTINE loadFlatTopVelocities(p, VR,VZ,VTHET, meanv, spanv) USE basic, ONLY: vnorm USE constants, ONLY: kb REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT):: meanv(3), spanv(3) ! Initial velocities distribution meanv=meanv/vnorm !thermal velocity spanv=spanv/vnorm CALL loduni(3,VZ) CALL loduni(5,VR) CALL loduni(7,VTHET) VR=(VR*2-1)*spanv(1)+meanv(1) VTHET=(VTHET*2-1)*spanv(2)+meanv(2) VZ=(VZ*2-1)*spanv(3)+meanv(3) p%Davidson=.false. END SUBROUTINE loadFlatTopVelocities !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Load slices of particles defined by axial and radial limits ! !> @param [inout] p particle memory to load into !> @param [in] lodr sampling function definig the particle distribution in r !> @param [in] ra lower radial limit of the slice !> @param [in] rb upper radial limit of the slice !> @param [in] z array giving the axial limits of each slice (slice i is betwwen z(i-1) and z(i)) !> @param [in] blocksize array containing the number of particles for each slice !--------------------------------------------------------------------------- SUBROUTINE loadPartslices(p, lodr, ra, rb, z, blocksize) USE 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%pos(3,blockstart:blockend)) p%pos(3,blockstart:blockend)= (z(n-1)+p%pos(3,blockstart:blockend)*(z(n)-z(n-1)))/rnorm CALL lodr(2, p%pos(1,blockstart:blockend), ra(n), rb(n)) p%pos(1,blockstart:blockend)=p%pos(1,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 INTEGER:: iiee_id !< index of species to add particles to for IIEE INTEGER:: neuttype_id !< index of neutral gas producing ions INTEGER:: material_id !< index determining the type of electrode material LOGICAL:: zero_vel !< logical to chose wether or not el. are gen. with non 0 init. vel. REAL(kind=db):: mass=me REAL(kind=db):: charge=-elchar REAL(kind=db):: weight=1.0 REAL(kind=db):: qmratioscale REAL(kind=db):: meanv(3) !< mean velocity in each direction for velocitytype 3 [m/s] REAL(kind=db):: spanv(3) !< pos/neg extent of velocity in each direction for velocitytype 3 [m/s] CHARACTER(len=256) :: header=' ' !< header of csv file section REAL(kind=db):: H0=3.2e-14 !< Total energy [J] REAL(kind=db):: P0=8.66e-25 !< Canonical angular momentum [kg m^2/s] REAL(kind=db):: temperature=10000 !< temperature in kelvins real(kind=db):: n0 !< density factor LOGICAL :: is_test !< Defines if particle are saved on ittracer or not LOGICAL :: is_field !< Defines if particle contributes to Poisson solver LOGICAL :: calc_moments !< Defines if moments matrix must be calculated each it2d CHARACTER(len=16) :: partformat = 'slices' INTEGER:: i, ierr, openerr NAMELIST /partsload/ nblock, mass, charge, weight, npartsalloc, velocitytype, & & radialtype, temperature, H0, P0, is_test, n0, partformat, meanv, spanv, & & calc_moments, qmratioscale, is_field, iiee_id, neuttype_id, material_id, zero_vel ! Set defaults qmratioscale=1.0 weight=1.0 meanv=0 spanv=0 mass=me charge=-elchar calc_moments=.false. is_test=.false. is_field=.true. iiee_id = -1 neuttype_id=1 material_id=1 zero_vel = .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.'slices') THEN IF( nblock .ge. 1) THEN ALLOCATE(z(0:nblock),ra(nblock),rb(nblock), npartsslice(nblock)) DO WHILE(header(1:8) .ne. '//slices') READ(lu_partfile,'(a)') header END DO DO i=1,nblock READ(lu_partfile,*) z(i-1),ra(i),rb(i),npartsslice(i) END DO READ(lu_partfile,*) z(nblock) CALL creat_parts(p,max(npartsalloc,sum(npartsslice))) p%Nploc=sum(npartsslice) p%Nptot=p%Nploc IF( allocated(VR) ) THEN DEALLOCATE(VR,VZ,VTHET) end if if(.not. allocated(VR)) THEN ALLOCATE(VR(p%Nploc)) ALLOCATE(VZ(p%Nploc)) ALLOCATE(VTHET(p%Nploc)) END IF p%m=mass p%q=charge p%weight=weight p%qmRatio=charge/mass*qmratioscale p%is_test=is_test p%is_field=is_field p%calc_moments=calc_moments p%Newindex=sum(npartsslice) p%iiee_id = iiee_id p%neuttype_id = neuttype_id p%material_id = material_id p%zero_vel = zero_vel SELECT CASE(radialtype) CASE(1) ! 1/R distribution in R CALL loadPartslices(p, lodunir, ra, rb, z, npartsslice) CASE(2) ! flat top distribution in R CALL loadPartslices(p, lodlinr, ra, rb, z, npartsslice) CASE(3) ! 1/R^2 distribution in R CALL loadPartslices(p, lodinvr, ra, rb, z, npartsslice) CASE(4) ! gaussian distribution in R CALL loadPartslices(p, lodgausr, ra, rb, z, npartsslice) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT SELECT CASE(velocitytype) CASE(1) ! Gaussian with temperature CALL loadGaussianVelocities(p, VR, VZ, VTHET, temperature) CASE(2) ! Davidson magnetic mirror high wr equilibrium CALL loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) CASE(3) ! flat top velocity CALL loadFlatTopVelocities(p, VR, VZ, VTHET, meanv, spanv) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of velocity distribution:", velocitytype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT END IF END IF ! The plasma cloud is defined as a set individual particles IF( trim(partformat) .eq. 'parts' ) THEN IF( nblock .ge. 1) THEN !Allocate necessary memory CALL creat_parts(p,max(npartsalloc,nblock)) IF( allocated(VR) ) THEN DEALLOCATE(VR,VZ,VTHET) end if if(.not. allocated(VR)) THEN ALLOCATE(VR(nblock)) ALLOCATE(VZ(nblock)) ALLOCATE(VTHET(nblock)) END IF ! Read the particles from the file DO WHILE(header(1:8) .ne. '//parts') READ(lu_partfile,'(a)') header END DO DO i=1,nblock READ(lu_partfile,*) p%pos(1,i),p%pos(2,i),p%pos(3,i), VR(i), VTHET(i), VZ(i) END DO p%Nploc=nblock p%Nptot=p%Nploc p%m=mass p%q=charge p%Newindex=nblock p%weight=weight p%qmRatio=charge/mass*qmratioscale p%is_test=is_test p%is_field=is_field p%calc_moments=calc_moments p%iiee_id = iiee_id p%neuttype_id = neuttype_id p%material_id = material_id p%zero_vel = zero_vel !normalizations p%pos(1,:)=p%pos(1,:)/rnorm p%pos(3,:)=p%pos(3,:)/rnorm !p%z=p%z/rnorm VR=VR/vnorm VTHET=VTHET/vnorm VZ=VZ/vnorm END IF END IF CLOSE(unit=lu_partfile) END SUBROUTINE !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Increase the number of macroparticles by separating each previous macroparticles into !> samplefactor new macroparticles of equally divided weight. The new sub particles are distributed !> uniformly in space to maintain the density and other moments. ! !> @param [in] samplefactor multiplicator of the number of macroparticles. !> @param [in] p particles type to increase. !--------------------------------------------------------------------------- SUBROUTINE upsample(p, samplefactor) USE basic, ONLY : nplasma, dr, dz INTEGER, INTENT(IN) ::samplefactor TYPE(particles), INTENT(INOUT):: p INTEGER:: i, j, currentindex REAL(kind=db), DIMENSION(p%Nploc) :: spreaddir ! random direction for the spread of each initial macro particle REAL(kind=db) :: dir ! Direction in which the particle is moved REAL(kind=db) :: dl ! Particle displacement used for ! Load and scale the direction angle for spreading the new particles CALL loduni(2, spreaddir) spreaddir=spreaddir*2*pi/samplefactor dl=min(minval(dz,1,dz.GT.0),minval(dr,1,dr.GT.0))/10 DO i=1,p%Nploc DO j=1,samplefactor-1 currentindex=p%Nploc+(i-1)*(samplefactor-1)+j CALL move_part(p,i,currentindex) p%partindex(currentindex)=currentindex dir = spreaddir(i)+2*pi*j/samplefactor p%pos(1,currentindex)=p%pos(1,currentindex) + dl*cos(dir) p%pos(3,currentindex)=p%pos(3,currentindex) + dl*sin(dir) END DO p%partindex(i)=i p%pos(1,i)=p%pos(1,i) + dl*cos(spreaddir(i)) p%pos(3,i)=p%pos(3,i) + dl*sin(spreaddir(i)) END DO nplasma=nplasma*samplefactor p%weight=p%weight/samplefactor p%Nploc=p%Nploc*samplefactor p%Nptot=p%Nptot*samplefactor END SUBROUTINE upsample ! Taken from https://rosettacode.org/wiki/Sorting_algorithms/Radix_sort#Fortran ! No Copyright is exerted due to considerable prior art in the Public Domain. ! This Fortran version by Peter Kelly ~ peter.kelly@acm.org ! ! Permission is hereby granted, free of charge, to any person obtaining ! a copy of this software and associated documentation files (the ! "Software"), to deal in the Software without restriction, including ! without limitation the rights to use, copy, modify, merge, publish, ! distribute, sublicense, and/or sell copies of the Software, and to ! permit persons to whom the Software is furnished to do so, subject to ! the following conditions: ! The above copyright notice and this permission notice shall be ! included in all copies or substantial portions of the Software. ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. ! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY ! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, ! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE ! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ! ! Implementation of a classic Radix Sort LSD style :) SUBROUTINE LSDRADIXSORT(A , N) IMPLICIT NONE ! ! Dummy arguments ! INTEGER :: N INTEGER , target, DIMENSION(0:N - 1) :: A ! All arrays based off zero, one day I'll fix it INTENT (IN) N INTENT (INOUT) A ! ! Local variables ! INTEGER , DIMENSION(0:9) :: counts INTEGER :: digitplace INTEGER :: i INTEGER :: j INTEGER :: largestnum INTEGER, DIMENSION(0:N - 1) :: results ! digitplace = 1 ! Count of the keys largestnum = MAXVAL(A) DO WHILE ( (largestnum/digitplace)>0 ) counts = 0 ! Init the count array DO i = 0 , N - 1 , 1 J = (A(i)/digitplace) J = MODULO(j , 10) counts(j) = counts(j) + 1 END DO ! Change count(i) so that count(i) now contains actual position of this digit in result() ! Working similar to the counting sort algorithm DO i = 1 , 9 , 1 counts(i) = counts(i) + counts(i - 1) ! Build up the prefix sum END DO ! DO i = N - 1 , 0 , -1 ! Move from left to right j = (A(i)/digitplace) j = MODULO(j, 10) results(counts(j) - 1) = A(i) ! Need to subtract one as we are zero based but prefix sum is 1 based counts(j) = counts(j) - 1 END DO ! DO i = 0 , N - 1 , 1 ! Copy the semi-sorted data into the input A(i) = results(i) END DO ! digitplace = digitplace*10 END DO ! While loop RETURN END SUBROUTINE LSDRADIXSORT END MODULE beam diff --git a/src/fields_mod.f90 b/src/fields_mod.f90 index 85d2fc7..fd3cf03 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,1575 +1,1577 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: fields ! !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for initializing the magnetic field, solving the Poisson equation and computing the moments of the particles distribution function !------------------------------------------------------------------------------ MODULE fields USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, & & splrz, splrz_ext, nlperiod, phinorm, nlPhis, nrank, mpirank, mpisize, step, it2d, timera, potxt, erxt, ezxt, splrthet,splrthet_ext USE beam, ONLY: partslist USE bsplines USE mumps_bsplines use mpi Use omp_lib Use mpihelper, ONLY: db_type USE particletypes IMPLICIT NONE REAL(kind=db), allocatable, SAVE :: matcoef(:, :), phi_spline(:), vec1(:), vec2(:) REAL(kind=db), allocatable, SAVE :: loc_moments(:, :), loc_rhs(:),omp_loc_rhs(:,:), gradgtilde(:), fverif(:), ppformwork(:,:,:) INTEGER, SAVE:: loc_zspan TYPE(mumps_mat), SAVE :: femat !< Finite Element Method matrix for the full domain TYPE(mumps_mat), SAVE :: reduccedmat !< Finite Element Method matrix in the redduced web-spline sub-space TYPE(mumps_mat), SAVE :: fematmpi !< Finite Element Method matrix prepared for mpi parallelism INTEGER :: nbmoments = 10 !< number of moments to be calculated and stored INTEGER(kind=omp_lock_kind), Allocatable:: mu_lock(:) !< Stores the lock for fields parallelism real(kind=db), allocatable:: reducedsol(:) CONTAINS SUBROUTINE mag_init USE basic, ONLY: nr, nz,lu_in,cstep USE magnet ALLOCATE (Br((nr + 1)*(nz + 1)), Bz((nr + 1)*(nz + 1))) ALLOCATE (Athet((nr + 1)*(nz + 1))) ! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed) ! or load it from magnetfile if present call timera(0, "mag_init") CALL magnet_init(lu_in,cstep) call timera(1, "mag_init") end subroutine mag_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for solving Poisson and computes the magnetic field on the grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_init USE basic, ONLY: pot, nlperiod, nrank, rhs, volume, rgrid, Dimensionality,nthet,thetgrid USE bsplines USE geometry USE mumps_bsplines USE mpihelper INTEGER :: nrz(2), i, d2, k1, n1 select case(Dimensionality) case (1) ! Auxiliary vectors ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1))) DO i=0,nz vec1(i*(nr+1)+1:(i+1)*(nr+1))=rgrid!(0:nz) vec2(i*(nr+1)+1:(i+1)*(nr+1))=zgrid(i) END DO ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder, ngauss, rgrid, zgrid, splrz, nlppform=nlppform, period=nlperiod) ! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential CALL set_spline(femorder, ngauss, rgrid, zgrid, splrz_ext, nlppform=nlppform, period=nlperiod) !Allocate the work buffer to calculate the ppform d2 = splrz%sp2%dim k1 = splrz%sp1%order n1 = splrz%sp1%nints case(2) ALLOCATE(vec1((nthet+1)*(nr+1)),vec2((nr+1)*(nthet+1))) DO i=0,nthet vec1(i*(nr+1)+1:(i+1)*(nr+1))=rgrid!(0:nz) vec2(i*(nr+1)+1:(i+1)*(nr+1))=thetgrid(i) END DO ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder(2:3), ngauss(2:3), rgrid, thetgrid, splrthet, nlppform=nlppform, period=nlperiod(2:3)) ! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential CALL set_spline(femorder(2:3), ngauss(2:3), rgrid, thetgrid, splrthet_ext, nlppform=nlppform, period=nlperiod(2:3)) !Allocate the work buffer to calculate the ppform d2 = splrz%sp2%dim k1 = splrz%sp1%order n1 = splrz%sp1%nints end select ALLOCATE(ppformwork(d2,k1,n1)) ! Calculate dimension of splines nrz(1) = nz nrz(2) = nr CALL get_dim(splrz, nrank, nrz, femorder) ! Allocate necessary variables ALLOCATE (matcoef(nrank(1), nrank(2))) ALLOCATE (pot((nr + 1)*(nz + 1))) ALLOCATE (potxt((nr + 1)*(nz + 1))) ALLOCATE (Erxt((nr + 1)*(nz + 1))) ALLOCATE (Ezxt((nr + 1)*(nz + 1))) ALLOCATE (rhs(nrank(1)*nrank(2))) ALLOCATE (gradgtilde(nrank(1)*nrank(2))) gradgtilde = 0 ALLOCATE (phi_spline(nrank(1)*nrank(2))) ALLOCATE (volume(nrank(1)*nrank(2))) volume = 0 ALLOCATE (Er((nr + 1)*(nz + 1)), Ez((nr + 1)*(nz + 1))) ALLOCATE (mu_lock(nrank(1)*nrank(2))) do i = 1, nrank(1)*nrank(2) call omp_init_lock(mu_lock(i)) end do end SUBROUTINE fields_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the geometry definition and read it from the standard input !> Precomputes the LHS matrix to solve Poisson abd the RHS effect of the dirichlet boundaries ! !--------------------------------------------------------------------------- SUBROUTINE fields_start USE geometry USE basic, ONLY: nrank implicit none INTEGER:: i,j, k, ierr DOUBLE PRECISION:: val ! set up the geometry module for setting up non-conforming boundary conditions call timera(0, "geom_init") call geom_init(splrz, vec1, vec2) call timera(1, "geom_init") ! Initialisation of FEM matrix CALL init(nrank(1)*nrank(2), 2, femat) ! Calculate and factorise FEM matrix (depends only on mesh) CALL fematrixrz(femat) If (walltype .lt. 0) then allocate (fverif(nrank(1)*nrank(2))) fverif = 0 end if ! Compute the volume of the splines and gtilde for solving E using web-splines CALL comp_volume !$OMP PARALLEL Call comp_gradgtilde !$OMP END PARALLEL if (nlweb) then ! Calculate reduced matrix for use of web splines call timera(0, "reduce femat") call Reducematrix(femat, reduccedmat) call timera(1, "reduce femat") call factor(reduccedmat) !call init(reduccedmat%rank,2,fematmpi,comm_in=MPI_COMM_WORLD) !do i=1,reduccedmat%rank ! do k=reduccedmat%irow(i),reduccedmat%irow(i+1)-1 ! j=reduccedmat%cols(k) ! call putele(fematmpi, i, j, reduccedmat%val(k)) ! end do !end do else call factor(femat) !call init(femat%rank,2,fematmpi,comm_in=MPI_COMM_WORLD) !do i=1,femat%rank ! do k=femat%irow(i),femat%irow(i+1)-1 ! j=femat%cols(k) ! call putele(fematmpi, i, j, femat%val(k)) ! end do !end do end if !call factor(fematmpi) !WRITE(*,*) "fematmpi is factorised" !WRITE(*,*) "Copy and to_mat worked" !CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) call vacuum_field END SUBROUTINE fields_start !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Recomputes the vacuum electric field ! !--------------------------------------------------------------------------- subroutine vacuum_field Use geometry USE basic, ONLY: pot, rhs implicit none INTEGER:: i, iend ! Computes the externally imposed electric field !$OMP PARALLEL private(i) !$OMP DO SIMD do i=1,nrank(1)*nrank(2) rhs(i)=-gradgtilde(i) !rhs = -gradgtilde if (walltype .lt. 0) rhs(i) = rhs(i) + fverif(i) end do !$OMP END DO SIMD !$OMP END PARALLEL call poisson(splrz_ext,reducedsol) call poisson_com(splrz_ext,reducedsol) !$OMP PARALLEL private(i,iend) call Update_phi(splrz_ext) !$OMP BARRIER !$OMP DO ! On the root process, compute the electric field for diagnostic purposes DO i=1,size(pot),16 iend=min(size(pot),i+15) potxt(i:iend) = pot(i:iend) erxt(i:iend) = Er(i:iend) Ezxt(i:iend) = Ez(i:iend) END DO !$OMP END DO NOWAIT !$OMP END PARALLEL end subroutine !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for the communication of moments and rhs grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_comm_init(Zbounds) USE basic, ONLY: nrank USE mpihelper USE omp_lib INTEGER:: Zbounds(0:), nbthreads nbthreads = omp_get_max_threads() - loc_zspan = Zbounds(mpirank + 1) - Zbounds(mpirank) + femorder(1) + loc_zspan = Zbounds(mpirank + 1) - Zbounds(mpirank) + femorder(2) if (allocated(loc_moments)) deallocate (loc_moments) ALLOCATE (loc_moments(nbmoments, loc_zspan*nrank(1))) if (allocated(loc_rhs)) deallocate (loc_rhs) ALLOCATE (loc_rhs(loc_zspan*nrank(1))) if (allocated(omp_loc_rhs)) deallocate (omp_loc_rhs) ALLOCATE (omp_loc_rhs(nbthreads,loc_zspan*nrank(1))) IF (mpisize .gt. 1) THEN CALL init_overlaps(nrank, femorder, Zbounds(mpirank), Zbounds(mpirank + 1), nbmoments) END IF END SUBROUTINE fields_comm_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Construct the right hand side vector used in the FEM Poisson solver ! !> @param[in] plist list of the particles type storing the desired specie parameters ! !--------------------------------------------------------------------------- SUBROUTINE rhscon(plist) USE bsplines use mpi USE basic, ONLY: rhs, Zbounds USE beam, ONLY: particles USE mpihelper Use geometry Use omp_lib type(particles), INTENT(INOUT):: plist(:) INTEGER:: i,j,k IF (nlphis) then ! We calculate the self-consistent field !$OMP DO SIMD Do i=1,size(loc_rhs) loc_rhs(i)=0 end do !$OMP END DO SIMD ! Assemble rhs for each specie Do i = 1, size(plist, 1) if (plist(i)%is_field) CALL deposit_charge(plist(i), loc_rhs) END Do !!$OMP BARRIER !Communicate the overlaps if(mpisize .gt. 1) call rhs_overlap ! Add gradgtilde !$OMP DO SIMD Do i=0,size(loc_rhs)-1 j=i/nrank(1)+Zbounds(mpirank) k=mod(i,nrank(1))+1 loc_rhs(i+1)=loc_rhs(i+1)-gradgtilde(k+j*nrank(1)) end do !$OMP END DO SIMD !add the fverif source for test cases if (walltype .lt. 0)then !$OMP DO Do i=0,size(loc_rhs)-1 j=i/nrank(1)+Zbounds(mpirank) k=mod(i,nrank(1))+1 loc_rhs(i+1)=loc_rhs(i+1)+fverif(k+j*nrank(1)) end do !$OMP END DO end if ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL rhs_gather(rhs) ELSE !in serial copy loc_rhs to rhs !$OMP DO Do i=1,size(loc_rhs) rhs(i)=loc_rhs(i) end do !$OMP END DO END IF ELSE ! We only consider the externally imposed field !$OMP DO Do i=1,size(rhs) rhs(i)=-gradgtilde(i) end do !$OMP END DO END IF END SUBROUTINE rhscon !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculate the 0th 1st and 2nd order moments of the particle p and stores it in moment ! !> @param[in] p the particles type storing the desired specie parameters !> @param[out] moment the 2d array storing the calculated moments ! !--------------------------------------------------------------------------- SUBROUTINE momentsdiag(p) USE bsplines use mpi USE beam, ONLY: particles USE mpihelper Use geometry type(particles), INTENT(INOUT):: p !REAL(kind=db), INTENT(INOUT):: moment(:, :) !$OMP SINGLE loc_moments = 0 ! Reset the moments matrix ! Assemble rhs !$OMP END SINGLE IF (p%Nploc .ne. 0) THEN CALL deposit_moments(p, loc_moments) END IF !$OMP SINGLE if(.not. allocated(p%moments))THEN if(mpirank.eq.0)THEN Allocate(p%moments(nbmoments,nrank(1)*nrank(2))) else Allocate(p%moments(0,0)) end if end if !$OMP END SINGLE ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL moments_gather(p%moments) ELSE !$OMP SINGLE p%moments = loc_moments !$OMP END SINGLE NOWAIT END IF END SUBROUTINE momentsdiag !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles moments (n,v,v^2) from p on the grid ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] p_loc_moments local tensor used to store the moments of the given specie !--------------------------------------------------------------------------- SUBROUTINE deposit_moments(p, p_loc_moments) USE bsplines use mpi USE basic, ONLY: Zbounds USE beam, ONLY: particles USE mpihelper USE geometry USE omp_lib TYPE(particles), INTENT(IN):: p REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: p_loc_moments REAL(kind=db), DIMENSION(:, :), Allocatable:: omp_loc_moments INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft REAL(kind=db) :: vr, vthet, vz, coeff REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) INTEGER:: num_threads num_threads = omp_get_max_threads() nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = min(nbunch, 64) ! Particle bunch size used when calling basfun ! Assemble rhs IF (p%Nploc .gt. 0) THEN !!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(zleft,rleft,jw,it,iend,irow,jcol,mu,k,vr,vz,vthet,coeff,fun,fun2) ALLOCATE (zleft(nbunch), rleft(nbunch)) ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node !allocate(omp_loc_moments(size(p_loc_moments,1),size(p_loc_moments,2))) !omp_loc_moments=0 !$OMP DO DO i = 1, p%Nploc, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, p%Nploc) k = iend - i + 1 ! Localize the particle !CALL locintv(splrz%sp2, p%R(i:iend), rleft(1:k)) !CALL locintv(splrz%sp1, p%Z(i:iend), zleft(1:k)) rleft(1:k) = p%cellindex(1,i:iend) zleft(1:k) = p%cellindex(3,i:iend) ! Compute the value of the splines at the particles positions CALL basfun(p%pos(1,i:iend), splrz%sp1, fun(:, :, 1:k), rleft(1:k) + 1) CALL basfun(p%pos(3,i:iend), splrz%sp2, fun2(:, :, 1:k), zleft(1:k) + 1) DO k = 1, (iend - i + 1) DO jw = 1, (femorder(2) + 1) jcol = zleft(k) + jw - Zbounds(mpirank) DO it = 1, (femorder(1) + 1) irow = rleft(k) + it mu = irow + (jcol - 1)*(nrank(1)) coeff = p%weight*fun(it, 0, k)*fun2(jw, 0, k) ! Add contribution of particle nbunch to rhs grid point mu vr = 0.5*(p%U(1,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(1,i + k - 1)/p%Gammaold(i + k - 1)) vz = 0.5*(p%U(3,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(3,i + k - 1)/p%Gammaold(i + k - 1)) vthet = 0.5*(p%U(2,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(2,i + k - 1)/p%Gammaold(i + k - 1)) call omp_set_lock(mu_lock(mu)) p_loc_moments(1:10,mu)=p_loc_moments(1:10,mu)+coeff*(/1.0_db,vr,vthet,vz, vr*vr, vr*vthet, vr*vz, vthet**2, vthet*vz, vz**2/) call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !!$OMP END PARALLEL DO !$OMP END DO NOWAIT !Do i=1,size(p_loc_moments,2) ! call omp_set_lock(mu_lock(i)) ! p_loc_moments(:,i)=p_loc_moments(:,i)+omp_loc_moments(:,i) ! call omp_unset_lock(mu_lock(i)) !end do !!$OMP END CRITICAL(loc_moments_reduce) DEALLOCATE (fun, fun2, zleft, rleft) END IF END subroutine deposit_moments !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles charges (q) from p on the grid ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] p_loc_moments local tensor used to store the moments of the given specie !--------------------------------------------------------------------------- SUBROUTINE deposit_charge(p, p_loc_moments) USE bsplines use mpi USE constants USE basic, ONLY: Zbounds, rnorm, phinorm USE beam, ONLY: particles USE mpihelper USE geometry USE omp_lib TYPE(particles), INTENT(IN):: p REAL(kind=db), DIMENSION(:), INTENT(INOUT):: p_loc_moments INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) INTEGER:: num_threads, curr_thread real(kind=db):: contrib, chargecoeff num_threads = omp_get_max_threads() nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = min(nbunch, 16) ! Particle bunch size used when calling basfun chargecoeff = p%weight*p%q/(2*pi*eps_0*phinorm*rnorm) ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (p%Nploc .gt. 0) THEN !!!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(i,zleft, rleft, jw, it, iend, irow, jcol, mu, k, fun, fun2, contrib) ALLOCATE (zleft(nbunch), rleft(nbunch)) ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node !allocate(omp_loc_moments(size(p_loc_moments))) !omp_loc_moments=0 zleft=0 rleft=0 curr_thread=omp_get_thread_num()+1 !$OMP DO SIMD Do i=1,size(loc_rhs) omp_loc_rhs(:,i)=0 end do !$OMP END DO SIMD !$OMP DO DO i = 1, p%Nploc, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, p%Nploc) k = iend - i + 1 ! Localize the particle rleft(1:k) = p%cellindex(1,i:iend) zleft(1:k) = p%cellindex(3,i:iend) ! Compute the value of the splines at the particles positions CALL basfun(p%pos(1,i:iend), splrz%sp1, fun, rleft(1:k) + 1) CALL basfun(p%pos(3,i:iend), splrz%sp2, fun2, zleft(1:k) + 1) !CALL geom_weight(p%Z(i:iend),p%R(i:iend),wgeom) DO k = 1, (iend - i + 1) DO jw = 1, (femorder(2) + 1) jcol = zleft(k) + jw - Zbounds(mpirank) DO it = 1, (femorder(1) + 1) irow = rleft(k) + it mu = irow + (jcol - 1)*(nrank(1)) ! Add contribution of particle k to rhs grid point mu !contrib = fun(it, 0, k)*fun2(jw, 0, k)*p%geomweight(0,i + k - 1)*chargecoeff !!$OMP ATOMIC UPDATE - omp_loc_rhs(curr_thread,mu) = omp_loc_rhs(curr_thread,mu) + fun(it, 0, k)*fun2(jw, 0, k)*p%geomweight(0,i + k - 1)*chargecoeff - !!$OMP END ATOMIC + contrib=fun(it, 0, k) + contrib=contrib*fun2(jw, 0, k) + contrib=contrib*p%geomweight(0,i + k - 1)*chargecoeff + omp_loc_rhs(curr_thread,mu) = omp_loc_rhs(curr_thread,mu) + contrib !!$OMP END ATOMIC END DO END DO END DO END DO !$OMP END DO - DEALLOCATE (fun, fun2, zleft, rleft) !$OMP DO SIMD Do i=1,size(p_loc_moments) Do k=1,num_threads p_loc_moments(i)=p_loc_moments(i)+omp_loc_rhs(k,i) end do end do !$OMP END DO SIMD + DEALLOCATE (fun, fun2, zleft, rleft) END IF END subroutine deposit_charge !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers for the overlap grid points !> ! !--------------------------------------------------------------------------- SUBROUTINE rhs_overlap USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc INTEGER:: ierr, i, j !$OMP MASTER !WRITE(*,*) mpirank, "wE communicate overlap rhs" CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs, nrank, femorder, loc_zspan - femorder(2)) !$OMP END MASTER !$OMP BARRIER IF (mpirank .gt. 0) THEN !$OMP DO SIMD DO i = 1, femorder(2)*nrank(1) loc_rhs(i) = loc_rhs(i)& & + rhsoverlap_buffer(i) END DO !$OMP END DO SIMD END IF !!$OMP BARRIER END SUBROUTINE rhs_overlap !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers to reduce the result on the host ! !--------------------------------------------------------------------------- SUBROUTINE rhs_gather(rhs) USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:), INTENT(INOUT):: rhs INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) INTEGER:: overlap_type INTEGER:: rcvoverlap_type displs = Zbounds(0:mpisize - 1)*nrank(1) counts = (Zbounds(1:mpisize) - Zbounds(0:mpisize - 1))*nrank(1) counts(mpisize) = counts(mpisize) + femorder(1)*nrank(1) ! Set communication vector type overlap_type = rhsoverlap_type rcvoverlap_type = rcvrhsoverlap_type !$OMP MASTER IF (mpirank .eq. 0) THEN rhs = 0 END IF CALL MPI_GATHERV(loc_rhs, counts(mpirank + 1), db_type, & & rhs, counts, displs, db_type, 0, MPI_COMM_WORLD, ierr) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE rhs_gather !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers for the overlap grid points !> and reduce the result on the host ! !--------------------------------------------------------------------------- SUBROUTINE moments_gather(moment) USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: moment INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) displs = Zbounds(0:mpisize - 1)*nrank(1)*10 counts = (Zbounds(1:mpisize) - Zbounds(0:mpisize - 1))*10*nrank(1) counts(mpisize) = counts(mpisize) + femorder(2)*10*nrank(1) !$OMP MASTER CALL momentsoverlapcomm(mpirank, leftproc, rightproc, loc_moments, nrank, femorder, loc_zspan - femorder(2)) !$OMP END MASTER !$OMP BARRIER IF (mpirank .gt. 0) THEN !!$OMP PARALLEL DO SIMD DEFAULT(SHARED) private(i) !$OMP DO SIMD DO i = 1, nrank(1)*femorder(2) loc_moments(1:nbmoments, i) = loc_moments(1:nbmoments, i)& & + momentsoverlap_buffer((i-1)*nbmoments+1:i*nbmoments) END DO !$OMP END DO SIMD END IF !$OMP MASTER ! Set communication vector type IF (mpirank .eq. 0) THEN moment = 0 END IF CALL MPI_GATHERV(loc_moments, counts(mpirank + 1), db_type, & & moment, counts, displs, db_type, 0, MPI_COMM_WORLD, ierr) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE moments_gather !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Solves Poisson equation using FEM. Distributes the result on all MPI workers. ! !--------------------------------------------------------------------------- SUBROUTINE poisson(splinevar,reducedsol) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend real(kind=db), allocatable::reducedrhs(:) real(kind=db), allocatable,intent(inout):: reducedsol(:) allocate (reducedrhs(nrank(1)*nrank(2))) allocate (reducedsol(nbreducedspline)) !reduccedmat%mumps_par%ICNTL(11)=1 if (nlweb .and. mpirank.eq.0) then ! we use the web-spline reduction for stability reducedrhs = vmx(etilde, rhs) Call bsolve(reduccedmat, reducedrhs(1:nbreducedspline), reducedsol) else if(mpirank.eq.0) then CALL bsolve(femat, rhs, phi_spline) end if END SUBROUTINE poisson !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Solves Poisson equation using FEM. Distributes the result on all MPI workers. ! !--------------------------------------------------------------------------- SUBROUTINE poisson_com(splinevar,reducedsol) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend real(kind=db), allocatable:: reducedsol(:), tempcol(:) allocate (tempcol(nrank(1)*nrank(2))) !$OMP MASTER if (nlweb) then ! we use the web-spline reduction for stability CALL MPI_Bcast(reducedsol, nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr) else CALL MPI_Bcast(phi_spline, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr) end if if (nlweb) then ! we use the web-spline reduction for stability tempcol = 0 tempcol(1:nbreducedspline) = reducedsol phi_spline = vmx(etildet, tempcol) end if deallocate(reducedsol) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE poisson_com SUBROUTINE poisson_mpi(splinevar) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend real(kind=db), allocatable::reducedrhs(:) real(kind=db), allocatable:: reducedsol(:), tempcol(:) allocate (reducedrhs(nrank(1)*nrank(2))) allocate (reducedsol(nbreducedspline)) allocate (tempcol(nrank(1)*nrank(2))) if (nlweb) then ! we use the web-spline reduction for stability reducedrhs = vmx(etilde, rhs) Call bsolve(fematmpi, reducedrhs(1:nbreducedspline), reducedsol) tempcol = 0 tempcol(1:nbreducedspline) = reducedsol !phi_spline = 0 phi_spline = vmx(etildet, tempcol) else CALL bsolve(fematmpi, rhs, phi_spline) !CALL MPI_Bcast(phi_spline, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr) end if END SUBROUTINE poisson_mpi !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Updates the splinevar variable with the new phi coefficients and calculates !> Phi Er and Ez on the grid ! !--------------------------------------------------------------------------- SUBROUTINE Update_phi(splinevar) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr, i, j, iend !$OMP DO SIMD collapse(2) Do j=1,nrank(2) Do i=1,nrank(1) matcoef(i,j) = phi_spline((j-1)*nrank(1)+i) END DO END DO !$OMP END DO SIMD ! update the ppform coefficients CALL updt_ppform2d(splinevar, matcoef) !$OMP BARRIER IF (mpirank .eq. 0 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN !$OMP DO ! On the root process, compute the electric field for diagnostic purposes DO i=1,size(pot),16 iend=min(size(pot),i+15) CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), pot(i:iend), (/0, 0/)) CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Er(i:iend), (/1, 0/)) CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Ez(i:iend), (/0, 1/)) Ez(i:iend) = -pot(i:iend)*gridwdir(1,i:iend) - Ez(i:iend)*gridwdir(0,i:iend) - gtilde(1,i:iend) Er(i:iend) = -pot(i:iend)*gridwdir(2,i:iend) - Er(i:iend)*gridwdir(0,i:iend) - gtilde(2,i:iend) pot(i:iend) = pot(i:iend)*gridwdir(0,i:iend) + gtilde(0,i:iend) END DO !$OMP END DO NOWAIT END IF END SUBROUTINE Update_phi !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the electric fields and potential at the particles position for particles !> between positions nstart and nend in the list ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] nstart starting index for the particle list !> @param[in] nend ending index for the particle list !--------------------------------------------------------------------------- SUBROUTINE EFieldscompatparts(p, nstart, nend) Use beam, ONLY: particles Use geometry Use splinebound TYPE(particles), INTENT(INOUT):: p INTEGER, OPTIONAL::nstart, nend INTEGER:: i, iend, nst, nnd INTEGER:: nbunch INTEGER:: num_threads Real(kind=db), ALLOCATABLE:: erext(:), ezext(:), gtildeloc(:, :) if (.not. present(nstart)) nst = 1 if (.not. present(nend)) nnd = p%Nploc !num_threads = omp_get_max_threads() !nbunch = (nnd - nst + 1)/num_threads ! Particle bunch size used when calling basfun !nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = 64 ! Particle bunch size used when calling basfun Allocate (erext(nbunch), ezext(nbunch), gtildeloc(0:2,0:nbunch - 1)) ! Evaluate the electric potential and field at the particles position !$OMP DO SIMD DO i = nst, nnd, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, nnd) CALL speval(splrz, p%pos(1,i:iend), p%pos(3,i:iend),p%cellindex(1,i:iend),p%cellindex(3,i:iend), p%pot(i:iend), p%E(1,i:iend), p%E(2,i:iend)) CALL speval(splrz_ext, p%pos(1,i:iend), p%pos(3,i:iend),p%cellindex(1,i:iend),p%cellindex(3,i:iend), p%potxt(i:iend)) Call total_gtilde(p%pos(1,i:iend), p%pos(3,i:iend), gtildeloc(:,0:iend - i),p%geomweight(:,i:iend)) p%E(2,i:iend) = -p%E(2,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(1,i:iend) - gtildeloc(1,0:iend - i) p%E(1,i:iend) = -p%E(1,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(2,i:iend) - gtildeloc(2,0:iend - i) p%pot(i:iend) = p%geomweight(0,i:iend)*p%pot(i:iend) + gtildeloc(0,0:iend - i) p%potxt(i:iend) = p%geomweight(0,i:iend)*p%potxt(i:iend) + gtildeloc(0,0:iend - i) END DO !$OMP END DO SIMD NOWAIT END SUBROUTINE EFieldscompatparts !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrixrz(mat) USE bsplines USE geometry USE omp_lib USE sparse type(mumps_mat):: mat REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) REAL(kind=db) :: contrib INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :) integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize kterms=8 If (allocated(fun)) deallocate (fun) If (allocated(fun2)) deallocate (fun2) ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2))) If (allocated(wgeom)) deallocate (wgeom) ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2)) ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) !Pointers on the order of derivatives call timera(0, "fematrix") ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO CALL coefeq(splrz%sp1%knots(0:1), idert, iderw, iderg, coefs, kterms) ! Assemble FEM matrix !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, fun, fun2,iid,jid), collapse(2) DO j = 1, nz ! Loop on z position DO i = 1, nr ! Loop on r position !! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) iid=i jid=j if (gausssize .gt. 1) then !If (allocated(wgeom)) deallocate (wgeom) !ALLOCATE (wgeom(0:2,gausssize)) CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) End if DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) DO iw = 1, (1 + femorder(1))*(femorder(2) + 1) irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1 mu2 = irow2 + (jcol2 - 1)*nrank(1) contrib=0.0_db DO igauss = 1, gausssize ! Loop on gaussian weights and positions DO iterm = 1, kterms ! Loop on the two integration dimensions contrib = contrib+wgeom(iderg(iterm, 1),igauss)*wgeom(iderg(iterm, 2),igauss)* & & fun(f(jt, 1), idert(iterm, 1),igauss)*fun(f(iw, 1), idert(iterm, 2),igauss)* & & fun2(f(jt, 2), iderw(iterm, 1),igauss)*fun2(f(iw, 2), iderw(iterm, 2),igauss)* & & wgauss(igauss)*xgauss(igauss, 1) END DO end do call omp_set_lock(mu_lock(mu)) CALL updt_sploc(mat%mat%row(mu), mu2, contrib) call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !$OMP End parallel do DEALLOCATE (f, aux) DEALLOCATE (idert, iderw, coefs, fun, fun2) call timera(1, "fematrix") END SUBROUTINE fematrixrz !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrixrtheta(mat) USE bsplines USE geometry USE omp_lib USE sparse type(mumps_mat):: mat REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) REAL(kind=db) :: contrib INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :) integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize kterms=8 If (allocated(fun)) deallocate (fun) If (allocated(fun2)) deallocate (fun2) ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2))) If (allocated(wgeom)) deallocate (wgeom) ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2)) ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) !Pointers on the order of derivatives call timera(0, "fematrix") ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO CALL coefeq(splrz%sp2%knots(0:1), idert, iderw, iderg, coefs, kterms) ! Assemble FEM matrix !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, fun, fun2,iid,jid), collapse(2) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position !! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) iid=i jid=j if (gausssize .gt. 1) then !If (allocated(wgeom)) deallocate (wgeom) !ALLOCATE (wgeom(0:2,gausssize)) CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) End if DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) DO iw = 1, (1 + femorder(1))*(femorder(2) + 1) irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1 mu2 = irow2 + (jcol2 - 1)*nrank(1) contrib=0.0_db DO igauss = 1, gausssize ! Loop on gaussian weights and positions DO iterm = 1, kterms ! Loop on the two integration dimensions contrib = contrib+wgeom(iderg(iterm, 1),igauss)*wgeom(iderg(iterm, 2),igauss)* & & fun(f(jt, 1), idert(iterm, 1),igauss)*fun(f(iw, 1), idert(iterm, 2),igauss)* & & fun2(f(jt, 2), iderw(iterm, 1),igauss)*fun2(f(iw, 2), iderw(iterm, 2),igauss)* & & wgauss(igauss)*xgauss(igauss, 1) END DO end do call omp_set_lock(mu_lock(mu)) CALL updt_sploc(mat%mat%row(mu), mu2, contrib) call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !$OMP End parallel do DEALLOCATE (f, aux) DEALLOCATE (idert, iderw, coefs, fun, fun2) call timera(1, "fematrix") END SUBROUTINE fematrixrtheta !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the volume of the splines cells needed to display the density in post-processing !--------------------------------------------------------------------------- SUBROUTINE comp_volume USE bsplines USE geometry USE basic, ONLY: Volume REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :), ftestpt(:, :) Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms Real(kind=db)::newcontrib call timera(0, "comp_volume") ALLOCATE (fun(1:femorder(1) + 1, 0:1), fun2(1:femorder(2) + 1, 0:1))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines nterms = 4 Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms)) ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO volume = 0 if (walltype .lt. 0) fverif = 0 ! Assemble Volume matrix !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib), collapse(2) DO j = 1, nz ! Loop on z position DO i = 1, nr ! Loop on r position ! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) If (allocated(wgeom)) deallocate (wgeom) if (gausssize .gt. 0) then ALLOCATE (wgeom(0:2,size(xgauss, 1))) CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom) End if if (walltype .lt. 0) then If (allocated(ftestpt)) deallocate (ftestpt) ALLOCATE (ftestpt(0:0,size(xgauss, 1))) CALL ftest(xgauss(:, 1), xgauss(:, 2), ftestpt) end if DO igauss = 1, gausssize ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss, 1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j) !CALL coefeqext(xgauss(igauss, :), idt, idw, idg, idp, coefs) DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) newcontrib = 2*pi*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)*wgauss(igauss)*xgauss(igauss, 1)!*wgeom(igauss,0) !$OMP ATOMIC UPDATE volume(mu) = volume(mu) + newcontrib !$OMP END ATOMIC if (walltype .lt. 0) THEN newcontrib = ftestpt(0,igauss)*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)& &*wgeom(0,igauss)*wgauss(igauss)*xgauss(igauss, 1) !$OMP ATOMIC UPDATE fverif(mu) = fverif(mu) + newcontrib !$OMP END ATOMIC end if END DO END DO END DO END DO !$OMP END PARALLEL DO !DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE (f, aux) DEALLOCATE (fun, fun2) call timera(1, "comp_volume") END SUBROUTINE comp_volume !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the gradient of the gtilde function for the web-spline method needed to correctly apply the dirichlet boundary conditions !--------------------------------------------------------------------------- SUBROUTINE comp_gradgtilde USE bsplines USE geometry REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :,:), fun2(:, :,:), gtildeintegr(:, :) Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms Real(kind=db)::newcontrib !call timera(0, "comp_gradgtilde") ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines nterms = 4 Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms)) ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) If (allocated(gtildeintegr)) deallocate (gtildeintegr) ALLOCATE (gtildeintegr(0:2,3*ngauss(1)*ngauss(2))) ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO CALL coefeqext(splrz%sp1%knots(0:1), idt, idw, idg, idp, coefs) !$OMP DO SIMD do j=1,size(gradgtilde) gradgtilde(j) = 0 END DO !$OMP END DO SIMD !$OMP BARRIER ! Assemble gradgtilde matrix !$OMP DO collapse(2) DO j = 1, nz ! Loop on z position DO i = 1, nr ! Loop on r position ! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) iid=i jid=j if (gausssize .gt. 1) then !If (allocated(wgeom)) deallocate (wgeom) !ALLOCATE (wgeom(0:2,gausssize)) CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) Call total_gtilde(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), gtildeintegr(:,1:gausssize),wgeom(:,1:gausssize)) End if DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) newcontrib = 0.0_db DO igauss = 1, gausssize ! Loop on gaussian weights and positions Do iterm = 1, nterms newcontrib = newcontrib + wgeom( idg(iterm),igauss)*gtildeintegr( idp(iterm),igauss)* & & fun(f(jt, 1), idt(iterm),igauss)*fun2(f(jt, 2), idw(iterm),igauss)* & & wgauss(igauss)*xgauss(igauss, 1) End do end do !$OMP ATOMIC UPDATE gradgtilde(mu) = gradgtilde(mu) + newcontrib !$OMP END ATOMIC END DO END DO END DO !$OMP END DO !DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE (f, aux) DEALLOCATE (fun, fun2) !call timera(1, "comp_gradgtilde") END SUBROUTINE comp_gradgtilde !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Imposes the dirichlet boundary conditions on the FEM matrix for the case where we use regular splines ( not web-splines). !--------------------------------------------------------------------------- SUBROUTINE fe_dirichlet REAL(kind=db), ALLOCATABLE :: arr(:) INTEGER :: i ALLOCATE (arr(nrank(1)*nrank(2))) DO i = 1, nrank(1) IF (rgrid(0) .ne. 0.0_db) THEN arr = 0; arr(i) = 1; CALL putrow(femat, i, arr) END IF arr = 0; arr(nrank(1)*nrank(2) + 1 - i) = 1; CALL putrow(femat, nrank(1)*nrank(2) + 1 - i, arr) END DO DEALLOCATE (arr) END SUBROUTINE fe_dirichlet !________________________________________________________________________________ SUBROUTINE coefeq(x, idt, idw, idg, c, kterms) REAL(kind=db), INTENT(in) :: x(:) INTEGER, INTENT(out) :: idt(:, :), idw(:, :), idg(:, :),kterms REAL(kind=db), INTENT(out) :: c(:) kterms=8 c = x(2) idt(1, 1) = 0 idt(1, 2) = 0 idw(1, 1) = 0 idw(1, 2) = 0 idg(1, 1) = 1 idg(1, 2) = 1 idt(2, 1) = 0 idt(2, 2) = 0 idw(2, 1) = 0 idw(2, 2) = 1 idg(2, 1) = 1 idg(2, 2) = 0 idt(3, 1) = 0 idt(3, 2) = 0 idw(3, 1) = 1 idw(3, 2) = 0 idg(3, 1) = 0 idg(3, 2) = 1 idt(4, 1) = 0 idt(4, 2) = 0 idw(4, 1) = 1 idw(4, 2) = 1 idg(4, 1) = 0 idg(4, 2) = 0 idt(5, 1) = 0 idt(5, 2) = 0 idw(5, 1) = 0 idw(5, 2) = 0 idg(5, 1) = 2 idg(5, 2) = 2 idt(6, 1) = 0 idt(6, 2) = 1 idw(6, 1) = 0 idw(6, 2) = 0 idg(6, 1) = 2 idg(6, 2) = 0 idt(7, 1) = 1 idt(7, 2) = 0 idw(7, 1) = 0 idw(7, 2) = 0 idg(7, 1) = 0 idg(7, 2) = 2 idt(8, 1) = 1 idt(8, 2) = 1 idw(8, 1) = 0 idw(8, 2) = 0 idg(8, 1) = 0 idg(8, 2) = 0 END SUBROUTINE coefeq SUBROUTINE coefeqext(x, idt, idw, idg, idp, c) REAL(kind=db), INTENT(in) :: x(:) INTEGER, INTENT(out) :: idp(:), idt(:), idw(:), idg(:) REAL(kind=db), INTENT(out) :: c(:) c(1) = x(2) idp(1) = 1 idg(1) = 1 idt(1) = 0 idw(1) = 0 c(2) = x(2) idp(2) = 1 idg(2) = 0 idt(2) = 0 idw(2) = 1 c(3) = x(2) idp(3) = 2 idg(3) = 2 idt(3) = 0 idw(3) = 0 c(4) = x(2) idp(4) = 2 idg(4) = 0 idt(4) = 1 idw(4) = 0 END SUBROUTINE coefeqext !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Free the memory used by the fields module !--------------------------------------------------------------------------- SUBROUTINE clean_fields Use bsplines USE basic, ONLY: rhs INTEGER:: i do i = 1, nrank(1)*nrank(2) call omp_destroy_lock(mu_lock(i)) end do DEALLOCATE (mu_lock) DEALLOCATE (matcoef) DEALLOCATE (pot) DEALLOCATE (rhs) DEALLOCATE (loc_rhs) DEALLOCATE (loc_moments) DEALLOCATE (phi_spline) DEALLOCATE (Br, Bz) DEALLOCATE (Er, Ez) DEALLOCATE (vec1, vec2) Call DESTROY_SP(splrz) Call DESTROY_SP(splrz_ext) END SUBROUTINE clean_fields SUBROUTINE updt_sploc(arow, j, val) ! ! Update element j of row arow or insert it in an increasing "index" ! USE sparse TYPE(sprow), TARGET :: arow INTEGER, INTENT(in) :: j DOUBLE PRECISION, INTENT(in) :: val ! TYPE(elt), TARGET :: pre_root TYPE(elt), POINTER :: t, p ! if(val.eq.0) return pre_root%next => arow%row0 ! pre_root is linked to the head of the list. t => pre_root DO WHILE (ASSOCIATED(t%next)) p => t%next IF (p%index .EQ. j) THEN p%val = p%val + val RETURN END IF IF (p%index .GT. j) EXIT t => t%next END DO ALLOCATE (p) p = elt(j, val, t%next) t%next => p ! arow%nnz = arow%nnz + 1 arow%row0 => pre_root%next ! In case the head is altered END SUBROUTINE updt_sploc SUBROUTINE updt_ppform2d(sp,c) use bsplines TYPE(spline2d), INTENT(inout) :: sp DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: c !DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:) INTEGER:: m,mm INTEGER :: d1, d2, k1, k2, n1, n2 d1 = sp%sp1%dim d2 = sp%sp2%dim k1 = sp%sp1%order k2 = sp%sp2%order n1 = sp%sp1%nints n2 = sp%sp2%nints !$OMP SINGLE IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs) ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2))) !$OMP END SINGLE !ALLOCATE(work(d2,k1,n1)) !$OMP DO DO m=1,SIZE(c,2) CALL topp0(sp%sp1, c(:,m), ppformwork(m,:,:)) sp%bcoefs(:,m)=c(:,m) END DO !$OMP END DO NOWAIT !$OMP SINGLE IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform) ALLOCATE(sp%ppform(k1,n1,k2,n2)) !$OMP END SINGLE !$OMP DO DO mm=1,SIZE(ppformwork,3) DO m=1,SIZE(ppformwork,2) CALL topp0(sp%sp2, ppformwork(:,m,mm), sp%ppform(m,mm,:,:)) END DO END DO !$OMP END DO !DEALLOCATE(work) end subroutine updt_ppform2d !=========================================================================== SUBROUTINE topp0(sp, c, ppform) ! ! Compute PPFORM of a fuction defined by the spline SP ! and spline coefficients C(1:d) ! use bsplines TYPE(spline1d), INTENT(in) :: sp DOUBLE PRECISION, INTENT(in) :: c(:) DOUBLE PRECISION, INTENT(out) :: ppform(0:,:) INTEGER :: p, nints, i, j, k ! p = sp%order - 1 nints = sp%nints ! ppform = 0.0d0 DO i=1,nints ! on each knot interval DO j=1,p+1 ! all spline in interval i DO k=0,p ! k_th derivatives ppform(k,i) = ppform(k,i) + sp%val0(k,j,i)*c(j+i-1) END DO END DO END DO ! END SUBROUTINE topp0 !+ END MODULE fields diff --git a/src/ion_induced_mod.f90 b/src/ion_induced_mod.f90 index 64a6936..c49920d 100644 --- a/src/ion_induced_mod.f90 +++ b/src/ion_induced_mod.f90 @@ -1,447 +1,534 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: IIEE ! !> @author !> S. Guinchard - EPFL/SPC ! !> Last modif. !> 11/28 2022 ! ! DESCRIPTION: !> Module handling ion induced electron emissions (IIEE) !> following Schou's model (for Kinetic emissions) and Auger neutralisation !> for potential emissions !------------------------------------------------------------------------------ MODULE iiee USE particletypes USE constants USE basic USE materials !#include "mkl_vsl.f90" ! for random # generators using MKL intel library IMPLICIT NONE CONTAINS !--------------------------------------------------------------------------- ! SUBROUTINE ion_induced(pion, losthole, pelec, nblostparts) ! ! ! DESCRIPTION !> function to determine the number of electrons !> to add to a given species as a function of th number !> of lost ions ! !-------------------------------------------------------------------------- SUBROUTINE ion_induced(pion, losthole, pelec, nblostparts) USE geometry TYPE(particles), INTENT(INOUT):: pion, pelec !< ion and electrons parts REAL(KIND = db), DIMENSION(3) :: last_pos !< last position for lost ion (revert push) REAL(KIND = db), DIMENSION(3) :: normal_dir !< normal direction vector (normalised) INTEGER, DIMENSION(pion%Nploc):: losthole !< indices of lost ions INTEGER ::i,j, nblostparts, Nploc, Nploc_old, Nploc_init!< loop indices and #lost particles INTEGER :: parts_size_increase, nbadded INTEGER :: neuttype_id !< neutral gas type_id INTEGER :: material_id !< electrode material type_id INTEGER :: gen_el, kmax !< # of electrons generated, max# possibly gen. elec. REAL(KIND = db) :: lambda !< Poisson param. to gen elec. (yield) REAL(KIND = db) :: kappa, theta, Emax !< gamma distribution parameters REAL(KIND = db) :: Ekin, Eem !< kinetic energy of lost particles (yield param) and of emitted electrons kmax = 12 !> Max num. elec. to be generated (Poisson) kappa = 4.0 !> kappa param. (Gamma) theta = 0.5 !> theta param. (Gamma) Emax = 25 !> Max value for el. (Gamma) neuttype_id = pion%neuttype_id !> temporarily stored in particle type material_id = pion%material_id !> temporarily stored in particle type IF(pelec%Nploc + 2*nblostparts .gt. size(pelec%pos,2)) THEN parts_size_increase=Max(floor(0.1*size(pelec%pos,2)),2*nblostparts) CALL change_parts_allocation(pelec, parts_size_increase) END IF Nploc_init = pelec%Nploc DO i=1,nblostparts Ekin = compute_Ekin( pion%U(:,losthole(i)), pion) IF (pion%neuttype_id .eq. 1) THEN ! Yield for H2 lambda = 2.0*compute_yield(Ekin/2.0, neuttype_id, material_id) ! computed from yield by protons ELSE ! Yield for other neutral lambda = compute_yield(Ekin, neuttype_id, material_id) ! needs to be changed for other gases END IF nbadded = gen_elec(lambda, kmax) Nploc_old = pelec%Nploc pelec%Nploc = pelec%Nploc + nbadded Nploc = pelec%Nploc last_pos = revert_push(pion, losthole(i)) pelec%nbadded = pelec%nbadded+nbadded normal_dir = find_normal(last_pos) DO j=1,nbadded pelec%pos(1,Nploc_old+j) = last_pos(1) pelec%pos(2,Nploc_old+j) = last_pos(2) pelec%pos(3,Nploc_old+j) = last_pos(3) pelec%newindex = pelec%newindex + 1 pelec%partindex(Nploc_old+j) = pelec%newindex IF(pelec%zero_vel .eqv. .false.) THEN Eem = gen_E_gamma(kappa, theta, Emax) !> generate an energy value following gamma distribution pelec%U(1,Nploc_old+j) = compute_Vnorm(Eem, pelec)* normal_dir(1) !> Vr pelec%U(3,Nploc_old+j) = compute_Vnorm(Eem, pelec)* normal_dir(3) !> Vz pelec%U(2,Nploc_old+j) = compute_Vnorm(Eem, pelec)* normal_dir(2) !> Vthet ELSE pelec%U(1,Nploc_old+j) = 0.0 pelec%U(3,Nploc_old+j) = 0.0 pelec%U(2,Nploc_old+j) = 0.0 END IF END DO END DO END SUBROUTINE ion_induced !--------------------------------------------------------------------------- ! FUNCTION compute_Ekin(velocity, p) ! ! ! DESCRIPTION !> Computes the kinetic energy in MeV of a particle given its 3-vel. components !-------------------------------------------------------------------------- FUNCTION compute_Ekin(velocity, p) RESULT(Ekin) TYPE(particles), INTENT(INOUT):: p REAL(KIND = db), DIMENSION(3) :: velocity REAL(KIND = db) :: Ekin Ekin = 5E-7 * p%m * vlight**2 /elchar * (velocity(1)**2 + velocity(2)**2 + velocity(3)**2) END FUNCTION compute_Ekin !--------------------------------------------------------------------------- ! FUNCTION compute_Vnorm(Ekin,p) ! ! ! DESCRIPTION !> Computes the normal velocity of an incident electron emitted !> with energy Ekin !-------------------------------------------------------------------------- FUNCTION compute_Vnorm(Ekin, p) RESULT(Vnorm) REAL(KIND = db) :: Ekin, Vnorm !> Ekin of emitted electron, Normal. corres. veloc. TYPE(particles) :: p !> electrons Vnorm = sqrt(2/p%m * Ekin * elchar) / vlight !> * elchar to get the enery in J and Vnorm in m/s END FUNCTION compute_Vnorm !--------------------------------------------------------------------------- ! FUNCTION fin_normal(last_position) ! ! ! DESCRIPTION !> Computes the normal velocity of an incident electron emitted !> with energy Ekin !-------------------------------------------------------------------------- FUNCTION find_normal(last_position) RESULT(normal_dir) USE geometry REAL(KIND = db), DIMENSION(3) :: last_position !> Last pos. to eval. geom. weight at REAL(KIND = db), DIMENSION(3) :: normal_dir !> Normal direction vector (Result) REAL(KIND = db), DIMENSION(3) :: weight !> Geom. weight at last pos. REAL(KIND = db) :: norm !> To normalise normal vect. call geom_weight(last_position(1), last_position(3), weight) norm = sqrt(weight(2)**2 + weight(3)**2) normal_dir(1) = 1/norm * weight(3) !> Normal along r normal_dir(2) = 0.0 !> Normal along theta normal_dir(3) = 1/norm * weight(2) !> Normal along z END FUNCTION find_normal !--------------------------------------------------------------------------- ! FUNCTION revert_push(pion, partid) ! ! ! DESCRIPTION !> reverts Buneman algorithm over one time step !> to obtain ion position right before recapture !> and hence new electron positions ! !-------------------------------------------------------------------------- FUNCTION revert_push(pion, partid) USE basic, ONLY: dt, tnorm REAL(KIND=db), DIMENSION(3):: revert_push TYPE(particles), INTENT(INOUT):: pion !> species: ions INTEGER :: partid !> id of particle to reverse position revert_push(1) = abs(pion%pos(1,partid) - pion%U(1,partid)*dt) revert_push(2) = pion%pos(2,partid) -1/pion%pos(1,partid)* pion%U(2,partid)*dt revert_push(3) = pion%pos(3,partid) -pion%U(3,partid)*dt END FUNCTION revert_push !--------------------------------------------------------------------------- ! FUNCTION eval_polynomial(coefficients, valeur) ! ! ! DESCRIPTION !> Evaluate a polynomial at a given point !> with its coefficients provided in an array !> s.t lowest order coeff = 1st element ! !-------------------------------------------------------------------------- REAL(KIND = db) FUNCTION eval_polynomial(coefficients, valeur) REAL(KIND = db), DIMENSION(:) :: coefficients !< polynomial (e.g fitted yield) coeffs REAL(KIND = db) :: valeur !< point where to evaluate polyn INTEGER :: ii eval_polynomial = 0 DO ii=1, size(coefficients) eval_polynomial = eval_polynomial+coefficients(ii)*valeur**(ii-1) END DO END FUNCTION eval_polynomial !--------------------------------------------------------------------------- ! FUNCTION compute_yield(energy, neuttype_id, material_id) ! ! ! DESCRIPTION !> Gives the theoretical value for the electron yield !> as a function of the energy of the incident ion and !> the type of neutral gas ! !-------------------------------------------------------------------------- REAL(KIND = db) FUNCTION compute_yield(energy, neuttype_id, material_id) !add material id asap REAL(KIND = db) :: energy INTEGER :: neuttype_id, material_id REAL(KIND = db) :: Lambda_exp Lambda_exp = 1E-3 SELECT CASE(material_id) CASE(1) !304 stainless steel IF(energy.le. 1E-3 ) THEN SELECT CASE(neuttype_id) CASE(1) compute_yield = eval_polynomial(coefficients_1H_SS, energy) CASE(2) compute_yield = eval_polynomial(coefficients_1He_SS, energy) CASE(3) compute_yield = eval_polynomial(coefficients_1Ne_SS, energy) CASE DEFAULT compute_yield = eval_polynomial(coefficients_1H_SS, energy) END SELECT ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_2_SS,energy) + compute_yield = Lambda_exp * eval_polynomial(coefficients_2H_SS,energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_3_SS,energy) + compute_yield = Lambda_exp * eval_polynomial(coefficients_3H_SS,energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_4_SS,energy) + compute_yield = Lambda_exp * eval_polynomial(coefficients_4H_SS,energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_5_SS,energy) + compute_yield = Lambda_exp * eval_polynomial(coefficients_5H_SS,energy) END IF CASE(2) ! Copper IF(energy.le. 1E-3 ) THEN SELECT CASE(neuttype_id) CASE(1) compute_yield = eval_polynomial(coefficients_1H_Cu, energy) CASE(2) compute_yield = eval_polynomial(coefficients_1He_Cu, energy) CASE(3) compute_yield = eval_polynomial(coefficients_1Ne_Cu, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_1Ar_Cu, energy) CASE DEFAULT compute_yield = eval_polynomial(coefficients_1H_Cu, energy) END SELECT + ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_2_Cu,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_2H_Cu,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_2He_Cu, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_2Ar_Cu, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_2H_Cu,energy) + END SELECT + ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_3_Cu,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_3H_Cu,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_3He_Cu, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_3Ar_Cu, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_3H_Cu,energy) + END SELECT + ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_4_Cu,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_4H_Cu,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_4He_Cu, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_4Ar_Cu, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_4H_Cu,energy) + + END SELECT ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_5_Cu,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_5H_Cu,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_5He_Cu, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_5Ar_Cu, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_5H_Cu,energy) + + END SELECT END IF + CASE(3) ! Alumium IF(energy.le. 1E-3 ) THEN SELECT CASE(neuttype_id) CASE(1) compute_yield = eval_polynomial(coefficients_1H_Al, energy) CASE(2) compute_yield = eval_polynomial(coefficients_1He_Al, energy) CASE(3) compute_yield = eval_polynomial(coefficients_1Ne_Al, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_1Ar_Al, energy) CASE DEFAULT compute_yield = eval_polynomial(coefficients_1H_Al, energy) END SELECT + ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_2_Al,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_2H_Al,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_2He_Al, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_2Ar_Al, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_2H_Al,energy) + + END SELECT ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_3_Al,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_3H_Al,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_3He_Al, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_3Ar_Al, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_3H_Al,energy) + + END SELECT ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_4_Al,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_4H_Al,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_4He_Al, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_4Ar_Al, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_4H_Al,energy) + + END SELECT ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN - compute_yield = Lambda_exp * eval_polynomial(coefficients_5_Al,energy) + SELECT CASE(neuttype_id) + CASE(1) + compute_yield = Lambda_exp * eval_polynomial(coefficients_5H_Al,energy) + CASE(2) + compute_yield = eval_polynomial(coefficients_5He_Al, energy) + CASE(4) + compute_yield = eval_polynomial(coefficients_5Ar_Al, energy) + CASE DEFAULT + compute_yield = Lambda_exp * eval_polynomial(coefficients_5H_Al,energy) + + END SELECT END IF END SELECT END FUNCTION compute_yield !--------------------------------------------------------------------------- ! FUNCTION gen_elec(lambda, kmax) ! ! ! DESCRIPTION !> Gives random values distributed !> following a Poisson distrib. of parameter !> lambda = yield(E) for incomin ion energy E !-------------------------------------------------------------------------- INTEGER FUNCTION gen_elec(lambda, kmax) USE random_mod REAL(KIND = db) :: lambda !< Lambda parameter for Poisson distribution REAL(KIND = db) :: nb_alea(1:1) !< random number unif. generated in [0,1] INTEGER :: kmax !< max number possible from Poisson REAL(KIND = db) :: CumulPoisson !< Flag to ensure CDF ~ 1 INTEGER :: i, ii !< loop indices REAL(KIND = db), DIMENSION(kmax) :: vect, SumPart !< terms, partial sums for CDF !> Compute probabilities for each int. value and CDF values DO i = 1,kmax vect(i) = exp(-lambda)*lambda**(i-1)/factorial_fun(i-1); SumPart(i) = sum(vect(1:i)); END DO !> CDF expected to sum to 1 CumulPoisson = sum(vect) !> Generate poisson distrib. int. (see Matlab. code for convg.) call random_array(nb_alea,1,ran_index(1),ran_array(:,1)) DO ii = 1,size(SumPart)-1 IF (nb_alea(1) .lt. SumPart(1)) THEN gen_elec = 0 ELSE IF ((SumPart(ii).le.nb_alea(1)) .and. (nb_alea(1) .lt. SumPart(ii+1))) THEN gen_elec = ii END IF END DO ! Below: see Intel oneAPI Math Kernel Library - Fortran ! to optimise running speed when generating random numbers ! TO DO : change if enough time !----------------------------------------------------------- ! ! INTEGER, INTENT(IN) :: method ! TYPE (VSL_STREAM_STATE), INTENT(IN) :: stream ! INTEGER, INTENT(IN) :: n = 1 ! INTEGER, INTENT(IN) :: r ! status = virngpoisson( method, stream, n, r, lambda ) ! gen_elec = r ! !------------------------------------------------------------ END FUNCTION gen_elec !--------------------------------------------------------------------------- ! FUNCTION gen_E_gamma(kappa, theta, Emax) ! ! ! DESCRIPTION !> Gives random values distributed !> following a Gamma distrib. of parameters (kappa, theta) !> in [0, Emax] eV and peaked at E=2eV !-------------------------------------------------------------------------- FUNCTION gen_E_gamma(kappa, theta, Emax) RESULT(E_el) USE random_mod USE incomplete_gamma REAL(KIND = db) :: E_el, Emin, Emax REAL(KIND = db) :: kappa, theta !> parameters to shape Gamma_distr REAL(KIND = db) :: nb_alea(1:1) REAL(KIND = db), DIMENSION(1000) :: Einc, CDF, Diff INTEGER :: ii !> loop index INTEGER :: ifault INTEGER :: posid Emin = 0.00 Emax = 25.0 kappa = 4.0 theta = 0.5 DO ii= 1,size(Einc) Einc(ii) = ( Emin + (ii-1)*(Emax-Emin)/(size(Einc)-1) ) CDF(ii) = gamain(Einc(ii)/theta, kappa, ifault) END DO !> Generate Gamma distrib. int. (see Matlab. code for convg.) call random_array(nb_alea,1,ran_index(1),ran_array(:,1)) Diff = abs(nb_alea(1) - CDF) posid = minloc(Diff, dim = 1) E_el = Einc(posid) END FUNCTION gen_E_gamma !--------------------------------------------------------------------------- ! FUNCTION factorial_fun(n) ! ! ! DESCRIPTION !> Gives the factorial of an integer !-------------------------------------------------------------------------- INTEGER FUNCTION factorial_fun(n) INTEGER :: ii INTEGER :: n factorial_fun = 1 IF (n .lt. 0) THEN RETURN ELSE IF (n .eq. 0) THEN factorial_fun = 1 ELSE DO ii=1,n factorial_fun = factorial_fun*ii END DO END IF END FUNCTION factorial_fun END MODULE iiee diff --git a/src/materials_mod.f90 b/src/materials_mod.f90 index 094938b..29dd05b 100644 --- a/src/materials_mod.f90 +++ b/src/materials_mod.f90 @@ -1,69 +1,136 @@ !------------------------------------------------------------------------------ -! EPFL/Swiss Plasma Center +! EPFL/Swiss Plasma Cente !------------------------------------------------------------------------------ ! ! MODULE: materials ! !> @author / date !> S. Guinchard - EPFL/SPC / 11-25-22 ! ! DESCRIPTION: !> Module storing coefficients corresponding to polynomial !> fits of energy loss function dE/dx in several materials !> N.B. These coefficients fit dE/dx in MeV/cm for protons !> only and should be adapted for other ion species !------------------------------------------------------------------------------ MODULE materials USE constants IMPLICIT NONE ! ----------------- 304 STAINLESS STEEL ---------------- !> All the coefficients below were obtained by piecewise !> fit of dE/dx curve for 304 stainless steel REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1H_SS = (/0.0, 2.9778E02 /) REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1He_SS = (/0.1273010048, 1.70478995200000E02 /) REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ne_SS = (/0.0518524160, 2.45927583999999E02 /) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2_SS = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2H_SS = & (/1.834520315818557E02,1.320304216355084E05,-8.583700602370013E06,3.526140145560557E08/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3_SS = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3H_SS = & (/2.471679999999999E02,9.695466666666670E04,-2.475200000000003E06, 2.325333333333340E07/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4_SS = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4H_SS = & (/2.533904454349683E03,-1.766016382825937E05,8.202640024592019E06, -1.125320217235288E08/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_5_SS = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_5H_SS = & (/8.786057142856745E02,2.856595238095524e+04,-1.834285714286391E05, 3.333333333338620E05/) + ! --------------------- COPPER ------------------------- - !> All the coefficients below were obtained by piecewise - !> fit of dE/dx curve for Copper - REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1H_Cu = (/0.0, 2.812300E02 /) - REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1He_Cu = (/0.1273010048, 1.539289952000001E02 /) - REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ne_Cu = (/0.0518524160, 2.293775840000000E02 /) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2_Cu = & + !> All the coefficients for H were obtained by piecewise + !> fit of dE/dx curve for Copper + + !> 1 Hydrogen + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1H_Cu = & + (/0.0, 2.812300E02 /) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2H_Cu = & (/1.703970762187365E02, 1.273920899063999E05, -8.680492697427949E06, 3.535110244304671E08/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3_Cu = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3H_Cu = & (/2.618095714285664E02, 8.520814285714373E04, -2.094971428571471E06, 2.064000000000054E07/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4_Cu = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4H_Cu = & (/1.030650000000107E03, -6.231904761918026E03, 1.468571428571969E06, -2.506666666667393E07/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_5_Cu = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_5H_Cu = & (/7.117599999999986E02, 3.874166666666669E04, -5.099999999999992E05, 2.733333333333320E06/) - ! --------------------- ALUMINUM ---------------------- + !> All the coefficients below were obtained by piecewise + !> fit of yield curve for Copper + + !> 2 Helium (https://doi.org/10.1103/PhysRevB.19.121) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1He_Cu = & + (/0.3165714286, 0.0/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2He_Cu = & + (/1.83738612473721E-01, 7.95824573697734E01, -3.9809623919645E03, 2.18406914272366E05/) + REAL(KIND = db), DIMENSION(3), PARAMETER :: coefficients_3He_Cu = & + (/2.000000000000011E-02, 9.799999999999999E01, -2.0E03/) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_4He_Cu = & + (/6.6E-01, 2.6E01/) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_5He_Cu = & + (/6.6E-01, 2.6E01/) + + !> 3 Neon + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ne_Cu = & + (/0.0518524160, 2.293775840000000E02 /) + + !> 4 Argon (https://doi.org/10.1016/0039-6028(79)90341-8) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ar_Cu = & + (/0.1108571429, 0.0/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2Ar_Cu = & + (/1.22823377416242E-02, 6.22563244459554E01, 7.0358339950776E03, -4.26694208031656E05/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3Ar_Cu = & + (/-2.5992776033708E-02, 1.07801648650678E02, -1.41855662651229E03, 4.91523169780497E03/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4Ar_Cu = & + (/-3.87791890832554E-01, 1.48442167310384E02, -3.1298539515514E03, 3.39747092112027E04/) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_5Ar_Cu = & + (/5.78008566722373E-01, 5.29272726756834E01/) + + + + ! --------------------- ALUMINUM ---------------------- + !> All the coefficients for H were obtained by piecewise !> fit of dE/dx curve for Aluminum - REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1H_Al = (/0.0, 2.100900000000000E02 /) - REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1He_Al = (/0.1273010048, 82.788995200000016 /) - REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ne_Al = (/0.0518524160, 1.582375839999999E02 /) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2_Al = & - (/1.224483183711612E02, 9.731277035468460E04, -6.878546755544925E06, 2.795627748449915E08 /) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3_Al = & + + !> 1 Hydrogen + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1H_Al = & + (/0.0, 2.100900000000000E02 /) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2H_Al = & + (/1.224483183711612E02, 9.731277035468460E04, -6.878546755544925E06, 2.795627748449915E08 /) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3H_Al = & (/1.998525714285755E02, 6.229714285714197E04, -1.531771428571365E06, 1.567999999999851E07/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4_Al = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4H_Al = & (/8.135440000000424E02, -1.333600000000502E04, 1.553600000000198E06, -2.624000000000260E07/) - REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_5_Al = & + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_5H_Al = & (/4.052485714285687E02, 3.974642857142875E04, -6.631428571428614E05, 3.800000000000037E06 /) + !> All the coefficients below were obtained by piecewise + !> fit of yield curve for Aluminum + + !> 2 Helium (https://doi.org/10.1103/PhysRevB.19.121) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1He_Al = & + (/0.1757264957, 54.2735043/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2He_Al = & + (/1.76659521215925E-01, 5.34536358214894E01, -6.19216675102752E02, -2.93593275345582E04/) + REAL(KIND = db), DIMENSION(3), PARAMETER :: coefficients_3He_Al = & + (/2.550000000000002E-01, 4.14999999999998E01, -4.999999999999995E02/) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_4He_Al = & + (/4.55E-01, 2.15E01/) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_5He_Al = & + (/4.7E-01, 2.1E01/) + + !> 3 Neon + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ne_Al = & + (/0.0518524160, 1.582375839999999E02 /) + + !> 4 Argon (https://doi.org/10.1016/0039-6028(79)90341-8) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_1Ar_Al = & + (/0.0526495726, 3.98687652E01/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_2Ar_Al = & + (/9.54897034753516E-02, -1.76674600320705E01, 1.54604387031526E04, -7.64344295477741E05/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_3Ar_Al = & + (/-3.26838161662861E-01, 1.46366287920992E02, -5.18224285283153E03, 8.11411852216318E04/) + REAL(KIND = db), DIMENSION(4), PARAMETER :: coefficients_4Ar_Al = & + (/5.51191348971975E-02, 7.50010618137324E01, -1.14344429560049E03, 9.72498529527904E03/) + REAL(KIND = db), DIMENSION(2), PARAMETER :: coefficients_5Ar_Al = & + (/5.25429038367765E-01, 3.37724550898206E01/) END MODULE materials diff --git a/src/restart.f90 b/src/restart.f90 index 6bcc68f..d6680ce 100644 --- a/src/restart.f90 +++ b/src/restart.f90 @@ -1,35 +1,35 @@ SUBROUTINE restart ! ! Additional data specific for a restart run ! use BASIC use mpi use constants IMPLICIT NONE Real(kind=db):: file_version=1.0 ! ! Local vars and arrays !________________________________________________________________________________ ! WRITE(*,'(a/)') '=== Define additional data for a restart run ===' !________________________________________________________________________________ - IF (mpirank .eq. 0) then + IF (mpirank .eq. 0 .and. newres.eq..false.) then CALL openf(resfile, fidres, real_prec='d') CALL getatt(fidres,"/data","file_version",file_version) if (file_version .lt. 2.0) then WRITE(*,*) " " WRITE(*,*) " " WRITE(*,*) " " WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" WRITE(*,*) "ERROR: This simulation was run with FENNECS version 1.0!" WRITE(*,*) "You need to use newres=1 to continue the simulation" WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" WRITE(*,*) " " WRITE(*,*) " " WRITE(*,*) " " call closef(fidres) CALL MPI_abort(MPI_COMM_WORLD, -1, ierr) end if end if ! END SUBROUTINE restart diff --git a/src/stepon.f90 b/src/stepon.f90 index 39d72d6..40e7ab4 100644 --- a/src/stepon.f90 +++ b/src/stepon.f90 @@ -1,140 +1,141 @@ SUBROUTINE stepon ! ! Advance one time step ! USE basic USE constants USE fields USE maxwsrce USE celldiag USE neutcol USE sort Use psupply use omp_lib USE beam IMPLICIT NONE INTEGER:: i !$OMP PARALLEL DEFAULT(shared), private(i) DO i=1,nbspecies ! Boundary conditions for plasma particles outside the plasma region CALL bound(partslist(i)) END DO !$OMP BARRIER DO i=1,nbspecies ! Localisation of particles in cells (calculation of the r and z indices) call boundary_loss(partslist(i)) END DO !$OMP BARRIER + ! 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 !$OMP BARRIER !$OMP SINGLE ! The particles are injected by the source CALL maxwsrce_inject(time) !$OMP END SINGLE !$OMP BARRIER !$OMP END PARALLEL ! We sort the particles according to their linear index !IF(modulo(step,100) .eq. 0) THEN ! DO i=1,nbspecies ! ! Boundary conditions for plasma particles outside the plasma region ! call gridsort(partslist(i), 1, partslist(i)%Nploc) ! END DO !end if ! update the power supply voltage if necessary call psupply_step(the_ps,partslist,cstep) !$OMP PARALLEL ! Assemble right hand side of Poisson equation CALL rhscon(partslist) !$OMP END PARALLEL if (.not. nlfreezephi) THEN ! Solve Poisson equation !!$OMP MASTER CALL poisson(splrz,reducedsol) !!$OMP END MASTER end if !$OMP PARALLEL DEFAULT(shared), private(i) DO i=1,nbspecies ! Compute the magnetic field at the particle position call comp_mag_p(partslist(i)) END DO if (.not. nlfreezephi) THEN call poisson_com(splrz,reducedsol) CALL Update_phi(splrz) end if DO i=1,nbspecies ! Compute the electric field at the particle position CALL EFieldscompatparts(partslist(i)) END DO !$OMP BARRIER DO i=1,nbspecies ! 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 !$OMP BARRIER - ! Cell diag quantities - IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN - CALL celldiag_save(time, fidres) - END IF - ! 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 !$OMP BARRIER !$OMP MASTER ! Save variables to file CALL diagnose(step) !$OMP END MASTER !IF (modulo(step,itparts).eq. 0 .or. modulo(step,ittracer).eq. 0 .or. modulo(step,itrestart).eq. 0 .or. nlend) THEN !$OMP BARRIER !END IF Do i=1,nbspecies ! Calculate new positions of particles at time t+delta t CALL push(partslist(i)) END DO !$OMP END PARALLEL ! We recalculate the mpi axial boundaries and we adapt them if necessary IF(modulo(step,100) .eq. 0) THEN CALL calc_Zbounds(partslist(1),Zbounds, femorder) CALL fields_comm_init(Zbounds) CALL maxwsrce_calcfreq(Zbounds) Do i=1,nbspecies if( partslist(i)%Nploc*2 .lt. size(partslist(i)%pot,1)) then call change_parts_allocation(partslist(i),-(size(partslist(i)%pot,1)-partslist(i)%Nploc*2)) end if end do END IF END SUBROUTINE stepon