diff --git a/src/auxval.f90 b/src/auxval.f90 index dde46cb..50c2e28 100644 --- a/src/auxval.f90 +++ b/src/auxval.f90 @@ -1,114 +1,114 @@ SUBROUTINE auxval ! USE constants USE basic USE fields USE beam ! ! Set auxiliary values ! IMPLICIT NONE ! INTEGER:: i !________________________________________________________________________________ IF(mpirank .eq. 0) WRITE(*,'(a/)') '=== Set auxiliary values ===' !________________________________________________________________________________ ! ! Total number of intervals nr=sum(nnr) ! Normalisation constants qsim=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)*n0*elchar/nplasma msim=abs(qsim)/elchar*me vnorm=vlight tnorm=1/sqrt(abs(n0)*elchar*elchar/me/eps_0) rnorm=vnorm*tnorm bnorm=B0 enorm=vlight*bnorm phinorm=enorm*rnorm ! Normalised boundary conditions potinn=potinn/phinorm potout=potout/phinorm ! Characteristic frequencies and normalised volume omegac=sign(elchar,qsim)/me*B0 omegap=sqrt(elchar**2*abs(n0)/me/eps_0) V=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)/rnorm/rnorm/rnorm IF(mpirank .eq. 0) THEN IF(omegap.GT.omegac) THEN WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegap/omegac', omegap, omegac, omegap/omegac ELSE WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegac/omegap', omegap, omegac, omegac/omegap END IF END IF ! Construction of the mesh rgrid in r and zgrid in z and its normalisation CALL mesh rgrid=rgrid/rnorm zgrid=zgrid/rnorm dz=dz/rnorm dr=dr/rnorm ! Define Zindex bounds for the MPI process in the case of sequential run ALLOCATE( Zbounds(0:mpisize), Nplocs_all(0:mpisize-1)) - Zbounds(0) = 0 + Zbounds=0 Zbounds(mpisize) = nz IF(mpisize.gt.1) THEN leftproc=mpirank-1 IF(mpirank .eq. 0 .AND. partperiodic) THEN leftproc = mpisize-1 ELSE IF (mpirank .eq. 0) THEN leftproc=-1 END IF rightproc = mpirank+1 IF(mpirank .eq. mpisize-1 .AND. partperiodic) THEN rightproc = 0 ELSE IF (mpirank .eq. mpisize-1) THEN rightproc=-1 END IF ELSE leftproc=-1 rightproc=-1 END IF WRITE(*,*) "mpirank, left, right", mpirank, leftproc, rightproc ! Initialize the neighbors communicators !CALL init_mpicomms(mpirank, mpisize, leftproc, rightproc) ! Auxiliary vectors ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1))) DO i=0,nr vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz) vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i) END DO !________________________________________________________________________________ CONTAINS !________________________________________________________________________________ SUBROUTINE mesh USE basic, ONLY: zgrid, rgrid, nr, nnr, nz, dr, dz, lz, radii INTEGER :: j ALLOCATE(zgrid(0:nz),rgrid(0:nr)) dz=(lz(2)-lz(1))/nz DO j=0,nz zgrid(j)=j*dz+lz(1) END DO dr(1)=(radii(2)-radii(1))/nnr(1) DO j=0,nnr(1) rgrid(j)=radii(1)+j*dr(1) END DO dr(2)=(radii(3)-radii(2))/nnr(2) DO j=1,nnr(2) rgrid(nnr(1)+j)=radii(2)+j*dr(2) END DO dr(3)=(radii(4)-radii(3))/nnr(3) DO j=1,nnr(3) rgrid(nnr(1)+nnr(2)+j)=radii(3)+j*dr(3) END DO END SUBROUTINE mesh !________________________________________________________________________________ END SUBROUTINE auxval diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index 33003f1..40a3e68 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,323 +1,323 @@ MODULE basic ! USE hashtable USE constants USE bsplines USE mumps_bsplines USE futils USE mpihelper IMPLICIT NONE ! ! Basic module for time dependent problems ! CHARACTER(len=80) :: label1, label2, label3, label4 ! ! BASIC Namelist ! LOGICAL :: nlres = .FALSE. !< Restart flag LOGICAL :: nlsave = .TRUE. !< Checkpoint (save) flag LOGICAL :: newres=.FALSE. !< New result HDF5 file LOGICAL :: nlxg=.FALSE. !< Show graphical interface Xgrafix LOGICAL :: nlmaxwellsource = .FALSE. !< Activate the maxwell source INTEGER :: nrun=1 !< Number of time steps to run REAL(kind=db) :: job_time=3600.0 !< Time allocated to this job in seconds REAL(kind=db) :: tmax=100000.0 !< Maximum simulation time REAL(kind=db) :: extra_time=60.0 !< Extra time allocated REAL(kind=db) :: dt=1 !< Time step REAL(kind=db) :: time=0 !< Current simulation time (Init from restart file) ! ! Other basic global vars and arrays ! INTEGER :: jobnum !< Job number INTEGER :: step !< Calculation step of this run INTEGER :: cstep=0 !< Current step number (Init from restart file) LOGICAL :: nlend !< Signal end of run INTEGER :: ierr !< Integer used for MPI INTEGER :: it0d=1 !< Number of iterations between 0d values writes to hdf5 INTEGER :: it2d=1 !< Number of iterations between 2d values writes to hdf5 INTEGER :: itparts=1 !< Number of iterations between particles values writes to hdf5 INTEGER :: ittext=1 !< Number of iterations between text outputs in the console INTEGER :: itrestart=10000 !< Number of iterations between save of restart.h5 file INTEGER :: itgraph !< Number of iterations between graphical interface updates INTEGER :: mpirank !< MPIrank of the current processus INTEGER :: mpisize !< Size of the MPI_COMM_WORLD communicator INTEGER :: rightproc !< Rank of next processor in the z decomposition INTEGER :: leftproc !< Rank of previous processor in the z decomposition ! ! List of logical file units INTEGER :: lu_in = 90 !< File duplicated from STDIN INTEGER :: lu_stop = 91 !< stop file, see subroutine TESEND ! ! HDF5 file CHARACTER(len=64) :: resfile = "results.h5" !< Main result file CHARACTER(len=64) :: rstfile = "restart.h5" !< Restart file INTEGER :: fidres !< File ID for resfile INTEGER :: fidrst !< File ID for restart file TYPE(BUFFER_TYPE) :: hbuf0 !< Hashtable for 0d var ! ! Plasma parameters LOGICAL :: nlPhis= .TRUE. !< Calculate self consistent electric field flag LOGICAL :: nlclassical= .FALSE. !< If true, solves the equation of motion according to classical !! dynamics LOGICAL :: nlperiod(2)=(/.false.,.false./)!< Set periodic splines on or off LOGICAL :: partperiodic= .TRUE. !< Sets if the particles boundary conditions are periodic or open INTEGER :: nplasma !< Number of macro-particles on initialisation INTEGER :: npartsalloc = 0 !< Size of particle memory allocated at the begining of the simulation INTEGER :: nblock !< Number of slices in Z for stable distribution initialisation REAL(kind=db) :: potinn !< Electric potential at the inner metallic wall REAL(kind=db) :: potout !< Electric potential at the outer metallic wall REAL(kind=db) :: B0 !< Max magnitude of magnetic field REAL(kind=db), allocatable :: Bz(:), Br(:) !< Magnetic field components REAL(kind=db), allocatable :: Athet(:) !< Theta component of the magnetic vector potential TYPE(spline2d) :: splrz !< Spline at r and z REAL(kind=db), allocatable :: Ez(:), Er(:) !< Electric field components REAL(kind=db), allocatable :: pot(:) !< Electro static potential REAL(kind=db) :: radii(4) !< Inner and outer radius of cylinder and radii of fine mesh region REAL(kind=db) :: plasmadim(4) !< Zmin Zmax Rmin Rmax values for plasma particle loading INTEGER :: distribtype=1 !< Type of distribution function used to load the particles !!1: gaussian, 2: Stable as defined in 4.95 of Davidson REAL(kind=db) :: H0=0 !< Initial value of Hamiltonian for distribution 2 REAL(kind=db) :: P0=0 !< Initial canonical angular momentum for distribution 2 REAL(kind=db) :: temprescale = -1.0 !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling) INTEGER :: samplefactor =-1 !< Factor used for the up-sampling of the particles number REAL(kind=db) :: lz(2) !< Lower and upper cylinder limits in z direction REAL(kind=db) :: n0 !< Physical plasma density parameter - REAL(kind=db), ASYNCHRONOUS, DIMENSION(:,:), ALLOCATABLE, SAVE:: moments !< Moments of the distribution function evaluated every it2d + REAL(kind=db), DIMENSION(:,:), ALLOCATABLE, SAVE:: moments !< Moments of the distribution function evaluated every it2d REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: rhs !< right hand side of the poisson equation solver REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: volume !< Volume covered by each spline for density calculation INTEGER :: nz !< Number of grid intervals in z INTEGER :: nr !< Total number of grid intervals in r INTEGER :: nnr(3) !< Number of grid intervals in r in each subdomain REAL(kind=db) :: dz !< Cell size in z REAL(kind=db) :: dr(3) !< Cell size in r for each region REAL(kind=db), ALLOCATABLE :: zgrid(:) !< Nodes positions in longitudinal direction REAL(kind=db), ALLOCATABLE :: rgrid(:) !< Nodes positions in radial direction REAL(kind=db) :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm !< Normalization constants REAL(kind=db) :: qsim !< Charge of superparticles REAL(kind=db) :: msim !< Mass of superparticles INTEGER :: femorder(2) !< FEM order INTEGER :: ngauss(2) !< Number of gauss points LOGICAL :: nlppform =.TRUE. !< Argument of set_spline INTEGER, SAVE :: nrank(2) !< Number of splines in both directions REAL(kind=db) :: omegac !< Cyclotronic frequency REAL(kind=db) :: omegap !< Plasma frequency REAL(kind=db) :: V !< Normalised volume REAL(kind=db) :: temp !< Initial temperature of plasma REAL(kind=db) :: Rcurv !< Magnetic field curvature coefficient REAL(kind=db) :: Width !< Distance between two magnetic mirrors INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction TYPE(BASICDATA) :: bdata !< Structure used for the mpi communication of the run parameters CONTAINS ! !================================================================================ SUBROUTINE basic_data ! ! Define basic data ! use mpihelper IMPLICIT NONE ! ! Local vars and arrays CHARACTER(len=256) :: line, inputfilename ! NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg, & & nplasma, potinn, potout, B0, lz, n0, nz, nnr, femorder, ngauss, & & nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, & & resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0, partperiodic, & & temprescale, samplefactor, nlmaxwellsource, npartsalloc !________________________________________________________________________________ ! 1. Process Standard Input File ! IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING' STOP ENDIF CALL GET_COMMAND_ARGUMENT(1,inputfilename) OPEN(UNIT=lu_in,FILE=trim(inputfilename),ACTION='READ') IF(mpirank .eq. 0) THEN !________________________________________________________________________________ ! 1. Label the run ! READ(lu_in,'(a)') label1 READ(lu_in,'(a)') label2 READ(lu_in,'(a)') label3 READ(lu_in,'(a)') label4 ! WRITE(*,'(12x,a/)') label1(1:len_trim(label1)) WRITE(*,'(12x,a/)') label2(1:len_trim(label2)) WRITE(*,'(12x,a/)') label3(1:len_trim(label3)) WRITE(*,'(12x,a/)') label4(1:len_trim(label4)) !________________________________________________________________________________ ! 2. Read in basic data specific to run ! READ(lu_in,basic) WRITE(*,basic) ELSE READ(lu_in,basic) END IF ! Distribute run parameters to all MPI workers CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation IF(samplefactor .gt. 1 .and. .not. newres) THEN IF(mpirank.eq.0) WRITE(*,*)"To increase the number of particles, you need to create a new result file (set newres to 1)" CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) END IF IF (npartsalloc .lt. nplasma) THEN npartsalloc=nplasma END IF ! END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! ! Print date and time ! IMPLICIT NONE ! CHARACTER(len=*), INTENT(in) :: str ! ! Local vars and arrays CHARACTER(len=16) :: d, t, dat, functime !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12) ! END SUBROUTINE daytim !================================================================================ SUBROUTINE timera(cntrl, str, eltime) ! ! Timers (cntrl=0/1 to Init/Update) ! IMPLICIT NONE INTEGER, INTENT(in) :: cntrl CHARACTER(len=*), INTENT(in) :: str REAL(kind=db), OPTIONAL, INTENT(out) :: eltime ! INTEGER, PARAMETER :: ncmax=128 INTEGER, SAVE :: icall=0, nc=0 REAL(kind=db), SAVE :: startt0=0.0 CHARACTER(len=16), SAVE :: which(ncmax) INTEGER :: lstr, found, i REAL(kind=db) :: seconds REAL(kind=db), DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0 !________________________________________________________________________________ IF( icall .EQ. 0 ) THEN icall = icall+1 startt0 = seconds() END IF lstr = LEN_TRIM(str) IF( lstr .GT. 0 ) found = loc(str) !________________________________________________________________________________ ! SELECT CASE (cntrl) ! CASE(-1) ! Current wall time IF( PRESENT(eltime) ) THEN eltime = seconds() - startt0 ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0 END IF ! CASE(0) ! Init Timer IF( found .EQ. 0 ) THEN ! Called for the 1st time for 'str' nc = nc+1 which(nc) = str(1:lstr) found = nc END IF startt(found) = seconds() ! CASE(1) ! Update timer endt(found) = seconds() - startt(found) IF( PRESENT(eltime) ) THEN eltime = endt(found) ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found) END IF ! CASE(2) ! Update and reset timer endt(found) = endt(found) + seconds() - startt(found) startt(found) = seconds() IF( PRESENT(eltime) ) THEN eltime = endt(found) END IF ! CASE(9) ! Display all timers IF( nc .GT. 0 ) THEN WRITE(*,'(a)') "Timer Summary" WRITE(*,'(a)') "=============" DO i=1,nc WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i) END DO END IF ! END SELECT ! CONTAINS INTEGER FUNCTION loc(str) CHARACTER(len=*), INTENT(in) :: str INTEGER :: i, ind loc = 0 DO i=1,nc ind = INDEX(which(i), str(1:lstr)) IF( ind .GT. 0 .AND. LEN_TRIM(which(i)) .EQ. lstr ) THEN loc = i EXIT END IF END DO END FUNCTION loc END SUBROUTINE timera !================================================================================ SUBROUTINE readbdata nlres = bdata%nlres nlsave = bdata%nlsave newres = bdata%newres nlxg = bdata%nlxg nlppform = bdata%nlppform nlPhis = bdata%nlPhis nlclassical = bdata%nlclassical nplasma = bdata%nplasma nz = bdata%nz it0d = bdata%it0d it2d = bdata%it2d itparts = bdata%itparts itgraph = bdata%itgraph distribtype = bdata%distribtype nblock = bdata%nblock nrun = bdata%nrun job_time = bdata%job_time extra_time = bdata%extra_time tmax = bdata%tmax dt = bdata%dt potinn = bdata%potinn potout = bdata%potout B0 = bdata%B0 n0 = bdata%n0 temp = bdata%temp Rcurv = bdata%Rcurv width = bdata%width H0 = bdata%H0 P0 = bdata%P0 femorder = bdata%femorder ngauss = bdata%ngauss nnr = bdata%nnr lz = bdata%lz radii = bdata%radii plasmadim = bdata%plasmadim resfile = bdata%resfile partperiodic = bdata%partperiodic samplefactor = bdata%samplefactor END SUBROUTINE readbdata !================================================================================ END MODULE basic diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index b5749eb..3fdf012 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,1619 +1,1620 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Guillaume Le Bars EPFL/SPC !> Patryk Kaminski EPFL/SPC !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> Module responsible for loading, advancing and computing the necessary diagnostics for the simulated particles. !------------------------------------------------------------------------------ MODULE beam ! USE constants USE mpi USE mpihelper USE basic, ONLY: mpirank, mpisize USE distrib IMPLICIT NONE !> Stores the particles properties for the run. TYPE particles INTEGER :: Nptot !< Local number of simulated particles REAL(kind=db) :: m !< Macroparticle mass REAL(kind=db) :: q !< Macroparticle charge REAL(kind=db) :: qmRatio !< Charge over mass ratio INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex !< Index in the electric potential grid for the R direction INTEGER, DIMENSION(:), ALLOCATABLE :: Zindex !< Index in the electric potential grid for the Z direction INTEGER, DIMENSION(:), ALLOCATABLE :: partindex !< Index of the particle to be able to follow it when it goes from one MPI host to the other REAL(kind=db), DIMENSION(:), ALLOCATABLE :: R !< radial coordinates of the particles REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Z !< longitudinal coordinates of the particles REAL(kind=db), DIMENSION(:), ALLOCATABLE :: THET !< azimuthal coordinates of the particles REAL(kind=db), DIMENSION(:), ALLOCATABLE :: WZ !< axial radial relative distances to the left grid line REAL(kind=db), DIMENSION(:), ALLOCATABLE :: WR !< radial relative distances to the bottom grid line REAL(kind=db), DIMENSION(:), ALLOCATABLE :: pot !< Electric potential REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Er !< Radial Electric field REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Ez !< Axial electric field REAL(kind=db), DIMENSION(:),POINTER:: UR !< normalized radial velocity at the current time step REAL(kind=db), DIMENSION(:),POINTER:: URold !< normalized radial velocity at the previous time step REAL(kind=db), DIMENSION(:),POINTER:: UTHET !< normalized azimuthal velocity at the current time step REAL(kind=db), DIMENSION(:),POINTER:: UTHETold !< normalized azimuthal velocity at the previous time step REAL(kind=db), DIMENSION(:),POINTER:: UZ !< normalized axial velocity at the current time step REAL(kind=db), DIMENSION(:),POINTER:: UZold !< normalized axial velocity at the previous time step REAL(kind=db), DIMENSION(:),POINTER:: Gamma !< Lorentz factor at the current time step REAL(kind=db), DIMENSION(:),POINTER:: Gammaold !< Lorentz factor at the previous time step END TYPE particles ! TYPE(particles) :: parts !< Storage for all the particles SAVE :: parts TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffers to send and receive particle from left and right processes ! Diagnostics (scalars) REAL(kind=db) :: ekin=0 !< Total kinetic energz (J) REAL(kind=db) :: epot=0 !< Total potential energy (J) REAL(kind=db) :: etot=0 !< Current total energy (J) REAL(kind=db) :: etot0=0 !< Initial total energy (J) REAL(kind=db) :: loc_etot0=0 !< theoretical local total energy (J) INTEGER, DIMENSION(3), SAVE :: ireducerequest=MPI_REQUEST_NULL !< MPI requests used for the IREDUCE in the diagnostic routine ! INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process LOGICAL:: collected=.false. !< Stores if the particles data have been collected to MPI root process during this timestep ! CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Allocate the memory for the particles variable storing the particles quantities. ! !> @param[inout] p the particles variable needing to be allocated. !> @param[in] nparts the maximum number of particles that will be stored in this variable !--------------------------------------------------------------------------- SUBROUTINE creat_parts(p, nparts) TYPE(particles) :: p INTEGER, INTENT(in) :: nparts IF (.NOT. ALLOCATED(p%Z) ) THEN p%Nptot = nparts ALLOCATE(p%Z(nparts)) ALLOCATE(p%R(nparts)) ALLOCATE(p%THET(nparts)) ALLOCATE(p%WZ(nparts)) ALLOCATE(p%WR(nparts)) ALLOCATE(p%UR(nparts)) ALLOCATE(p%UZ(nparts)) ALLOCATE(p%UTHET(nparts)) ALLOCATE(p%URold(nparts)) ALLOCATE(p%UZold(nparts)) ALLOCATE(p%UTHETold(nparts)) ALLOCATE(p%Gamma(nparts)) ALLOCATE(p%Rindex(nparts)) ALLOCATE(p%Zindex(nparts)) ALLOCATE(p%partindex(nparts)) ALLOCATE(p%pot(nparts)) ALLOCATE(p%Er(nparts)) ALLOCATE(p%Ez(nparts)) ALLOCATE(p%GAMMAold(nparts)) p%partindex=-1 p%URold=0 p%UZold=0 p%UTHETold=0 p%rindex=0 p%zindex=0 p%WR=0 p%WZ=0 p%UR=0 p%UZ=0 p%UTHET=0 p%Z=0 p%R=0 p%THET=0 p%Gamma=1 p%Er=0 p%Ez=0 p%pot=0 p%gammaold=1 END IF END SUBROUTINE creat_parts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Loads the particles at the beginning of the simulation and create the parts variable if necessary !--------------------------------------------------------------------------- SUBROUTINE load_parts USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical USE mpi INTEGER:: i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) parts%Nptot=nplasma ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V and uniform in R CALL loaduniformRZ(VR, VZ, VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 CALL loadDavidson(VR, VZ, VTHET, lodunir) CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2 CALL loadDavidson(VR, VZ, VTHET, lodinvr) CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R CALL loadDavidson(VR, VZ, VTHET, lodgausr) CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp CALL loadDavidson(VR, VZ, VTHET, lodunir) CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz @brief Checks for each particle if the z position is outside of the local/global simulation space. !> Depending on the boundary conditions, the leaving particles are sent to the correct neighbouring MPI process !> or deleted. ! !> @author Guillaume Le Bars EPFL/SPC !--------------------------------------------------------------------------- SUBROUTINE bound(p) USE basic, ONLY: zgrid, rnorm, nz, Zbounds, cstep, mpirank, partperiodic, step, leftproc, rightproc + IMPLICIT NONE type(particles), INTENT(INOUT):: p INTEGER :: i, rsendnbparts, lsendnbparts, nblostparts INTEGER, DIMENSION(p%Nptot) :: sendhole INTEGER, DIMENSION(p%Nptot) :: losthole LOGICAL:: leftcomm, rightcomm rsendnbparts=0 lsendnbparts=0 nblostparts=0 losthole=0 sendhole=0 ! We communicate with the left processus leftcomm = leftproc .ne. -1 .or. partperiodic ! We communicate with the right processus rightcomm = rightproc .ne. -1 .or. partperiodic IF (p%Nptot .NE. 0) THEN ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,p%Nptot IF(step .ne. 0) THEN ! If the particle is to the right of the local simulation space, it is sent to the right MPI process IF (p%Z(i) .ge. zgrid(Zbounds(mpirank+1))) THEN !$OMP CRITICAL (nbparts) IF(rightcomm) THEN rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i ELSE !WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in z" !WRITE(*,*) "!!!!!!!!!! Particle removed !!!!!!!!!!" nblostparts=nblostparts+1 losthole(nblostparts)=i END IF !$OMP END CRITICAL (nbparts) ! If the particle is to the left of the local simulation space, it is sent to the left MPI process ELSE IF (p%Z(i) .lt. zgrid(Zbounds(mpirank))) THEN !$OMP CRITICAL (nbparts) IF(leftcomm) THEN lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i ELSE !WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in z" !WRITE(*,*) "!!!!!!!!!! Particle removed !!!!!!!!!!" nblostparts=nblostparts+1 losthole(nblostparts)=i END IF !$OMP END CRITICAL (nbparts) END IF END IF ! The periodic boundary conditions in z are here applied DO WHILE (p%Z(i) .GT. zgrid(nz)) p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0) END DO DO WHILE (p%Z(i) .LT. zgrid(0)) p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0) END DO END DO !$OMP END PARALLEL DO END IF IF(mpisize .gt. 1 .and. step .ne. 0) THEN ! We send the particles leaving the local simulation space to the closest neighbour CALL particlescommunication(p, lsendnbparts, rsendnbparts, sendhole) END IF ! If the boundary conditions are not periodic, we delete the corresponding particles IF(nblostparts .gt. 0) THEN DO i=nblostparts,1,-1 !WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in z: ", p%Z(losthole(i))*rnorm, zgrid(0)*rnorm, zgrid(nz)*rnorm CALL delete_part(p, losthole(i)) END DO WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank END IF END subroutine bound !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Compute the grid cell indices for each particle as well as the distance weight Wr, Wz. !--------------------------------------------------------------------------- SUBROUTINE localisation(p) USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm type(particles), INTENT(INOUT):: p INTEGER :: i, nblostparts INTEGER, DIMENSION(p%Nptot) :: losthole REAL(kind=db):: invdz, invdr(3) losthole=0 nblostparts=0 invdz=1/dz invdr=1/dr IF (p%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,p%Nptot p%zindex(i)=(p%Z(i)-zgrid(0))*invdz IF (p%R(i) .GT. rgrid(0) .AND. p%R(i) .LT. rgrid(nnr(1))) THEN p%rindex(i)=(p%R(i)-rgrid(0))*invdr(1) ELSE IF(p%R(i) .GE. rgrid(nnr(1)) .AND. p%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN p%rindex(i)=(p%R(i)-rgrid(nnr(1)))*invdr(2)+nnr(1) ELSE IF(p%R(i) .GE. rgrid(nnr(1)+nnr(2)) .AND. p%R(i) .LT. rgrid(nr)) THEN p%rindex(i)=(p%R(i)-rgrid(nnr(1)+nnr(2)))*invdr(3)+nnr(1)+nnr(2) ELSE ! If the particle is outside of the simulation space in the r direction, it is deleted. !$OMP CRITICAL (lostparts) nblostparts=nblostparts+1 losthole(nblostparts)=i !$OMP END CRITICAL (lostparts) END IF END DO !$OMP END PARALLEL DO IF(nblostparts.gt.0) THEN DO i=nblostparts,1,-1 !WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in r: ", p%R(losthole(i))*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm CALL delete_part(p,losthole(i)) END DO WRITE(*,'(i6.2,a,i6.2)') nblostparts, " particles lost in r on process: ", mpirank END IF !$OMP PARALLEL DO SIMD DEFAULT(SHARED) DO i=1,p%Nptot p%WZ(i)=(p%Z(i)-zgrid(p%zindex(i)))/dz; p%WR(i)=(p%R(i)-rgrid(p%rindex(i)))/(rgrid(p%rindex(i)+1)-rgrid(p%rindex(i))); END DO !$OMP END PARALLEL DO SIMD END IF END SUBROUTINE localisation !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief General routine to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle. ! !--------------------------------------------------------------------------- SUBROUTINE comp_velocity(p) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : nlclassical type(particles), INTENT(INOUT):: p ! Store old Velocities CALL swappointer(p%UZold, p%UZ) CALL swappointer(p%URold, p%UR) CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (nlclassical) THEN CALL comp_velocity_fun(p, gamma_classical) ELSE CALL comp_velocity_fun(p, gamma_relativistic) END IF END SUBROUTINE comp_velocity !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine called by comp_velocity to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using the routine computing gamma as an input. This avoid the nlclassical flag check on each particle. ! !> @param[in] gamma the function used to compute the value of the lorentz factor \f$\gamma\f$ !--------------------------------------------------------------------------- SUBROUTINE comp_velocity_fun(p, gamma) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz interface - subroutine gamma(gam, UZ, UR, UTHET) + ELEMENTAL subroutine gamma(gam, UZ, UR, UTHET) USE constants REAL(kind=db), INTENT(IN):: UR,UZ,UTHET REAL(kind=db), INTENT(OUT):: gam end subroutine end interface type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau REAL(kind=db):: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2 INTEGER:: J1, J2, J3, J4 INTEGER:: i ! Normalized time increment tau=omegac/2/omegap*dt/tnorm IF (p%Nptot .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,p%Nptot J1=(p%rindex(i))*(nz+1) + p%zindex(i)+1 J2=(p%rindex(i))*(nz+1) + p%zindex(i)+2 J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1 J4=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+2 ! Interpolation for magnetic field BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) & & +p%WZ(i)*(1-p%WR(i))*Bz(J3) & & +(1-p%WZ(i))*p%WR(i)*Bz(J2) & & +p%WZ(i)*p%WR(i)*Bz(J1) BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) & & +p%WZ(i)*(1-p%WR(i))*Br(J3) & & +(1-p%WZ(i))*p%WR(i)*Br(J2) & & +p%WZ(i)*p%WR(i)*Br(J1) ! First half of electric pulse p%UZ(i)=p%UZold(i)+p%Ez(i)*tau p%UR(i)=p%URold(i)+p%ER(i)*tau CALL gamma(p%Gamma(i), p%UZ(i), p%UR(i), p%UTHETold(i)) ! Rotation along magnetic field ZBZ=tau*BRZ/p%Gamma(i) ZBR=tau*BRR/p%Gamma(i) ZPZ=p%UZ(i)-ZBR*p%UTHETold(i) !u'_{z} ZPR=p%UR(i)+ZBZ*p%UTHETold(i) !u'_{r} ZPTHET=p%UTHETold(i)+(ZBR*p%UZ(i)-ZBZ*p%UR(i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR p%UZ(i)=p%UZ(i)-ZBR2*ZPTHET !u+_{z} p%UR(i)=p%UR(i)+ZBZ2*ZPTHET !u+_{r} p%UTHET(i)=p%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration p%UZ(i)=p%UZ(i)+p%EZ(i)*tau p%UR(i)=p%UR(i)+p%ER(i)*tau ! Final computation of the Lorentz factor CALL gamma(p%Gamma(i), p%UZ(i), p%UR(i), p%UTHETold(i)) END DO !$OMP END PARALLEL DO SIMD END IF collected=.false. END SUBROUTINE comp_velocity_fun !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the classical simulations. !> This routine systematically returns 1.0 to treat the system according to classical dynamic. ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- ELEMENTAL SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET) !$OMP declare simd(gamma_classical) REAL(kind=db), INTENT(IN):: UR,UZ,UTHET REAL(kind=db), INTENT(OUT):: gamma gamma=1.0 END SUBROUTINE gamma_classical !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the relativistic simulations. !> This routine computes the Lorentz factor \f$\gamma=\sqrt{1+\mathbf{\gamma\beta}^2}\f$ ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- ELEMENTAL SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET) !$OMP declare simd(gamma_relativistic) REAL(kind=db), INTENT(IN):: UR,UZ,UTHET REAL(kind=db), INTENT(OUT):: gamma gamma=sqrt(1+UZ**2+UR**2+UTHET**2) END SUBROUTINE !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles position at time t+1 !> This routine computes the particles position at time t+1 according to the Bunemann algorithm. !--------------------------------------------------------------------------- SUBROUTINE push(p) Use basic, ONLY: dt, tnorm type(particles), INTENT(INOUT):: p REAL(kind=db):: XP, YP, COSA, SINA, U1, U2, ALPHA INTEGER :: i IF (p%Nptot .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(XP, YP, COSA, SINA, U1, U2, ALPHA) DO i=1,p%Nptot ! Local Cartesian coordinates XP=p%R(i)+dt/tnorm*p%UR(i)/p%Gamma(i) YP=dt/tnorm*p%UTHET(i)/p%Gamma(i) ! Conversion to cylindrical coordiantes p%Z(i)=p%Z(i)+dt/tnorm*p%UZ(i)/p%Gamma(i) p%R(i)=sqrt(XP**2+YP**2) ! Computation of the rotation angle IF (p%R(i) .EQ. 0) THEN COSA=1 SINA=0 ALPHA=0 ELSE COSA=XP/p%R(i) SINA=YP/p%R(i) ALPHA=asin(SINA) END IF ! New azimuthal position p%THET(i)=MOD(p%THET(i)+ALPHA,2*pi) ! Velocity in rotated reference frame U1=COSA*p%UR(i)+SINA*p%UTHET(i) U2=-SINA*p%UR(i)+COSA*p%UTHET(i) p%UR(i)=U1 p%UTHET(i)=U2 END DO !$OMP END PARALLEL DO SIMD END IF collected=.false. END SUBROUTINE push !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes several diagnostic quantities !> This routine computes the total kinetic and electric potential energy. It also computes the first and second order !> moments of the distribution function ! !--------------------------------------------------------------------------- SUBROUTINE diagnostics ! ! Compute energies ! USE constants, ONLY: vlight - USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, ierr, step, it0d, it2d, nlend, nr, itparts, nplasma, partperiodic, newres, Zbounds + USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, ierr, step, it0d, it2d, nlend, nr,& + & itparts, nplasma, partperiodic, newres, Zbounds INTEGER :: i, gridindex INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE) CALL MPI_WAITALL(3, ireducerequest, stat, ierr) ! Reset the quantities ekin=0 epot=0 etot=0 ! Computation of the kinetic and potential energy as well as fluid velocities and density !$OMP PARALLEL DO SIMD REDUCTION(+:epot, ekin) DEFAULT(SHARED) DO i=1,parts%Nptot ! Potential energy epot=epot+0.5*qsim*parts%pot(i)*phinorm ! Kinetic energy IF(.not. nlclassical) THEN ekin=ekin+msim*vlight**2*(parts%Gamma(i)-1) ELSE - ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i))) + ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))& + & +abs(parts%UZold(i)*parts%UZ(i))& + & +abs(parts%UTHETold(i)*parts%UTHET(i))) END IF END DO !$OMP END PARALLEL DO SIMD ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ.1 .or. (newres .AND. step .EQ. 1)) THEN ! Compute the local total energy loc_etot0 = epot+ekin END IF etot=loc_etot0 ! The computed energy is sent to the root process IF(mpisize .gt.1) THEN IF(mpirank .eq.0 ) THEN - CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ireducerequest(1), ierr) - CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ireducerequest(2), ierr) - CALL MPI_IREDUCE(etot, etot0, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + CALL MPI_IREDUCE(etot, etot0, 1, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) ELSE - CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + CALL MPI_IREDUCE(epot, epot, 1, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ireducerequest(1), ierr) - CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + CALL MPI_IREDUCE(ekin, ekin, 1, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ireducerequest(2), ierr) - CALL MPI_IREDUCE(etot, etot0, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + CALL MPI_IREDUCE(etot, etot0, 1, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) END IF ELSE etot0=loc_etot0 END IF ! Wait for the communication of the energy to be finished CALL MPI_WAITALL(3, ireducerequest(1:3), stat(1:3*MPI_STATUS_SIZE), ierr) !CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr) ! Compute the total energy etot=epot+ekin ! Send the local number of particles on each node to the root process IF(modulo(step,it0d).eq. 0 .and. mpisize .gt. 1) THEN Nplocs_all(mpirank)=parts%Nptot IF(mpirank .eq.0 ) THEN CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) END IF IF(mpirank .eq. 0 .AND. partperiodic .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN !WRITE(*,*) "Error particle lost on some processus" !CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) END IF END IF !Calls the communications to send the particles data to the root process for diagnostic file every itparts time steps IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN CALL collectparts(parts) END IF end subroutine diagnostics !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Collect the particles positions and velocities on the root process. !> If the collection has already been performed at this time step, the routine does nothing. ! !--------------------------------------------------------------------------- SUBROUTINE collectparts(p) USE basic, ONLY: mpirank, mpisize, ierr type(particles), INTENT(INOUT):: p INTEGER, DIMENSION(:), ALLOCATABLE :: displs INTEGER:: i IF(collected) RETURN ! exit subroutine if particles have already been collected during this time step IF(mpirank .ne. 0) THEN ALLOCATE(displs(0:mpisize-1)) displs=0 ! Send Particles informations to root process - CALL MPI_Gatherv(p%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%Z, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(p%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%R, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(p%THET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%THET, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(p%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UR, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(p%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UZ, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(p%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UTHET, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(p%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%pot, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%Z, Nplocs_all(mpirank), db_type, p%Z, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%R, Nplocs_all(mpirank), db_type, p%R, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%THET, Nplocs_all(mpirank), db_type, p%THET, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%UR, Nplocs_all(mpirank), db_type, p%UR, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%UZ, Nplocs_all(mpirank), db_type, p%UZ, Nplocs_all, displs, & + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%UTHET, Nplocs_all(mpirank), db_type, p%UTHET, Nplocs_all, displs, & + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(p%pot, Nplocs_all(mpirank), db_type, p%pot, Nplocs_all, displs, & + & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%Rindex, Nplocs_all(mpirank), MPI_INTEGER, p%Rindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%Zindex, Nplocs_all(mpirank), MPI_INTEGER, p%Zindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%partindex, Nplocs_all(mpirank), MPI_INTEGER, p%partindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) ELSE ALLOCATE(displs(0:mpisize-1)) displs(0)=0 Do i=1,mpisize-1 displs(i)=displs(i-1)+Nplocs_all(i-1) END DO ! Receive particle information from all processes - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%Z, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%R, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%THET, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UR, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UZ, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%UTHET, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, p%pot, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%Z, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%R, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%THET, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%UR, Nplocs_all, displs,& + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%UZ, Nplocs_all, displs, & + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%UTHET, Nplocs_all, displs, & + & db_type, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), db_type, p%pot, Nplocs_all, displs, & + & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%Rindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%Zindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, p%partindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) p%partindex(sum(Nplocs_all)+1:)=-1 END IF collected=.TRUE. END SUBROUTINE collectparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Computes the velocities at time t-1/2 to keep the second order precision in time on the velocity. ! !--------------------------------------------------------------------------- SUBROUTINE adapt_vinit(p) !! Computes the velocity at time -dt/2 from velocities computed at time 0 ! USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, ZBZ2, ZBR2, Vperp, v2 INTEGER :: J1, J2, J3, J4, i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ! In case Davidson distribution is used the longitudinal and radial velocities are adapted to take into account the ! electric potential. IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4) THEN ALLOCATE(VR(p%Nptot),VZ(p%Nptot),VTHET(p%Nptot)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=p%UTHET/p%Gamma*vnorm DO i=1,p%Nptot Vperp=sqrt(MAX(2*H0/me+2*elchar/me*p%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN p%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 p%Gamma(i)=sqrt(1/(1-v2/vnorm**2)) END IF p%UR(i)=p%Gamma(i)*VR(i)/vnorm p%UZ(i)=p%Gamma(i)*VZ(i)/vnorm p%UTHET(i)=p%Gamma(i)*VTHET(i)/vnorm END DO DEALLOCATE(VR,VZ,VTHET) END IF ! Normalised time increment tau=-omegac/2/omegap*dt/tnorm ! Store old Velocities CALL swappointer(p%UZold, p%UZ) CALL swappointer(p%URold, p%UR) CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (p%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,p%Nptot ! Compute the particle linear indices for the magnetic field interpolation J1=(p%rindex(i))*(nz+1)+p%zindex(i)+1 J2=J1+1 J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1 J4=J3+1 ! Interpolation for magnetic field BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) & & +p%WZ(i)*(1-p%WR(i))*Bz(J3) & & +(1-p%WZ(i))*p%WR(i)*Bz(J2) & & +p%WZ(i)*p%WR(i)*Bz(J1) BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) & & +p%WZ(i)*(1-p%WR(i))*Br(J3) & & +(1-p%WZ(i))*p%WR(i)*Br(J2) & & +p%WZ(i)*p%WR(i)*Br(J1) ! Half inverse Rotation along magnetic field ZBZ=tau*BRZ/p%Gammaold(i) ZBR=tau*BRR/p%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(p%UZold(i)-ZBR*p%UTHETold(i))/SQR !u-_{z} ZPR=(p%URold(i)+ZBZ*p%UTHETold(i))/SQR !u-_{r} ZPTHET=p%UTHETold(i)+(ZBR*p%UZold(i)-ZBZ*p%URold(i))/SQR !u-_{theta} p%UZ(i)=ZPZ p%UR(i)=ZPR p%UTHET(i)=ZPTHET ! half of decceleration p%UZ(i)=p%UZ(i)+p%Ez(i)*tau p%UR(i)=p%UR(i)+p%Er(i)*tau IF(.not. nlclassical) THEN p%Gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2) END IF END DO !$OMP END PARALLEL DO END IF END SUBROUTINE adapt_vinit !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief In the case of MPI parallelism, computes the indices of the particle assigned to the current process. !--------------------------------------------------------------------------- SUBROUTINE distribpartsonprocs ! Computes the start and end indices for the Z boundaries on local processus ! Computes the particle indices from initial particle loading vector, that stay in current process USE basic, ONLY: nz, ierr, Zbounds INTEGER:: k, i, nbparts REAL(kind=db):: idealnbpartsperproc INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k) INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector INTEGER:: partindexlast ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector INTEGER, DIMENSION(0:nz):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1) INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction INTEGER:: Zperproc, zindex partspercol=0 idealnbpartsperproc = FLOOR(REAL(parts%Nptot)/REAL(mpisize)) Zmin=MINVAL(parts%Zindex(1:parts%Nptot)) Zmax=MAXVAL(parts%Zindex(1:parts%Nptot)) Zperproc=(Zmax-Zmin)/mpisize IF(Zmax .eq. 0) Zmax=nz DO k=1,parts%Nptot zindex=parts%Zindex(k) partspercol(zindex)=partspercol(zindex)+1 END DO IF (Zperproc .eq. 0) THEN !! No particles are present initially Zperproc=nz/mpisize Zmin=0 DO k=1,mpisize-1 IF(k .lt. mpisize-1-MODULO(Zmax-Zmin,mpisize)) THEN Zbounds(k)=Zmin+k*Zperproc-1 ELSE Zbounds(k)=Zmin+k*Zperproc-1+k-mpisize+2+MODULO(Zmax-Zmin,mpisize) END IF END DO ELSE i=0 DO k=1,mpisize-1 nbparts=0 DO WHILE(nbparts<0.99*idealnbpartsperproc .and. i .lt. Zmax) nbparts=nbparts+partspercol(i) i=i+1 END DO Zbounds(k)=i END DO END IF partindexstart=1 totparts=0 partindexlast=parts%Nptot !DO k=0,mpisize-2 ! Nplocs_all(k)=idealnbpartsperproc !END DO ! NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize) !DO k=1, mpisize-1 ! Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc)) !END DO DO k=0,mpisize-1 Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1)) END DO ! Store local end of particle vector partindexlast=SUM(Nplocs_all(0:mpirank)) ! Store local start of particle vector IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1 IF(mpirank .ne. 0 .and. Nplocs_all(mpirank) .gt. 0) THEN DO k= partindexstart, partindexlast CALL move_part(parts,k,k-partindexstart+1) END DO END IF IF(INT(SUM(Nplocs_all)) .ne. parts%Nptot) THEN WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " Nptot:", parts%Nptot CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF parts%Nptot=Nplocs_all(mpirank) WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot END SUBROUTINE distribpartsonprocs !_______________________________________________________________________________ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Manage the particle communication between neighbours. !> This routine is responsible to receive the incoming particles from the MPI neighbours and to send its outgoing !> particles to these neighbours ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative !--------------------------------------------------------------------------- SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, sendholes) USE mpihelper, ONLY: particle_type USE basic, ONLY: rightproc, leftproc, ierr type(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(:) INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0 INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE) INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL ! Send and receive the number of particles to exchange CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, recvrequest(1), ierr) CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, recvrequest(2), ierr) CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, sendrequest(1), ierr) CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, sendrequest(2), ierr) CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr) CALL MPI_Wait(recvrequest(2), recvstatus(2), ierr) recvrequest=MPI_REQUEST_NULL lreceivednbparts=lrecvnbparts rreceivednbparts=rrecvnbparts ! Re/allocate enough memory to store the incoming particles IF( ALLOCATED(rrecvpartbuff) .AND. (size(rrecvpartbuff) .lt. rrecvnbparts) ) DEALLOCATE(rrecvpartbuff) IF(.not. ALLOCATED(rrecvpartbuff) .and. rrecvnbparts .gt. 0) ALLOCATE(rrecvpartbuff(INT(rreceivednbparts*1.3))) IF( ALLOCATED(lrecvpartbuff) .AND. (size(lrecvpartbuff) .lt. lrecvnbparts) ) DEALLOCATE(lrecvpartbuff) IF((.not. ALLOCATED(lrecvpartbuff) ).and. lrecvnbparts .gt. 0) ALLOCATE(lrecvpartbuff(INT(lreceivednbparts*1.3))) ! Receive particles from left and right processes to the corresponding buffers IF ( lrecvnbparts .gt. 0) THEN CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, recvrequest(1), ierr) END IF IF( rrecvnbparts .gt. 0) THEN CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, recvrequest(2), ierr) END IF ! Copy the leaving particles to the corresponding send buffers IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes) END IF CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr) CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr) ! Send the particles to the left and right neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, sendrequest(1), ierr) END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, sendrequest(2), ierr) END IF ! Receive the incoming parts in the receive buffers IF ( lreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(1) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF END IF IF ( rreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(2), recvstatus(MPI_STATUS_SIZE+1), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(2) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF END IF ! Copy the incoming particles from the receive buffers to the simulation parts variable CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes) ! Wait for the outgoing particles to be fully received by the neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr) END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(2), sendstatus(2), ierr) END IF ! ! END SUBROUTINE particlescommunication !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the receive buffers to the local simulation variable parts. !> The incoming particles will first be stored in the holes left by the outgoing particles, then they !> will be added at the end of the parts variable ! !> @param [in] rrecvnbparts number of particles received from the right neighbour (mpirank+1) !> @param [in] lrecvnbparts number of particles received from the left neighbour (mpirank-1) !> @param [in] sendnbparts total number of particles having left the local domain !> @param [in] sendholes array containing the indices of the particle having left the local domain in ascending order. !--------------------------------------------------------------------------- SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts, sendholes) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(:) INTEGER k,partpos, partdiff ! computes if we received less particles than we sent partdiff=sendnbparts-rrecvnbparts-lrecvnbparts ! First import the particles coming from the right IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN ! Fill the holes partpos=abs(sendholes(k)) ELSE ! Add at the end of parts p%Nptot=p%Nptot+1 partpos=p%Nptot END IF CALL Insertincomingpart(p, rrecvpartbuff, k, partpos) END DO END IF ! Then import the particles coming from the left IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN partpos=abs(sendholes(k+rrecvnbparts)) ELSE p%Nptot=p%Nptot+1 partpos=p%Nptot END IF CALL Insertincomingpart(p, lrecvpartbuff, k, partpos) END DO END IF ! If we received less particles than we sent, we fill the remaining holes with the particles from the end of the ! parts arrays IF(partdiff .gt. 0) THEN DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts CALL move_part(parts, parts%Nptot, abs(sendholes(k))) p%Nptot = p%Nptot-1 END DO END IF ! END SUBROUTINE Addincomingparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy one particle from the receive buffers to the local simulation variable parts. ! !> @param [in] buffer receive buffer containing the particles parameters to copy from !> @param [in] bufferindex particle index in the receive buffer !> @param [in] partsindex destination particle index in the local parts variable !--------------------------------------------------------------------------- SUBROUTINE Insertincomingpart(p, buffer, bufferindex, partsindex) USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer p%partindex(partsindex) = buffer(bufferindex)%partindex p%Rindex(partsindex) = buffer(bufferindex)%Rindex p%Zindex(partsindex) = buffer(bufferindex)%Zindex p%R(partsindex) = buffer(bufferindex)%R p%Z(partsindex) = buffer(bufferindex)%Z p%THET(partsindex) = buffer(bufferindex)%THET p%UZ(partsindex) = buffer(bufferindex)%UZ p%UR(partsindex) = buffer(bufferindex)%UR p%UTHET(partsindex) = buffer(bufferindex)%UTHET p%Gamma(partsindex) = buffer(bufferindex)%Gamma ! ! END SUBROUTINE Insertincomingpart !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy one particle from the local parts variable to the send buffer. ! !> @param [in] buffer send buffer to copy to !> @param [in] bufferindex particle index in the send buffer !> @param [in] partsindex origin particle index in the local parts variable !--------------------------------------------------------------------------- SUBROUTINE Insertsentpart(p, buffer, bufferindex, partsindex) USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer buffer(bufferindex)%partindex = p%partindex(partsindex) buffer(bufferindex)%Rindex = p%Rindex(partsindex) buffer(bufferindex)%Zindex = p%Zindex(partsindex) buffer(bufferindex)%R = p%R(partsindex) buffer(bufferindex)%Z = p%Z(partsindex) buffer(bufferindex)%THET = p%THET(partsindex) buffer(bufferindex)%UZ = p%UZ(partsindex) buffer(bufferindex)%UR = p%UR(partsindex) buffer(bufferindex)%UTHET = p%UTHET(partsindex) buffer(bufferindex)%Gamma = p%Gamma(partsindex) ! END SUBROUTINE Insertsentpart !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the local parts variable to the left and right send buffers. ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative !--------------------------------------------------------------------------- SUBROUTINE AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(:) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 ! Verify if the buffers are big enough to store the outgoing particles IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff) IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff) ! Allocate the buffers if needed IF(.not. ALLOCATED(rsendpartbuff) .and. rsendnbparts .gt. 0) ALLOCATE(rsendpartbuff(INT(rsendnbparts*1.3))) IF(.not. ALLOCATED(lsendpartbuff) .and. lsendnbparts .gt. 0) ALLOCATE(lsendpartbuff(INT(lsendnbparts*1.3))) ! Loop over the outgoing particles and fill the correct send buffer Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(sendholes(k)) IF(sendholes(k) .GT. 0) THEN rsendpos=rsendpos+1 CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos) ELSE IF(sendholes(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(p, lsendpartbuff, lsendpos, partpos) END IF END DO ! ! END SUBROUTINE AddPartSendBuffers !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Exchange two particles in the parts variable. ! !> @param [in] index1 index in parts of the first particle to exchange. !> @param [in] index2 index in parts of the second particle to exchange. !--------------------------------------------------------------------------- SUBROUTINE exchange_parts(p, index1, index2) TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index1, index2 REAL(kind=db):: R, Z, THET, UR, UZ, UTHET, Gamma INTEGER :: Rindex, Zindex, partindex !! Exchange particle at index1 with particle at index2 ! Store part at index1 in temporary value partindex= p%partindex(index1) Gamma = p%Gamma(index1) R = p%R(index1) Z = p%Z(index1) THET = p%THET(index1) UR = p%UR(index1) UTHET = p%UTHET(index1) UZ = p%UZ(index1) Rindex = p%Rindex(index1) Zindex = p%Zindex(index1) ! Move part at index2 in part at index 1 p%partindex(index1)= p%partindex(index2) p%Gamma(index1) = p%Gamma(index2) p%R(index1) = p%R(index2) p%Z(index1) = p%Z(index2) p%THET(index1) = p%THET(index2) p%UR(index1) = p%UR(index2) p%UTHET(index1) = p%UTHET(index2) p%UZ(index1) = p%UZ(index2) p%Rindex(index1) = p%Rindex(index2) p%Zindex(index1) = p%Zindex(index2) ! Move temporary values from part(index1) to part(index2) p%partindex(index2)= partindex p%Gamma(index2) = Gamma p%R(index2) = R p%Z(index2) = Z p%THET(index2) = THET p%UR(index2) = UR p%UTHET(index2) = UTHET p%UZ(index2) = UZ p%Rindex(index2) = Rindex p%Zindex(index2) = Zindex END SUBROUTINE exchange_parts SUBROUTINE change_parts_allocation(p, sizedifference) implicit none TYPE(particles), INTENT(INOUT):: p INTEGER,INTENT(IN) :: sizedifference CALL change_array_size_int(p%Rindex, sizedifference) CALL change_array_size_int(p%Zindex, sizedifference) CALL change_array_size_int(p%partindex, sizedifference) CALL change_array_size_dp(p%ER,sizedifference) CALL change_array_size_dp(p%EZ,sizedifference) CALL change_array_size_dp(p%pot,sizedifference) CALL change_array_size_dp(p%R,sizedifference) CALL change_array_size_dp(p%Z,sizedifference) CALL change_array_size_dp(p%THET,sizedifference) CALL change_array_size_dp(p%WR,sizedifference) CALL change_array_size_dp(p%WZ,sizedifference) CALL change_array_size_dp_ptr(p%UR,sizedifference) CALL change_array_size_dp_ptr(p%URold,sizedifference) CALL change_array_size_dp_ptr(p%UZ,sizedifference) CALL change_array_size_dp_ptr(p%UZold,sizedifference) CALL change_array_size_dp_ptr(p%UTHET,sizedifference) CALL change_array_size_dp_ptr(p%UTHETold,sizedifference) CALL change_array_size_dp_ptr(p%Gamma,sizedifference) CALL change_array_size_dp_ptr(p%Gammaold,sizedifference) p%Nptot=MIN(p%Nptot,size(p%R)+sizedifference) END SUBROUTINE change_parts_allocation SUBROUTINE change_array_size_dp(arr, sizedifference) implicit none REAL(kind=db), ALLOCATABLE, INTENT(INOUT):: arr(:) INTEGER, INTENT(IN):: sizedifference REAL(kind=db), ALLOCATABLE:: temp(:) INTEGER:: current_size, new_size if(allocated(arr)) THEN current_size=size(arr) new_size=current_size+sizedifference ALLOCATE(temp(new_size)) temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size)) DEALLOCATE(arr) CALL move_alloc(temp, arr) END IF END SUBROUTINE change_array_size_dp SUBROUTINE change_array_size_dp_ptr(arr, sizedifference) implicit none REAL(kind=db), POINTER, INTENT(INOUT):: arr(:) INTEGER, INTENT(IN):: sizedifference REAL(kind=db), POINTER:: temp(:) INTEGER:: current_size, new_size if(associated(arr)) THEN current_size=size(arr) new_size=current_size+sizedifference ALLOCATE(temp(new_size)) temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size)) DEALLOCATE(arr) arr=> temp END IF END SUBROUTINE change_array_size_dp_ptr SUBROUTINE change_array_size_int(arr, sizedifference) implicit none INTEGER, ALLOCATABLE, INTENT(INOUT):: arr(:) INTEGER, INTENT(IN):: sizedifference INTEGER, ALLOCATABLE:: temp(:) INTEGER:: current_size, new_size if(allocated(arr)) THEN current_size=size(arr) new_size=current_size+sizedifference ALLOCATE(temp(new_size)) temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size)) DEALLOCATE(arr) CALL move_alloc(temp,arr) END IF END SUBROUTINE change_array_size_int SUBROUTINE add_created_part(p, buffer) USE bsplines USE basic, ONLY: splrz, qsim, phinorm, nlclassical, msim IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:) INTEGER:: i, nptotinit, parts_size_increase, nb_added nptotinit=p%Nptot+1 nb_added=size(buffer,1) ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nptot + nb_added .gt. size(p%Z,1)) THEN parts_size_increase=Max(floor(0.1*size(p%Z,1)),nb_added) CALL change_parts_allocation(p, parts_size_increase) END IF DO i=1,nb_added p%Nptot=p%Nptot+1 CALL Insertincomingpart(p, buffer,i,p%Nptot) END DO ! Compute the energy added by these new particles for energy conservation diagnostic IF(nb_added .gt. 0) THEN CALL getgrad(splrz, p%Z(nptotinit:p%Nptot), p%R(nptotinit:p%Nptot), & & p%pot(nptotinit:p%Nptot), p%EZ(nptotinit:p%Nptot), p%ER(nptotinit:p%Nptot)) p%EZ(nptotinit:p%Nptot)=-p%Ez(nptotinit:p%Nptot) p%ER(nptotinit:p%Nptot)=-p%ER(nptotinit:p%Nptot) loc_etot0=loc_etot0+qsim*sum(p%pot(nptotinit:p%Nptot))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0+sum(msim*vlight**2*(p%Gamma(nptotinit:p%Nptot)-1)) ELSE loc_etot0=loc_etot0+0.5*msim*vlight**2*sum(p%UR(nptotinit:p%Nptot)**2 & & +p%UZ(nptotinit:p%Nptot)**2 & & +p%UTHET(nptotinit:p%Nptot)**2) END IF END IF END SUBROUTINE !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Delete particle at given index removing its energy from the system ! !> @param [in] index index of particle to be deleted !--------------------------------------------------------------------------- SUBROUTINE delete_part(p, index) !! This will destroy particle at index USE constants, ONLY: vlight USE basic, ONLY: qsim, phinorm, msim, nlclassical TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index IF(index .le. p%Nptot) THEN loc_etot0=loc_etot0-qsim*p%pot(index)*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0-msim*vlight**2*(p%Gamma(index)-1) ELSE loc_etot0=loc_etot0-0.5*msim*vlight**2*(abs(p%URold(index)*p%UR(index))+abs(p%UZold(index)*p%UZ(index))+abs(p%UTHETold(index)*p%UTHET(index))) END IF ! We fill the gap CALL move_part(p, p%Nptot, index) p%partindex(p%Nptot)=-1 ! Reduce the total number of simulated parts p%Nptot=p%Nptot-1 END IF END SUBROUTINE delete_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Move particle with index sourceindex to particle with index destindex. !> !WARNING! This will overwrite particle at destindex. ! !> @param [in] sourceindex index in parts of the particle to move. !> @param [in] destindex index in parts of the moved particle destination. !--------------------------------------------------------------------------- SUBROUTINE move_part(p, sourceindex, destindex) !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex TYPE(particles), INTENT(INOUT)::p ! Move part at sourceindex in part at destindex p%partindex(destindex)= p%partindex(sourceindex) p%Gamma(destindex) = p%Gamma(sourceindex) p%R(destindex) = p%R(sourceindex) p%Z(destindex) = p%Z(sourceindex) p%THET(destindex) = p%THET(sourceindex) p%UR(destindex) = p%UR(sourceindex) p%UTHET(destindex) = p%UTHET(sourceindex) p%UZ(destindex) = p%UZ(sourceindex) p%Rindex(destindex) = p%Rindex(sourceindex) p%Zindex(destindex) = p%Zindex(sourceindex) END SUBROUTINE move_part !________________________________________________________________________________ SUBROUTINE destroy_parts(p) TYPE(particles) :: p p%nptot=0 IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z) IF(ALLOCATED(p%R)) DEALLOCATE(p%R) IF(ALLOCATED(p%THET)) DEALLOCATE(p%THET) IF(ALLOCATED(p%WZ)) DEALLOCATE(p%WZ) IF(ALLOCATED(p%WR)) DEALLOCATE(p%WR) IF(ASSOCIATED(p%UR)) DEALLOCATE(p%UR) IF(Associated(p%URold)) DEALLOCATE(p%URold) IF(Associated(p%UZ)) DEALLOCATE(p%UZ) IF(Associated(p%UZold)) DEALLOCATE(p%UZold) IF(Associated(p%UTHET)) DEALLOCATE(p%UTHET) IF(Associated(p%UTHETold)) DEALLOCATE(p%UTHETold) IF(Associated(p%Gamma)) DEALLOCATE(p%Gamma) IF(Associated(p%Gammaold)) DEALLOCATE(p%Gammaold) IF(ALLOCATED(p%Rindex)) DEALLOCATE(p%Rindex) IF(ALLOCATED(p%Zindex)) DEALLOCATE(p%Zindex) IF(ALLOCATED(p%partindex)) DEALLOCATE(p%partindex) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE clean_beam ! CALL destroy_parts(parts) ! END SUBROUTINE clean_beam !________________________________________________________________________________ SUBROUTINE swappointer( pointer1, pointer2) REAL(kind=db), DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2 REAL(kind=db), DIMENSION(:), POINTER:: temppointer temppointer=>pointer1 pointer1=>pointer2 pointer2=>temppointer END SUBROUTINE swappointer !_______________________________________________________________________________ SUBROUTINE loaduniformRZ(VR,VZ,VTHET) USE basic, ONLY: plasmadim, rnorm, temp, vnorm USE constants, ONLY: me, kb, elchar REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) REAL(kind=db):: vth ! Initial distribution in z with normalisation CALL loduni(1,parts%Z(1:parts%Nptot)) parts%Z(1:parts%Nptot)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(1:parts%Nptot))/rnorm ! Initial distribution in r with normalisation CALL loduni(2,parts%R(1:parts%Nptot)) parts%R(1:parts%Nptot)=(plasmadim(3)+parts%R(1:parts%Nptot)*(plasmadim(4)-plasmadim(3)))/rnorm ! Initial velocities distribution vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity CALL lodgaus(3,VZ) CALL lodgaus(5,VR) CALL lodgaus(7,VTHET) VZ=VZ*vth VR=VR*vth VTHET=VTHET*vth END SUBROUTINE loaduniformRZ !_______________________________________________________________________________ SUBROUTINE loadDavidson(VR,VZ,VTHET,lodr) USE constants, ONLY: me, kb, elchar USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, B0, width, qsim, msim, vnorm, & & omegac, zgrid, nz, rnorm, V, n0, nblock, temp interface subroutine lodr(nbase,y,ra,rb) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rb, ra end subroutine end interface REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb - REAL(kind=db) :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL, vth + REAL(kind=db) :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, VOL, vth INTEGER :: i, j, n, blockstart, blockend, addedpart, remainparts INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize Allocate(ra(nblock),rb(nblock)) !r0=(plasmadim(4)+plasmadim(3))/2 r0=sqrt(4*H0/(me*omegac**2)) halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) DO n=1,nblock ! Compute limits in radius and load radii for each part Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1) deltar2=1-MirrorRatio*cos(2*pi*Le/width) rb(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2)) ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/H0*deltar2)) END DO VOL=SUM(2*pi*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) qsim=VOL*n0*elchar/nplasma msim=abs(qsim)/elchar*me V=VOL/rnorm**3 - !H0=me*omegac**2*r0**2*0.5 - gamma=1/sqrt(1-omegac**2*r0**2/vlight**2) - !P0=H0/omegac blockstart=1 blockend=0 ALLOCATE(blocksize(nblock)) WRITE(*,*) "blocksize: ", size(blocksize), nblock DO n=1,nblock blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock END DO remainparts=parts%Nptot-SUM(blocksize) addedpart=1 n=nblock/2 j=1 DO WHILE(remainparts .GT. 0) blocksize(n)=blocksize(n)+addedpart remainparts=remainparts-addedpart n=n+j j=-1*(j+SIGN(1,j)) END DO DO n=1,nblock blockstart=blockend+1 blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot) ! Initial distribution in z with normalisation between magnetic mirrors CALL loduni(1,parts%Z(blockstart:blockend)) parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm CALL lodr(2,parts%R(blockstart:blockend),ra(n), rb(n)) parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm END DO IF(distribtype .eq. 5) THEN ! Initial velocities distribution vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity CALL lodgaus(3,VZ) CALL lodgaus(5,VR) CALL lodgaus(7,VTHET) VZ=VZ*vth/4 VR=VR*8*vth VTHET=VTHET*8*vth ELSE ! Load velocities theta velocity ! Loading of r and z velocity is done in adapt_vinit to have ! access to parts%pot DO i=1,parts%Nptot ! Interpolation for Magnetic potential rg=parts%R(i)*rnorm zg=(parts%Z(i)-halfLz)*rnorm Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width)) Pcomp=P0/rg/me Acomp=-SIGN(elchar/me*Athetpos,qsim) VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp) !VTHET(i)=Pcomp+Acomp END DO VTHET=VTHET/vnorm VZ=0._db VR=0._db END IF END SUBROUTINE loadDavidson !=============================================================================== SUBROUTINE Temp_rescale USE basic, ONLY: temprescale, nlclassical, nz ! INTEGER:: i, gridindex ! REAL(kind=db) :: vr, vz, vthet ! DO i=1,parts%Nptot ! gridindex=parts%zindex(i)+1+parts%rindex(i)*nz ! vr=parts%UR(i)/parts%Gamma(i)-fluidur(gridindex) ! vz=parts%UZ(i)/parts%Gamma(i)-fluiduz(gridindex) ! vthet=parts%UTHET(i)/parts%Gamma(i)-fluiduthet(gridindex) ! parts%UR(i)=vr*temprescale+fluidur(gridindex) ! parts%UTHET(i)=vthet*temprescale+fluiduthet(gridindex) ! parts%UZ(i)=vz*temprescale+fluiduz(gridindex) ! IF(nlclassical) THEN ! parts%Gamma(i)=1.0 ! ELSE ! parts%Gamma(i)=1/sqrt(1-(parts%UR(i)**2+parts%UTHET(i)**2+parts%UZ(i)**2)) ! END IF ! parts%UR(i)=parts%UR(i)*parts%Gamma(i) ! parts%UTHET(i)=parts%UTHET(i)*parts%Gamma(i) ! parts%UZ(i)=parts%UZ(i)*parts%Gamma(i) ! END DO END SUBROUTINE Temp_rescale !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Increase the number of macroparticles by separating each previous macroparticles into !> samplefactor new macroparticles of equally divided weight. The new sub particles are distributed !> uniformly in space to maintain the density and other moments. ! !> @param [in] samplefactor multiplicator of the number of macroparticles. !> @param [in] p particles type to increase. !--------------------------------------------------------------------------- SUBROUTINE upsample(p, samplefactor) USE basic, ONLY : nplasma, qsim, msim, dt, dr, dz INTEGER, INTENT(IN) ::samplefactor TYPE(particles), INTENT(INOUT):: p INTEGER:: i, j, currentindex REAL(kind=db), DIMENSION(p%Nptot) :: spreaddir ! random direction for the spread of each initial macro particle REAL(kind=db) :: dir ! Direction in which the particle is moved REAL(kind=db) :: dl ! Particle displacement used for ! Load and scale the direction angle for spreading the new particles CALL loduni(2, spreaddir) spreaddir=spreaddir*2*pi/samplefactor dl=min(dz,minval(dr))/100 DO i=1,p%Nptot DO j=1,samplefactor-1 currentindex=p%Nptot+(i-1)*(samplefactor-1)+j CALL move_part(p,i,currentindex) p%partindex(currentindex)=currentindex dir = spreaddir(i)+2*pi*j/samplefactor p%R(currentindex)=p%R(currentindex) + dl*cos(dir) p%Z(currentindex)=p%Z(currentindex) + dl*sin(dir) END DO p%partindex(i)=i p%R(i)=p%R(i) + dl*cos(spreaddir(i)) p%Z(i)=p%Z(i) + dl*sin(spreaddir(i)) END DO nplasma=nplasma*samplefactor qsim=qsim/samplefactor msim=msim/samplefactor p%Nptot=p%Nptot*samplefactor END SUBROUTINE upsample END MODULE beam diff --git a/src/chkrst.f90 b/src/chkrst.f90 index 0fe535f..ba8df48 100644 --- a/src/chkrst.f90 +++ b/src/chkrst.f90 @@ -1,144 +1,144 @@ SUBROUTINE chkrst(flag) ! ! Process checkpoint/restart file ! USE basic USE futils USE beam USE fields IMPLICIT NONE INTEGER, INTENT(in) :: flag INTEGER :: remainingparts REAL(kind=db):: old_msim, old_qsim, old_n0 INTEGER:: partsrank, partsdim(1) ! Only process 0 should save on file ! ! Local vars and arrays !________________________________________________________________________________ ! SELECT CASE(flag) !________________________________________________________________________________ ! 1. Open and read restart file ! CASE(0) IF(mpirank .ne. 0) RETURN CALL openf(rstfile, fidrst,'r','d') CALL getatt(fidrst, '/Basic', 'cstep', cstep) CALL getatt(fidrst, '/Basic', 'time', time) CALL getatt(fidrst, '/Basic', 'msim', old_msim) CALL getatt(fidrst, '/Basic', 'qsim', old_qsim) CALL getatt(fidrst, '/Basic', 'n0', old_n0) qsim=old_qsim/old_n0*n0 msim=old_msim/old_n0*n0 CALL getatt(fidrst, '/Basic', 'V', V) ! Renormalize V V=V*sqrt(n0/old_n0)**3 CALL getatt(fidrst, '/var0d', 'etot0', loc_etot0) CALL getatt(fidrst, '/var0d', 'epot', epot) CALL getatt(fidrst, '/var0d', 'ekin', ekin) CALL getatt(fidrst, '/var0d', 'etot', etot) CALL getatt(fidrst, '/Parts', 'nplasma', nplasma) CALL getatt(fidrst, '/Parts', 'remainingparts', remainingparts) IF (samplefactor .gt. 1 ) THEN ! We increase the number of macro particles CALL creat_parts(parts,remainingparts*samplefactor) ELSE CALL getdims(fidrst, '/Parts/Z', partsrank, partsdim) CALL creat_parts(parts,partsdim(1)) END IF CALL getarr(fidrst, '/Parts/Z', parts%Z) CALL getarr(fidrst, '/Parts/R', parts%R) ! Renormalize R and Z parts%R=parts%R*sqrt(n0/old_n0) parts%Z=parts%Z*sqrt(n0/old_n0) CALL getarr(fidrst, '/Parts/THET', parts%THET) CALL getarr(fidrst, '/Parts/UZ', parts%UZ) CALL getarr(fidrst, '/Parts/UR', parts%UR) CALL getarr(fidrst, '/Parts/UTHET', parts%UTHET) CALL getarr(fidrst, '/Parts/Zindex', parts%Zindex) CALL getarr(fidrst, '/Parts/Rindex', parts%Rindex) CALL getarr(fidrst, '/Parts/partindex', parts%partindex) CALL getatt(fidrst, '/Parts', 'remainingparts', parts%Nptot) WRITE(*,*) "Read ", parts%Nptot, " particles out of ", nplasma IF(isdataset(fidrst,'/Parts/fluidur')) THEN CALL getarr(fidrst, '/Parts/GAMMA', parts%Gamma) !CALL getarr(fidrst, '/Parts/fluidur', moments(2,:)) !CALL getarr(fidrst, '/Parts/fluiduthet', moments(3,:)) !CALL getarr(fidrst, '/Parts/fluiduz', moments(4,:)) IF(temprescale .gt. 0) THEN CALL Temp_rescale END IF END IF CALL closef(fidrst) IF(samplefactor .gt. 1) THEN ! We increase the number of macro particles CALL upsample(parts, samplefactor) END IF WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!" !________________________________________________________________________________ ! 2. Create and write to restart file (DP reals) ! CASE(1) IF(mpisize .gt. 1) THEN call collectparts(parts) END IF IF(mpirank .ne. 0) RETURN IF( .NOT. nlsave ) RETURN CALL mv2bk(rstfile) CALL creatf(rstfile, fidrst, real_prec='d', desc='Restart file') CALL creatg(fidrst, '/Basic', 'Basic data') CALL attach(fidrst, '/Basic', 'cstep', cstep) CALL attach(fidrst, '/Basic', 'time', time) CALL attach(fidrst, '/Basic', 'jobnum', jobnum) CALL attach(fidrst, '/Basic', 'qsim', qsim) CALL attach(fidrst, '/Basic', 'msim', msim) CALL attach(fidrst, '/Basic', 'n0', n0) CALL attach(fidrst, '/Basic', 'V', V) ! ! 0D variables ! CALL creatg(fidrst, '/var0d', '0D variables') CALL attach(fidrst, '/var0d','etot0', etot0) CALL attach(fidrst, '/var0d','epot', epot) CALL attach(fidrst, '/var0d','ekin', ekin) CALL attach(fidrst, '/var0d','etot', etot) ! ! Parts ! CALL creatg(fidrst, '/Parts', 'Particles data') CALL attach(fidrst, '/Parts', 'nplasma', nplasma) IF(mpisize .gt. 1) THEN remainingparts=sum(Nplocs_all) ELSE remainingparts=parts%Nptot END IF CALL attach(fidrst, '/Parts', 'remainingparts', remainingparts) - CALL putarr(fidrst, '/Parts/Z', parts%Z(1:remainingparts)) - CALL putarr(fidrst, '/Parts/R', parts%R(1:remainingparts)) - CALL putarr(fidrst, '/Parts/THET', parts%THET(1:remainingparts)) - CALL putarr(fidrst, '/Parts/UZ', parts%UZ(1:remainingparts)) - CALL putarr(fidrst, '/Parts/UR', parts%UR(1:remainingparts)) - CALL putarr(fidrst, '/Parts/UTHET', parts%UTHET(1:remainingparts)) - CALL putarr(fidrst, '/Parts/GAMMA', parts%Gamma(1:remainingparts)) - CALL putarr(fidrst, '/Parts/Zindex', parts%Zindex(1:remainingparts)) - CALL putarr(fidrst, '/Parts/Rindex', parts%Rindex(1:remainingparts)) - CALL putarr(fidrst, '/Parts/partindex', parts%partindex(1:remainingparts)) + CALL putarr(fidrst, '/Parts/Z', parts%Z) + CALL putarr(fidrst, '/Parts/R', parts%R) + CALL putarr(fidrst, '/Parts/THET', parts%THET) + CALL putarr(fidrst, '/Parts/UZ', parts%UZ) + CALL putarr(fidrst, '/Parts/UR', parts%UR) + CALL putarr(fidrst, '/Parts/UTHET', parts%UTHET) + CALL putarr(fidrst, '/Parts/GAMMA', parts%Gamma) + CALL putarr(fidrst, '/Parts/Zindex', parts%Zindex) + CALL putarr(fidrst, '/Parts/Rindex', parts%Rindex) + CALL putarr(fidrst, '/Parts/partindex', parts%partindex) CALL putarr(fidrst, '/Parts/fluidur', moments(2,:)) CALL putarr(fidrst, '/Parts/fluiduthet', moments(3,:)) CALL putarr(fidrst, '/Parts/fluiduz', moments(4,:)) ! ! Fields ! CALL closef(fidrst) WRITE(*,'(3x,a)') "Writing to restart file "//TRIM(rstfile)//" completed!" ! END SELECT ! END SUBROUTINE chkrst diff --git a/src/constants.f90 b/src/constants.f90 index 23d5478..a1f77fa 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -1,15 +1,16 @@ MODULE constants ! ! Define some constants ! - INTEGER, PARAMETER :: db = SELECTED_REAL_KIND(15,307) + INTEGER, PARAMETER :: dprequestedprec = 15 + INTEGER, PARAMETER :: db = SELECTED_REAL_KIND(dprequestedprec) ! REAL(kind=db), PARAMETER :: vlight = 299792458.0_db ! c REAL(kind=db), PARAMETER :: vacimp = 376.73031346177066_db ! \mu_0*c REAL(kind=db), PARAMETER :: eev = 510998.89613320108_db ! m_e*c^2/e REAL(kind=db), PARAMETER :: me = 9.109383E-31 ! electron mass [kg] REAL(kind=db), PARAMETER :: pi = 3.1415926535897931_db REAL(kind=db), PARAMETER :: elchar = 1.60217662E-19_db ! electron charge [C] REAL(kind=db), PARAMETER :: eps_0 = 8.85418781762E-12_db ! electric constant [F/m] REAL(kind=db), PARAMETER :: kb = 1.38064852E-23_db ! Boltzman constant [J/K] END MODULE constants diff --git a/src/diagnose.f90 b/src/diagnose.f90 index 57ae8d2..b58eeea 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,240 +1,243 @@ SUBROUTINE diagnose(kstep) ! ! Diagnostics ! USE basic USE futils USE hashtable USE beam, ONLY : parts, epot, ekin, etot, etot0, ireducerequest, Nplocs_all USE xg, ONLY : initw, updt_xg_var USE fields, ONLY: phi_spline USE mpi IMPLICIT NONE ! INTEGER, INTENT(in) :: kstep ! ! Local vars and arrays INTEGER, PARAMETER :: BUFSIZE = 20 CHARACTER(len=128) :: str, fname INTEGER:: Ntotal ! Total number of simulated particles remaining in the simulation INTEGER:: partsrank, partsdim(2) + INTEGER:: Nplocs_all_save(64) ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN !________________________________________________________________________________ ! 1. Initial diagnostics IF( kstep .EQ. 0) THEN ! WRITE(*,'(a)') ' Initial diagnostics' ! ! 1.1 Initial run or when NEWRES set to .TRUE. ! IF( .NOT. nlres .OR. newres) THEN CALL creatf(resfile, fidres, 'Espic2d result file', 'd') WRITE(*,'(3x,a,a)') TRIM(resfile), ' created' ! ! Label the run IF( LEN_TRIM(label1).GT.0 ) CALL attach(fidres, "/", "label1", TRIM(label1)) IF( LEN_TRIM(label2).GT.0 ) CALL attach(fidres, "/", "label2", TRIM(label2)) IF( LEN_TRIM(label3).GT.0 ) CALL attach(fidres, "/", "label3", TRIM(label3)) IF( LEN_TRIM(label4).GT.0 ) CALL attach(fidres, "/", "label4", TRIM(label4)) ! ! Job number jobnum = 0 ! ! Data group CALL creatg(fidres, "/data", "data") CALL creatg(fidres, "/data/var0d", "0d history arrays") CALL creatg(fidres, "/data/var1d","1d history arrays") CALL creatg(fidres, "/data/part", "part phase space") CALL creatg(fidres, "/data/fields", "Electro static potential and Er Ez fields") ! ! File group CALL creatg(fidres, "/files", "files") CALL attach(fidres, "/files", "jobnum", jobnum) CALL putarr(fidres, "/data/var1d/rgrid", rgrid) CALL putarr(fidres, "/data/var1d/zgrid", zgrid) ! Part and fields vectors ! Initialize time-dependent particle variables CALL creatd(fidres, 1, SHAPE(parts%R), "/data/part/R", "R") CALL creatd(fidres, 1, SHAPE(parts%Z), "/data/part/Z", "Z") CALL creatd(fidres, 1, SHAPE(parts%Z), "/data/part/THET", "THET") CALL creatd(fidres, 1, SHAPE(parts%UZ), "/data/part/UR", "UZ") CALL creatd(fidres, 1, SHAPE(parts%UR), "/data/part/UZ", "UR") CALL creatd(fidres, 1, SHAPE(parts%UTHET), "/data/part/UTHET", "UTHET") CALL creatd(fidres, 1, SHAPE(parts%Rindex), "/data/part/Rindex", "Rindex") CALL creatd(fidres, 1, SHAPE(parts%Zindex), "/data/part/Zindex", "Zindex") CALL creatd(fidres, 1, SHAPE(parts%partindex), "/data/part/partindex", "partindex") CALL creatd(fidres, 1, SHAPE(parts%pot), "/data/part/pot", "pot") CALL creatd(fidres, 0, SHAPE(time), "/data/part/time", "time" ) CALL creatd(fidres, 0, SHAPE(Ntotal), "/data/part/Nparts", "number of remaining parts in the simulation space") CALL creatd(fidres, 1, SHAPE(pot), "/data/fields/pot", "pot") CALL creatd(fidres, 1, SHAPE(Er), "/data/fields/Er", "Er") CALL creatd(fidres, 1, SHAPE(Ez), "/data/fields/Ez", "Ez") CALL creatd(fidres, 1, SHAPE(phi_spline), "/data/fields/phi", "spline form of Phi") CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments") CALL creatd(fidres, 0, SHAPE(time), "/data/fields/time", "time" ) - IF(mpisize .gt. 1) CALL creatd(fidres, 1, SHAPE(Nplocs_all), "/data/part/Nplocs_all", "Nplocs_all") + CALL creatd(fidres, 1, (/ 64 /), "/data/part/Nplocs_all", "Nplocs_all") CALL putarr(fidres, "/data/fields/Br", Br) CALL putarr(fidres, "/data/fields/Bz", Bz) CALL putarr(fidres, "/data/fields/Athet", Athet) CALL putarr(fidres, "/data/fields/volume", Volume) ! ! 1.2 Restart run ! ELSE CALL cp2bk(resfile) ! backup previous result file CALL openf(resfile, fidres,'rw', 'd') WRITE(*,'(3x,a,a)') TRIM(resfile), ' open' CALL getatt(fidres, "/files", "jobnum", jobnum) jobnum = jobnum+1 WRITE(*,'(3x,a,i3)') "Current Job Number =", jobnum CALL attach(fidres, "/files", "jobnum", jobnum) END IF ! ! Add input namelist variables as attributes of /data/input WRITE(str,'(a,i2.2)') "/data/input.",jobnum CALL creatg(fidres, TRIM(str)) CALL attach(fidres, TRIM(str), "job_time", job_time) CALL attach(fidres, TRIM(str), "extra_time", extra_time) CALL attach(fidres, TRIM(str), "dt", dt) CALL attach(fidres, TRIM(str), "tmax", tmax) CALL attach(fidres, TRIM(str), "nrun", nrun) CALL attach(fidres, TRIM(str), "nlres", nlres) CALL attach(fidres, TRIM(str), "nlsave", nlsave) CALL attach(fidres, TRIM(str), "newres", newres) CALL attach(fidres, TRIM(str), "nz", nz) CALL putarr(fidres, TRIM(str)//"/lz", lz) CALL attach(fidres, TRIM(str), "nplasma", nplasma) CALL attach(fidres, TRIM(str), "potinn", potinn) CALL attach(fidres, TRIM(str), "potout", potout) CALL attach(fidres, TRIM(str), "B0", B0) CALL attach(fidres, TRIM(str), "Rcurv", Rcurv) CALL attach(fidres, TRIM(str), "width", width) CALL attach(fidres, TRIM(str), "n0", n0) CALL attach(fidres, TRIM(str), "temp", temp) CALL attach(fidres, TRIM(str), "it0d", it0d) CALL attach(fidres, TRIM(str), "it2d", it2d) CALL attach(fidres, TRIM(str), "itparts", itparts) CALL attach(fidres, TRIM(str), "nlclassical", nlclassical) CALL attach(fidres, TRIM(str), "nlPhis", nlPhis) CALL attach(fidres, TRIM(str), "qsim", qsim) CALL attach(fidres, TRIM(str), "msim", msim) CALL attach(fidres, TRIM(str), "V", V) CALL attach(fidres, TRIM(str), "startstep", cstep) CALL attach(fidres, TRIM(str), "H0", H0) CALL attach(fidres, TRIM(str), "P0", P0) CALL putarr(fidres, TRIM(str)//"/femorder", femorder) CALL putarr(fidres, TRIM(str)//"/ngauss", ngauss) CALL putarr(fidres, TRIM(str)//"/nnr", nnr) CALL putarr(fidres, TRIM(str)//"/radii", radii) CALL putarr(fidres, TRIM(str)//"/plasmadim", plasmadim) ! Save STDIN of this run WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum INQUIRE(unit=lu_in, name=fname) CALL putfile(fidres, TRIM(str), TRIM(fname)) CLOSE(lu_in) ! ! Initialize buffers for 0d history arrays CALL htable_init(hbuf0, BUFSIZE) CALL set_htable_fileid(hbuf0, fidres, "/data/var0d") ! ! ! Initialize Xgrafix ! IF(nlxg) THEN CALL initw END IF END IF !________________________________________________________________________________ ! 2. Periodic diagnostics IF( kstep .NE. -1) THEN ! IF (mpisize .gt. 1) THEN Ntotal=sum(Nplocs_all) ELSE Ntotal=parts%Nptot END IF ! IF(modulo(step,ittext).eq. 0 .or. nlend) THEN WRITE(*,'(a,1x,i8.8,a1,i8.8,20x,a,1pe10.3)') '*** Timestep (this run/total) =', & & step, '/', cstep, 'Time =', time WRITE(*,'(a,4(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Eerr, Nbparts/Ntotal', epot, ekin, etot, etot-etot0, Ntotal,'/',nplasma END IF !________________________________________________________________________________ ! ! 2.1 0d history arrays ! IF(modulo(step,it0d).eq. 0 .or. nlend) THEN CALL add_record(hbuf0, "time", "simulation time", time) CALL add_record(hbuf0, "epot", "potential energy", epot) CALL add_record(hbuf0, "ekin", "kinetic energy", ekin) CALL add_record(hbuf0, "etot", "total energy", etot) CALL add_record(hbuf0, "etot0", "theoretical total energy", etot0) CALL add_record(hbuf0, "nbparts", "number of remaining parts in the simulation space", REAL(Ntotal,kind=db)) CALL htable_endstep(hbuf0) END IF ! ! 2.2 1d profiles IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL append(fidres, "/data/fields/time", time) CALL append(fidres, "/data/fields/pot", pot) CALL append(fidres, "/data/fields/Er", Er) CALL append(fidres, "/data/fields/Ez", Ez) CALL append(fidres, "/data/fields/phi", phi_spline) CALL append(fidres, "/data/fields/moments", moments) END IF ! ! 2.3 2d profiles IF(modulo(step,itparts).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' CALL append(fidres, "/data/part/time", time) CALL append(fidres, "/data/part/Nparts", REAL(Ntotal,kind=db)) IF ( isdataset(fidres,'/data/part/R') ) THEN CALL getdims(fidres, '/data/part/R', partsrank, partsdim) partsdim(1)=min(size(parts%R,1), partsdim(1)) CALL append(fidres, "/data/part/R", parts%R(1:partsdim(1))) CALL append(fidres, "/data/part/Z", parts%Z(1:partsdim(1))) CALL append(fidres, "/data/part/THET", parts%THET(1:partsdim(1))) CALL append(fidres, "/data/part/UZ", parts%UZ(1:partsdim(1))) CALL append(fidres, "/data/part/UR", parts%UR(1:partsdim(1))) CALL append(fidres, "/data/part/UTHET", parts%UTHET(1:partsdim(1))) CALL append(fidres, "/data/part/pot", parts%pot(1:partsdim(1))) CALL append(fidres, "/data/part/Rindex", REAL(parts%Rindex(1:partsdim(1)),kind=db)) CALL append(fidres, "/data/part/Zindex", REAL(parts%Zindex(1:partsdim(1)),kind=db)) CALL append(fidres, "/data/part/partindex", REAL(parts%partindex(1:partsdim(1)),kind=db)) END IF - IF (mpisize .gt. 1) CALL append(fidres, "/data/part/Nplocs_all", REAL(Nplocs_all,kind=db)) + Nplocs_all_save=0 + Nplocs_all_save(1:mpisize)=Nplocs_all(0:mpisize-1) + CALL append(fidres, "/data/part/Nplocs_all", REAL(Nplocs_all_save,kind=db)) ! END IF ! ! 2.4 3d profiles ! ! ! 2.5 Xgrafix ! IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0) THEN call xgevent CALL updt_xg_var CALL xgupdate END IF !________________________________________________________________________________ ! 3. Final diagnostics ELSE ! ! Flush 0d history array buffers CALL htable_hdf5_flush(hbuf0) ! ! Close all diagnostic files CALL closef(fidres) !________________________________________________________________________________ END IF ! END SUBROUTINE diagnose diff --git a/src/extra.c b/src/extra.c index 76fe577..5bdf91c 100644 --- a/src/extra.c +++ b/src/extra.c @@ -1,36 +1,36 @@ -#include - -/* double seconds_(void) -{ - // Return elapsed wall time in s. - - struct timeval tp; - struct timezone tzp; - int i; - i = gettimeofday(&tp,&tzp); - return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 ); -} -/**********************************************/ -#include - -void quit_(); -void dump_(char *filename, int *l); - -void Dump(filename) -char *filename; -{ - /* The user's dump routine should go here. */ - int l = strlen(filename); - // dump_(filename, &l); - -} /* End DUMP */ - -/**********************************************/ - -void Quit() -{ - /* The user's quit routine should go here. */ - - // quit_(); - -} /* End QUIT */ +#include + +/* double seconds_(void) +{ + // Return elapsed wall time in s. + + struct timeval tp; + struct timezone tzp; + int i; + i = gettimeofday(&tp,&tzp); + return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 ); +} +/**********************************************/ +#include + +void quit_(); +void dump_(char *filename, int *l); + +void Dump(filename) +char *filename; +{ + /* The user's dump routine should go here. */ + int l = strlen(filename); + // dump_(filename, &l); + +} /* End DUMP */ + +/**********************************************/ + +void Quit() +{ + /* The user's quit routine should go here. */ + + // quit_(); + +} /* End QUIT */ diff --git a/src/fields_mod.f90 b/src/fields_mod.f90 index 402bd12..c10f004 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,706 +1,710 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for initializing the magnetic field and solve the Poisson equation !------------------------------------------------------------------------------ MODULE fields USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, splrz, nlperiod, phinorm, nlPhis, nrank, mpirank, mpisize, step, it2d USE beam, ONLY: parts USE bsplines USE mumps_bsplines USE mpi Use omp_lib +Use mpihelper, ONLY: db_type IMPLICIT NONE REAL(kind=db), allocatable, SAVE :: matcoef(:,:), phi_spline(:), vec1(:), vec2(:) REAL(kind=db), allocatable, SAVE :: loc_moments(:,:) INTEGER, SAVE:: loc_zspan TYPE(mumps_mat), SAVE :: femat ! Finite Element Method matrix CONTAINS !--------------------------------------------------------------------------- !> @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, moments, leftproc, rightproc, Zbounds, volume USE bsplines USE mumps_bsplines USE mpihelper INTEGER :: nrz(2) ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod) ! Calculate dimension of splines nrz(1)=nz nrz(2)=nr CALL get_dim(splrz, nrank, nrz, femorder) ALLOCATE(matcoef(nrank(1),nrank(2))) ALLOCATE(pot((nr+1)*(nz+1))) ALLOCATE(rhs(nrank(1)*nrank(2))) IF (mpirank .eq. 0) THEN ALLOCATE(moments(10,nrank(1)*nrank(2))) ELSE ALLOCATE(moments(0,0)) END IF ALLOCATE(phi_spline(nrank(1)*nrank(2))) ALLOCATE(volume(nrank(1)*nrank(2))) ALLOCATE(Br((nr+1)*(nz+1)),Bz((nr+1)*(nz+1))) ALLOCATE(Athet((nr+1)*(nz+1))) ALLOCATE(Er((nr+1)*(nz+1)),Ez((nr+1)*(nz+1))) ! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed) CALL magnet ! Calculate and factorise FEM matrix (depends only on mesh) CALL fematrix CALL factor(femat) CALL comp_volume loc_zspan=Zbounds(mpirank+1)-Zbounds(mpirank)+femorder(1) ALLOCATE(loc_moments(10,loc_zspan*nrank(2))) IF(mpisize .gt. 1) THEN CALL init_overlaps(nrank, femorder, Zbounds(mpirank), Zbounds(mpirank+1), loc_moments, mpirank, leftproc, rightproc) END IF END SUBROUTINE fields_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 ! !--------------------------------------------------------------------------- SUBROUTINE rhscon(p) USE bsplines USE mumps_bsplines USE mpi USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma, moments, nlend, step, it2d USE beam, ONLY: particles USE mpihelper type(particles), INTENT(INOUT):: p REAL(kind=db) :: chargecoeff loc_moments=0 ! Reset the moments matrix rhs=0 ! reset rhs chargecoeff=0.5/pi*omegap/omegac*V/real(nplasma) ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (p%Nptot .ne. 0) THEN IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL deposit_moments(p) ELSE CALL deposit_charge(p) END IF END IF ! If we are using MPI parallelism, reduce the rhs on the root process IF(mpisize .gt. 1) THEN CALL fields_comm ELSE moments=loc_moments END IF IF(mpirank.eq.0 )THEN IF(nlphis) THEN ! We consider self-consistent interactions rhs=moments(1,1:nrank(1)*nrank(2))*chargecoeff ! Normalized charge density simulated by each macro particle ELSE rhs=0 END IF END IF ! Apply the Dirichlet boundary conditions on the coaxial insert (if present) IF(rgrid(0).ne.0.0_db) rhs(1:nrank(1))=potinn ! Apply the Dirichlet boundary conditions on the cylinder walls rhs((nrank(2)-1)*nrank(1)+1:nrank(2)*nrank(1))=potout END SUBROUTINE rhscon !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles moments (q,v,v^2) from p on the grid ! !--------------------------------------------------------------------------- SUBROUTINE deposit_moments(p) USE bsplines USE mumps_bsplines USE mpi USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma, moments, nlend, step, it2d, Zbounds USE beam, ONLY: particles USE mpihelper TYPE(particles), INTENT(IN):: p INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER,DIMENSION(:),ALLOCATABLE::zleft, rleft REAL(kind=db) :: vr, vthet, vz, coeff, chargecoeff REAL(kind=db), ALLOCATABLE :: fun(:,:,:), fun2(:,:,:) INTEGER:: num_threads num_threads=omp_get_max_threads() nbunch=p%Nptot/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 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 chargecoeff=0.5/pi*omegap/omegac*V/real(nplasma) ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (p%Nptot .ne. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(i, zleft, rleft, jw, it, iend, irow, jcol, mu, k, vr, vz, vthet, coeff) PRIVATE(fun, fun2) DO i=1,p%Nptot,nbunch ! Avoid segmentation fault by accessing non relevant data iend=min(i+nbunch-1,p%Nptot) 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)) ! Compute the value of the splines at the particles positions CALL basfun(p%Z(i:iend), splrz%sp1, fun, zleft(1:k)+1) CALL basfun(p%R(i:iend), splrz%sp2, fun2, rleft(1:k)+1) !$OMP SIMD DO k=1,(iend-i+1) DO jw=1,(femorder(2)+1) DO it=1,(femorder(1)+1) irow=zleft(k)+it-Zbounds(mpirank) jcol=rleft(k)+jw mu=irow+(jcol-1)*(loc_zspan) coeff=fun(it,0,k)*fun2(jw,0,k) ! Add contribution of particle nbunch to rhs grid point mu vr=parts%UR(i+k-1)/parts%Gamma(i+k-1) vz=parts%UZ(i+k-1)/parts%Gamma(i+k-1) vthet=parts%UTHET(i+k-1)/parts%Gamma(i+k-1) !$OMP ATOMIC UPDATE loc_moments(1, mu) = loc_moments(1, mu) + coeff !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(2, mu) = loc_moments(2, mu) + coeff*vr !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(3, mu) = loc_moments(3, mu) + coeff*vthet !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(4, mu) = loc_moments(4, mu) + coeff*vz !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(5, mu) = loc_moments(5, mu) + coeff*vr*vr !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(6, mu) = loc_moments(6, mu) + coeff*vr*vthet !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(7, mu) = loc_moments(7, mu) + coeff*vr*vz !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(8, mu) = loc_moments(8, mu) + coeff*vthet*vthet !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(9, mu) = loc_moments(9, mu) + coeff*vthet*vz !$OMP END ATOMIC !$OMP ATOMIC UPDATE loc_moments(10,mu) = loc_moments(10,mu) + coeff*vz*vz !$OMP END ATOMIC END DO END DO END DO !$OMP END SIMD END DO !$OMP END PARALLEL DO END IF DEALLOCATE(fun,fun2, zleft, rleft) 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 ! !--------------------------------------------------------------------------- SUBROUTINE deposit_charge(p) USE bsplines USE mumps_bsplines USE mpi USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma, moments, nlend, step, it2d, Zbounds USE beam, ONLY: particles USE mpihelper TYPE(particles), INTENT(IN):: p INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER,DIMENSION(:),ALLOCATABLE::zleft, rleft REAL(kind=db) :: vr, vthet, vz, coeff, chargecoeff REAL(kind=db), ALLOCATABLE :: fun(:,:,:), fun2(:,:,:) INTEGER:: num_threads num_threads=omp_get_max_threads() nbunch=p%Nptot/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 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 chargecoeff=0.5/pi*omegap/omegac*V/real(nplasma) ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (p%Nptot .ne. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(i, zleft, rleft, jw, it, iend, irow, jcol, mu, k, vr, vz, vthet, coeff) PRIVATE(fun, fun2) DO i=1,p%Nptot,nbunch ! Avoid segmentation fault by accessing non relevant data iend=min(i+nbunch-1,p%Nptot) 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)) ! Compute the value of the splines at the particles positions CALL basfun(p%Z(i:iend), splrz%sp1, fun, zleft(1:k)+1) CALL basfun(p%R(i:iend), splrz%sp2, fun2, rleft(1:k)+1) !$OMP SIMD DO k=1,(iend-i+1) DO jw=1,(femorder(2)+1) DO it=1,(femorder(1)+1) irow=zleft(k)+it-Zbounds(mpirank) jcol=rleft(k)+jw mu=irow+(jcol-1)*(loc_zspan) ! Add contribution of particle nbunch to rhs grid point mu !$OMP ATOMIC update loc_moments(1, mu) = loc_moments(1, mu) + fun(it,0,k)*fun2(jw,0,k) !$OMP END ATOMIC END DO END DO END DO !$OMP END SIMD END DO !$OMP END PARALLEL DO END IF DEALLOCATE(fun,fun2, zleft, rleft) END subroutine deposit_charge SUBROUTINE fields_comm USE mpihelper USE Basic, ONLY: moments, Zbounds, mpirank, nlend, it2d, step, leftproc, rightproc INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) INTEGER:: overlap_type INTEGER:: rcvoverlap_type displs=Zbounds(0:mpisize-1) counts=Zbounds(1:mpisize)-Zbounds(0:mpisize-1) counts(mpisize)=counts(mpisize)+femorder(1) IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL start_momentscomm(mpirank, leftproc, rightproc, loc_moments, nrank, femorder, loc_zspan-femorder(1)) IF(leftproc .lt. mpirank) THEN - !$OMP PARALLEL DO SIMD DEFAULT(SHARED), PRIVATE(i) + !$OMP PARALLEL DO SIMD DEFAULT(SHARED) DO j=1,femorder(1) DO i=1,nrank(2) loc_moments(1:10,(i-1)*loc_zspan+j) = loc_moments(1:10,(i-1)*loc_zspan+j)& & + momentsoverlap_buffer(10*(nrank(2)*(j-1)+i-1)+1:10*(nrank(2)*(j-1)+i)) END DO END DO !$OMP END PARALLEL DO SIMD END IF ! Set communication vector type overlap_type=momentsoverlap_type rcvoverlap_type=rcvmomentsoverlap_type ELSE CALL start_rhscomm(mpirank, leftproc, rightproc, loc_moments, nrank, femorder, loc_zspan-femorder(1)) IF(leftproc .lt. mpirank) THEN - !$OMP PARALLEL DO SIMD DEFAULT(SHARED), PRIVATE(i) + !$OMP PARALLEL DO SIMD DEFAULT(SHARED) DO j=1,femorder(1) DO i=1,nrank(2) loc_moments(1, (i-1)*loc_zspan+j) = loc_moments(1,(i-1)*loc_zspan+j)& & + rhsoverlap_buffer(nrank(2)*(j-1)+i) END DO END DO !$OMP END PARALLEL DO SIMD END IF ! Set communication vector type overlap_type=rhsoverlap_type rcvoverlap_type=rcvrhsoverlap_type END IF - CALL MPI_GATHERV(loc_moments(1,1), counts(mpirank+1), overlap_type, & + IF(mpirank.eq.0) THEN + moments=0 + END IF + CALL MPI_GATHERV(loc_moments, counts(mpirank+1), overlap_type, & & moments, counts, displs, rcvoverlap_type, 0, MPI_COMM_WORLD, ierr) END SUBROUTINE fields_comm !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Solves Poisson equation using FEM. Distributes the result on all MPI workers and interpolate the electric forces !> for each particle. ! !--------------------------------------------------------------------------- SUBROUTINE poisson USE basic, ONLY: rhs, nrank, pot USE bsplines USE mumps_bsplines USE futils INTEGER:: ierr, jder(2) jder(1)=0 jder(2)=0 ! Only the root process solves Poisson IF(mpirank.eq.0) CALL bsolve(femat,rhs,phi_spline) ! Broadcast the solution to all MPI workers - CALL MPI_Bcast(phi_spline(1), nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Bcast(phi_spline(1), nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr) matcoef= reshape(phi_spline, (/nrank(1),nrank(2)/)) ! Evaluate the electric potential on the grid points CALL gridval(splrz, vec1, vec2, pot, jder, matcoef) IF( mpirank .eq. 0 .and. modulo(step,it2d) .eq. 0) THEN ! On the root process, compute the electric field for diagnostic purposes CALL gridval(splrz, vec1, vec2, Ez, (/1,0/)) CALL gridval(splrz, vec1, vec2, Er, (/0,1/)) Ez=-Ez Er=-Er END IF END SUBROUTINE poisson SUBROUTINE EForcescomp(p) Use beam, ONLY: particles TYPE(particles), INTENT(INOUT):: p INTEGER:: i, iend INTEGER:: nbunch=64 INTEGER:: num_threads num_threads=omp_get_max_threads() nbunch=p%Nptot/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 ! Evaluate the electric potential and field at the particles position !$OMP PARALLEL DO SIMD DEFAULT(SHARED), PRIVATE(iend) DO i=1,p%Nptot,nbunch ! Avoid segmentation fault by accessing non relevant data iend=min(i+nbunch-1,p%Nptot) CALL getgrad(splrz, p%Z(i:iend), p%R(i:iend), p%pot(i:iend), p%EZ(i:iend), p%ER(i:iend)) p%EZ(i:iend)=-p%Ez(i:iend) p%ER(i:iend)=-p%ER(i:iend) END DO !$OMP END PARALLEL DO SIMD END SUBROUTINE EForcescomp !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrix USE bsplines USE mumps_bsplines REAL(kind=db), ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:) REAL(kind=db), ALLOCATABLE :: f(:,:), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:,:), fun2(:,:) REAL(kind=db) :: contrib INTEGER, ALLOCATABLE :: idert(:), iderw(:) INTEGER :: i, j, k, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms kterms=2 !Number of terms in weak form 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 ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms)) !Pointers on the order of derivatives ! 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 ! Initialisation of FEM matrix CALL init(nrank(1)*nrank(2),2,femat) ! Assemble FEM matrix 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 get_gauss(splrz%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(splrz%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO DO igauss=1,ngauss(1)*ngauss(2) ! 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 coefeq(xgauss(igauss,:), idert, iderw, coefs) DO iterm=1,kterms ! Loop on the two integration dimensions DO jt=1,(1+femorder(1))*(femorder(2)+1) DO iw=1,(1+femorder(1))*(femorder(2)+1) irow=i+f(jt,1)-1; jcol=j+f(jt,2)-1 irow2=i+f(iw,1)-1; jcol2=j+f(iw,2)-1 mu=irow+(jcol-1)*nrank(1) mu2=irow2+(jcol2-1)*nrank(1) contrib=fun(f(jt,1),idert(iterm))* fun(f(iw,1),idert(iterm))* & & fun2(f(jt,2),iderw(iterm))*fun2(f(iw,2),iderw(iterm))* & & wgauss(igauss)*xgauss(igauss,2) CALL updtmat(femat, mu, mu2, contrib) END DO END DO END DO END DO END DO END DO ! Impose Dirichlet boundary conditions on FEM matrix CALL fe_dirichlet DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE(f, aux) DEALLOCATE(idert, iderw, coefs, fun,fun2) END SUBROUTINE fematrix SUBROUTINE comp_volume USE bsplines USE mumps_bsplines USE basic, ONLY: Volume REAL(kind=db), ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:) REAL(kind=db), ALLOCATABLE :: f(:,:), aux(:) REAL(kind=db), ALLOCATABLE :: fun(:,:), fun2(:,:) INTEGER :: i, j, k, jt, irow, jcol, mu, igauss ALLOCATE(fun(1:femorder(1)+1,0:0),fun2(1:femorder(2)+1,0:0))!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 ! 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 ! Assemble Volume matrix 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 get_gauss(splrz%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(splrz%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO DO igauss=1,ngauss(1)*ngauss(2) ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss,1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss,2), splrz%sp2, fun2, j) 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) volume(mu)=volume(mu)+2*pi*fun(f(jt,1),0) * fun2(f(jt,2),0)* wgauss(igauss)*xgauss(igauss,2) END DO END DO END DO END DO DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE(f, aux) DEALLOCATE(fun,fun2) END SUBROUTINE comp_volume !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Imposes the dirichlet boundary conditions on the FEM matrix. !--------------------------------------------------------------------------- SUBROUTINE fe_dirichlet REAL(kind=db), ALLOCATABLE :: arr(:) INTEGER :: i ALLOCATE(arr(nrank(1)*nrank(2))) DO i=1,nrank(1)*nrank(2) CALL getrow(femat,i,arr) END DO 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, c) REAL(kind=db), INTENT(in) :: x(2) INTEGER, INTENT(out) :: idt(:), idw(:) REAL(kind=db), INTENT(out) :: c(SIZE(idt)) c(1) = x(2) idt(1) = 0 idw(1) = 1 c(2) = x(2) idt(2) = 1 idw(2) = 0 END SUBROUTINE coefeq !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the magnetic field on the grid according to a magnetic mirror. !--------------------------------------------------------------------------- SUBROUTINE magnet USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm USE constants, ONLY: Pi REAL(kind=db) :: rg, zg,halfLz, MirrorRatio INTEGER :: i, rindex halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) DO i=1,(nr+1)*(nz+1) rindex=(i-1)/(nz+1) rg=rgrid(rindex) zg=zgrid(i-rindex*(nz+1)-1)-halfLz Br(i)=-B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)/bnorm Bz(i)=B0*(1-MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm Athet(i)=0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm)) END DO END SUBROUTINE magnet !________________________________________________________________________________ !Modified Bessel functions of the first kind of the zero order FUNCTION bessi0(x) REAL(kind=db) :: bessi0,x REAL(kind=db) :: ax REAL(kind=db) p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2 ,0.916281d-2, & & -0.2057706d-1, 0.2635537d-1,-0.1647633d-1, 0.392377d-2 / if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) endif return END FUNCTION bessi0 !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) REAL(kind=db) :: bessi1,x REAL(kind=db) :: ax REAL(kind=db) p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 !--------------------------------------------------------------------------- !> @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, moments DEALLOCATE(matcoef) DEALLOCATE(pot) DEALLOCATE(rhs) DEALLOCATE(phi_spline) DEALLOCATE(moments) DEALLOCATE(Br,Bz) DEALLOCATE(Er,Ez) DEALLOCATE(vec1,vec2) Call DESTROY_SP(splrz) END SUBROUTINE clean_fields END MODULE fields diff --git a/src/gcc.mk b/src/gcc.mk new file mode 100644 index 0000000..9b87a2d --- /dev/null +++ b/src/gcc.mk @@ -0,0 +1,13 @@ +CC = gcc +CFLAGS = +FC = mpif90 +#FFLAGS = -g -traceback -check bounds -mkl=cluster +RELEASEFLAGS = -fopenmp -O2 -march=native +DEBUGFLAGS = -g -fopenmp -fbacktrace \ + -fcheck=all -Wall +PROFILEFLAGS= -g + +F90 = $(FC) +F90FLAGS= $(RELEASEFLAGS) +LDFLAGS = +LIBS = diff --git a/src/intel.mk b/src/intel.mk index b85ebb7..bb1ba26 100644 --- a/src/intel.mk +++ b/src/intel.mk @@ -1,14 +1,15 @@ CC = icc CFLAGS = FC = mpiifort #FFLAGS = -g -traceback -check bounds -mkl=cluster RELEASEFLAGS = -mkl=cluster -qopenmp -O3 -xHost -simd -DEBUGFLAGS = -g -O2 -xHost -qopenmp -traceback -mkl=cluster -simd\ +DEBUGFLAGS = -g -O0 -xHost -qopenmp -traceback -mkl=cluster -simd\ -check all -check bounds -check noarg_temp_created \ - -warn all -ftrapuv -fpe3 -debug extended + -warn all -ftrapuv -fpe0 -debug extended \ + -check uninit -debug all PROFILEFLAGS= -g F90 = $(FC) F90FLAGS= $(RELEASEFLAGS) LDFLAGS = LIBS = diff --git a/src/main.f90 b/src/main.f90 index dc49431..cc7a906 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,85 +1,90 @@ PROGRAM main ! ! Skeleton for a time dependent program ! Note: Even in this sequential version, MPI is required ! because of FUTILS (more specifcally because ! of the HASTABLE module)! ! USE basic USE mpi USE bsplines USE mumps_bsplines USE futils IMPLICIT NONE INTEGER:: required, provided ! ! required=MPI_THREAD_FUNNELED CALL mpi_init_thread(required,provided,ierr) + IF(provided .lt. required) CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) + IF(ierr .ne. 0) CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) + IF(ierr .ne. 0) CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, mpisize, ierr) + IF(ierr .ne. 0) CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) + !-------------------------------------------------------------------------------- ! 1. Prologue CALL timera(0, 'Prologue') CALL daytim('Start at') ! ! Define data specific to run ! CALL basic_data !Definition of global variables and input paramaters loading step=0 ! IF( .NOT. nlres ) THEN CALL newrun !not implemented yet ELSE CALL restart !not implemented yet END IF ! ! Compute auxilliary values ! CALL auxval !time independent values ! ! Initial conditions ! IF( .NOT. nlres ) THEN CALL inital !plasma initialisation ELSE CALL resume !loads restart.h5 file END IF ! ! Start or restart the run ! CALL start !not implemented yet ! ! Initial diagnostocs ! CALL diagnose(0) CALL timera(1, 'Prologue') !-------------------------------------------------------------------------------- ! 2. Time stepping CALL timera(0, 'Main loop') ! DO step = step+1 cstep = cstep+1 time = time+dt CALL tesend CALL stepon CALL diagnose(step) IF(modulo(step,itrestart) .eq. 0) CALL chkrst(1) IF( nlend ) EXIT END DO CALL timera(1, 'Main loop') !-------------------------------------------------------------------------------- ! 9. Epilogue CALL timera(0, 'Epilogue') ! CALL diagnose(-1) CALL endrun IF(mpirank .eq. 0) THEN CALL timera(1, 'Epilogue') CALL timera(9, '') CALL timera(-1, '') CALL daytim('Done at ') END IF CALL mpi_finalize(ierr) END PROGRAM main diff --git a/src/maxwellsource_mod.f90 b/src/maxwellsource_mod.f90 index cc34b2e..cf7914c 100644 --- a/src/maxwellsource_mod.f90 +++ b/src/maxwellsource_mod.f90 @@ -1,148 +1,148 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: maxwellsource ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Adds particle in the simulation according to a maxwellian distribution in velocity !> and a uniform distribution in space !------------------------------------------------------------------------------ MODULE maxwellsource ! USE constants USE mpi USE mpihelper USE basic, ONLY: mpirank, mpisize, vnorm, rnorm, Zbounds, zgrid, & & nlclassical, nlmaxwellsource USE beam USE distrib IMPLICIT NONE PRIVATE REAL(kind=db), SAVE :: frequency !< Number of macro particles added per second over the spawning region REAL(kind=db), SAVE :: temperature !< temperature used for the Maxwellian velocity distribution REAL(kind=db), SAVE :: rlimits(2) !< radial limits in which particles will be spawned REAL(kind=db), SAVE :: zlimits(2) !< axial limits in which particles will be spawned REAL(kind=db), SAVE :: last_t=0 !< last time when a macro-particle was added to the simulation REAL(kind=db), SAVE :: vth !< Normalized thermal velocity computed from temperature REAL(kind=db), SAVE :: loc_zlimits(2)!< Local axial limits in which particles will be spawned used for current mpi process REAL(kind=db), SAVE :: loc_rlimits(2)!< Local radial limits in which particles will be spawned used for current mpi process REAL(kind=db), SAVE :: time_start=-1.0 !< time at which the source is turned on REAL(kind=db), SAVE :: time_end=-1.0 !< time at which the source is turned off NAMELIST /maxwellsourceparams/ frequency, temperature, rlimits, zlimits, time_start, time_end PUBLIC:: maxwellsource_init, maxwellsource_inject contains subroutine maxwellsource_init(lu_in, time) implicit none INTEGER, INTENT(IN)::lu_in REAL(kind=db), INTENT(IN):: time REAL(kind=db):: surface REAL(kind=db):: remsurface READ(lu_in, maxwellsourceparams) IF(mpirank .eq. 0) THEN WRITE(*, maxwellsourceparams) ELSE END IF vth=sqrt(3*kb*temperature/me)/vnorm !thermal velocity last_t=time surface=(rlimits(2)-rlimits(1))*(zlimits(2)-zlimits(1)) loc_rlimits=rlimits/rnorm IF (zlimits(1) .lt. zgrid(Zbounds(mpirank))*rnorm) THEN remsurface=(rlimits(2)-rlimits(1))*(zgrid(Zbounds(mpirank))*rnorm-zlimits(1)) frequency=frequency*(1-remsurface/surface) loc_zlimits(1)=zgrid(Zbounds(mpirank)) ELSE loc_zlimits(1)=zlimits(1)/rnorm END IF IF (zlimits(2) .gt. zgrid(Zbounds(mpirank+1))*rnorm) THEN remsurface=(rlimits(2)-rlimits(1))*(zlimits(2)-zgrid(Zbounds(mpirank+1))*rnorm) frequency=frequency*(1-remsurface/surface) loc_zlimits(2)=zgrid(Zbounds(mpirank+1)) ELSE loc_zlimits(2)=zlimits(2)/rnorm END IF WRITE(*,*) "init frequency : ", frequency End subroutine maxwellsource_init subroutine maxwellsource_inject(time) implicit none REAL(kind=db), INTENT(IN) :: time Type(particle), ALLOCATABLE, Dimension(:):: newparts INTEGER:: npartsadd INTEGER:: i REAL(kind=db), ALLOCATABLE:: VR(:),VZ(:),VTHET(:),R(:),Z(:) ! check if source is on IF(.not. maxwellsource_on(time)) THEN RETURN END IF ! Number of particles to add at this time step npartsadd=floor((time-last_t)*frequency) IF (npartsadd.gt. 0) THEN ALLOCATE(newparts(npartsadd)) ALLOCATE(VR(npartsadd),VZ(npartsadd),VTHET(npartsadd),R(npartsadd),Z(npartsadd)) ! Initial velocities distribution CALL lodgaus(0,VZ) CALL lodgaus(0,VR) CALL lodgaus(0,VTHET) CALL lodlinr(0,R,loc_rlimits(1),loc_rlimits(2)) CALL loduni (0,Z) DO i=1,npartsadd newparts(i)%UR=VR(i)*vth newparts(i)%UZ=VZ(i)*vth newparts(i)%UTHET=VTHET(i)*vth newparts(i)%R=R(i) newparts(i)%Z=Z(i)*(loc_zlimits(2)-loc_zlimits(1)) + loc_zlimits(1) newparts(i)%THET=0 IF (nlclassical) THEN newparts(i)%gamma=1.0 ELSE newparts(i)%gamma=sqrt(1/(1-newparts(i)%UZ**2-newparts(i)%UR**2-newparts(i)%UTHET**2)) newparts(i)%UR=VR(i)*newparts(i)%gamma newparts(i)%UZ=VZ(i)*newparts(i)%gamma newparts(i)%UTHET=VTHET(i)*newparts(i)%gamma END IF END DO CALL add_created_part(parts, newparts) last_t=time END IF end subroutine maxwellsource_inject logical function maxwellsource_on(time) - real(kind=db), intent(in):: time + REAL(kind=db), intent(in):: time maxwellsource_on=.false. IF (time_start .lt. 0 .and. time_end .lt. 0) THEN maxwellsource_on = .true. ELSE IF (time .gt. time_start .and. time .lt. time_end ) THEN maxwellsource_on = .true. END IF maxwellsource_on=nlmaxwellsource .and. maxwellsource_on end function End Module maxwellsource diff --git a/src/mpihelper_mod.f90 b/src/mpihelper_mod.f90 index 54cb44f..e0bb82e 100644 --- a/src/mpihelper_mod.f90 +++ b/src/mpihelper_mod.f90 @@ -1,379 +1,375 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: mpihelper ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for setting up the MPI variables used in the communications. !------------------------------------------------------------------------------ MODULE mpihelper USE constants USE mpi IMPLICIT NONE !> Structure containing the simulation parameters read from the input file TYPE BASICDATA LOGICAL :: nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical INTEGER :: nplasma, nz, it0d, it2d, itparts, & & itgraph, distribtype, nblock, nrun REAL(kind=db) :: job_time, extra_time, tmax, dt, & & potinn, potout, B0, n0, temp, & & Rcurv, width, H0, P0 INTEGER, DIMENSION(2) :: femorder, ngauss INTEGER, DIMENSION(3) :: nnr REAL(kind=db), DIMENSION(2):: lz REAL(kind=db), DIMENSION(4):: radii, plasmadim CHARACTER(len=64) :: resfile= "results.h5" LOGICAL :: partperiodic INTEGER :: samplefactor END TYPE BASICDATA !> Structure containing a single particle position and velocity used in MPI communications. TYPE particle INTEGER :: Rindex, Zindex, partindex REAL(kind=db) :: R, Z, THET,& & UZ, UR, UTHET, & & Gamma END TYPE particle INTEGER, SAVE :: basicdata_type !< Stores the MPI data type used for communicating basicdata INTEGER, SAVE :: particle_type !< Stores the MPI data type used for particles communication INTEGER, SAVE :: rhsoverlap_type !< Stores the MPI data type used for the communication of a rhs column + INTEGER, SAVE :: db_type !< Stores the MPI data type used for the communication of a REAL(kind=db) INTEGER, SAVE :: momentsoverlap_type !< Stores the MPI data type used for the communication of a column of a grid variable INTEGER, SAVE :: rcvrhsoverlap_type !< Stores the MPI data type used for the receive communication of a rhs column INTEGER, SAVE :: rcvmomentsoverlap_type !< Stores the MPI data type used for the receive communication of a column of a grid variable + INTEGER, SAVE:: db_sum_op !< Store the MPI sum operation for db_type REAL(kind=db), ALLOCATABLE, SAVE:: rhsoverlap_buffer(:) !< buffer used for storing the rhs ghost cells !< received from the left or right MPI process REAL(kind=db), ALLOCATABLE, SAVE:: momentsoverlap_buffer(:) !< buffer used for storing the moments ghost cells !< received from the left or right MPI process - INTEGER, SAVE:: momentsoverlap_requests(2) = MPI_REQUEST_NULL - INTEGER, SAVE:: rhsoverlap_requests(2) = MPI_REQUEST_NULL + !INTEGER, SAVE:: momentsoverlap_requests(2) = MPI_REQUEST_NULL + !INTEGER, SAVE:: rhsoverlap_requests(2) = MPI_REQUEST_NULL INTEGER:: rhsoverlap_tag= 20 INTEGER:: momentsoverlap_tag= 30 CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Initialize the MPI types used for inter process communications ! !--------------------------------------------------------------------------- SUBROUTINE init_mpitypes - + IMPLICIT NONE + INTEGER:: ierr + CALL MPI_TYPE_CREATE_F90_REAL(dprequestedprec,MPI_UNDEFINED,db_type,ierr) + CALL MPI_Type_commit(db_type,ierr) + CALL MPI_Op_Create(DB_sum, .true., db_sum_op, ierr) CALL init_particlempi CALL init_basicdatampi END SUBROUTINE init_mpitypes + SUBROUTINE DB_sum(INVEC, INOUTVEC, LEN, TYPE) + REAL(kind=db):: INVEC(LEN), INOUTVEC(LEN) + INTEGER:: LEN, TYPE + INTEGER:: i -!--------------------------------------------------------------------------- + Do i=1,LEN + INOUTVEC(i)=INVEC(i)+INOUTVEC(i) + END DO + + END SUBROUTINE DB_sum + +!-------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Initialize the MPI communicators used for allreduce between neighbors ! !> @param[in] mpirank rank of the MPI process in MPI_COMM_WORLD ! !--------------------------------------------------------------------------- SUBROUTINE init_overlaps(nrank, femorder, zlimleft, zlimright, moments, mpirank, left_proc, right_proc) REAL(kind=db), ASYNCHRONOUS:: moments(:,:) INTEGER, INTENT(IN):: nrank(:), femorder(:), left_proc, right_proc, zlimright, mpirank, zlimleft - INTEGER:: ierr - ALLOCATE(rhsoverlap_buffer(nrank(2)*(femorder(1)))) - ALLOCATE(momentsoverlap_buffer(10*nrank(2)*(femorder(1)))) + ALLOCATE(rhsoverlap_buffer(nrank(2)*femorder(1))) + ALLOCATE(momentsoverlap_buffer(10*nrank(2)*femorder(1))) WRITE(*,*)"nrank, femorder", nrank, femorder ! Initialize the MPI column overlap type for rhs - CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(1), 1, 10, MPI_DOUBLE_PRECISION, rhsoverlap_type) + CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(1), 1, 10, db_type, rhsoverlap_type) ! Initialize the MPI grid col type - CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(1), 10, 1, MPI_DOUBLE_PRECISION, momentsoverlap_type) + CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(1), 10, 1, db_type, momentsoverlap_type) ! Initialize the MPI receive column overlap type for rhs - CALL init_coltypempi(nrank(2), nrank(1), 1, 10, MPI_DOUBLE_PRECISION, rcvrhsoverlap_type) + CALL init_coltypempi(nrank(2), nrank(1), 1, 10, db_type, rcvrhsoverlap_type) ! Initialize the MPI receive grid col type - CALL init_coltypempi(nrank(2), nrank(1), 10, 1, MPI_DOUBLE_PRECISION, rcvmomentsoverlap_type) - - ! If the processor on the right has actually higher z positions - IF(right_proc .gt. mpirank) THEN - CALL MPI_SEND_INIT(moments(1,zlimright+1), femorder(1), rhsoverlap_type, right_proc, rhsoverlap_tag, & - & MPI_COMM_WORLD, rhsoverlap_requests(1), ierr ) - IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in rhs send_init" - CALL MPI_SEND_INIT(moments(1,zlimright+1), femorder(1), momentsoverlap_type, right_proc, momentsoverlap_tag, & - & MPI_COMM_WORLD, momentsoverlap_requests(1), ierr ) - IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in moments send_init" - END IF + CALL init_coltypempi(nrank(2), nrank(1), 10, 1, db_type, rcvmomentsoverlap_type) - ! If the processor on the left has actually lower z positions - IF(left_proc .lt. mpirank) THEN - CALL MPI_RECV_INIT(rhsoverlap_buffer, nrank(2)*(femorder(1)), MPI_DOUBLE_PRECISION, left_proc, rhsoverlap_tag, & - & MPI_COMM_WORLD, rhsoverlap_requests(2), ierr ) - IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in rhs recv_init" - CALL MPI_RECV_INIT(momentsoverlap_buffer, 10*nrank(2)*(femorder(1)), MPI_DOUBLE_PRECISION, left_proc, & - & momentsoverlap_tag, MPI_COMM_WORLD, momentsoverlap_requests(2), ierr ) - IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in moments recv_init" - END IF END SUBROUTINE init_overlaps -SUBROUTINE start_persistentcomm(requests, mpirank, leftproc, rightproc, moments, nrank, femorder, zlimright) +SUBROUTINE start_persistentcomm(requests, mpirank, leftproc, rightproc) INTEGER, ASYNCHRONOUS, INTENT(INOUT):: requests(:) INTEGER, INTENT(IN):: mpirank, leftproc, rightproc -REAL(kind=db), DIMENSION(:,:), INTENT(INOUT):: moments -INTEGER, INTENT(IN):: nrank(2), femorder(2), zlimright INTEGER:: ierr INTEGER:: stats(MPI_STATUS_SIZE,2) LOGICAL:: completed=.false. IF(leftproc .lt. mpirank) THEN ! Start to receive CALL MPI_START(requests(2),ierr) IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in recv_init" END IF IF(rightproc .gt. mpirank) THEN ! Start to send CALL MPI_START(requests(1),ierr) IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in send_init" END IF IF(leftproc .lt. mpirank) THEN ! Start to receive completed=.FALSE. DO WHILE(.not. completed) CALL MPI_TEST(requests(2), completed,stats(:,2),ierr) END DO WRITE(*,*)"status 2", completed, stats(:,2) !CALL MPI_WAIT(requests(2),stats(:,2),ierr) !WRITE(*,*)"status 2", stats(:,2) IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in recv_init" END IF IF(rightproc .gt. mpirank) THEN ! Start to send completed=.FALSE. DO WHILE(.not. completed) CALL MPI_TEST(requests(1), completed,stats(:,1),ierr) END DO !CALL MPI_WAIT(requests(1),stats(:,1),ierr) IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in send_init" END IF END SUBROUTINE start_persistentcomm SUBROUTINE start_rhscomm(mpirank, leftproc, rightproc, moments, nrank, femorder, zlimright) INTEGER, INTENT(IN):: mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:,:), INTENT(INOUT):: moments INTEGER, INTENT(IN):: nrank(2), femorder(2), zlimright +INTEGER, SAVE:: rhsoverlap_requests(2) = MPI_REQUEST_NULL INTEGER:: ierr INTEGER:: stats(MPI_STATUS_SIZE,2) -LOGICAL:: completed=.false. IF(rightproc .gt. mpirank) THEN CALL MPI_ISEND(moments(1,zlimright+1), femorder(1), rhsoverlap_type, rightproc, rhsoverlap_tag, & & MPI_COMM_WORLD, rhsoverlap_requests(1), ierr ) END IF ! If the processor on the left has actually lower z positions IF(leftproc .lt. mpirank) THEN - CALL MPI_IRECV(rhsoverlap_buffer, nrank(2)*(femorder(1)), MPI_DOUBLE_PRECISION, leftproc, rhsoverlap_tag, & + CALL MPI_IRECV(rhsoverlap_buffer, nrank(2)*(femorder(1)), db_type, leftproc, rhsoverlap_tag, & & MPI_COMM_WORLD, rhsoverlap_requests(2), ierr ) END IF CALL MPI_WAITALL(2,rhsoverlap_requests,stats, ierr) END SUBROUTINE start_rhscomm SUBROUTINE start_momentscomm(mpirank, leftproc, rightproc, moments, nrank, femorder, zlimright) INTEGER, INTENT(IN):: mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:,:), INTENT(INOUT):: moments INTEGER, INTENT(IN):: nrank(2), femorder(2), zlimright +INTEGER, SAVE:: momentsoverlap_requests(2) = MPI_REQUEST_NULL INTEGER:: ierr INTEGER:: stats(MPI_STATUS_SIZE,2) IF(rightproc .gt. mpirank) THEN CALL MPI_ISEND(moments(1,zlimright+1), femorder(1), momentsoverlap_type, rightproc, momentsoverlap_tag, & & MPI_COMM_WORLD, momentsoverlap_requests(1), ierr ) END IF ! If the processor on the left has actually lower z positions IF(leftproc .lt. mpirank) THEN - CALL MPI_IRECV(momentsoverlap_buffer, 10*nrank(2)*(femorder(1)), MPI_DOUBLE_PRECISION, leftproc, momentsoverlap_tag, & + CALL MPI_IRECV(momentsoverlap_buffer, 10*nrank(2)*(femorder(1)), db_type, leftproc, momentsoverlap_tag, & & MPI_COMM_WORLD, momentsoverlap_requests(2), ierr ) END IF CALL MPI_WAITALL(2,momentsoverlap_requests,stats, ierr) END SUBROUTINE start_momentscomm !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Initialize the particle MPI type used for inter process communications and publish it to !> the process in the communicator ! !--------------------------------------------------------------------------- SUBROUTINE init_particlempi() INTEGER :: nblock = 10 INTEGER:: blocklength(10) INTEGER(kind=MPI_ADDRESS_KIND):: displs(10), displ0 INTEGER, DIMENSION(10) :: types TYPE(particle) :: part INTEGER:: ierr CALL mpi_get_address(part%Rindex, displs(1), ierr) CALL mpi_get_address(part%Zindex, displs(2), ierr) CALL mpi_get_address(part%partindex, displs(3), ierr) types(1:3)=MPI_INTEGER CALL mpi_get_address(part%R, displs(4), ierr) CALL mpi_get_address(part%Z, displs(5), ierr) CALL mpi_get_address(part%THET, displs(6), ierr) CALL mpi_get_address(part%UZ, displs(7), ierr) CALL mpi_get_address(part%UR, displs(8), ierr) CALL mpi_get_address(part%UTHET, displs(9), ierr) CALL mpi_get_address(part%GAMMA, displs(10), ierr) - types(4:10)=MPI_DOUBLE_PRECISION + types(4:10)=db_type blocklength(1:10) = 1 CALL mpi_get_address(part, displ0, ierr) displs=displs-displ0 CALL MPI_Type_create_struct(nblock, blocklength, displs, types, particle_type, ierr) CALL MPI_Type_commit(particle_type,ierr) END SUBROUTINE init_particlempi !--------------------------------------------------------------------------- !> @author Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Initialize the column MPI type used for inter process communications and publish it to !> the processes in the communicator (can be rhs or grid quantities) ! !> @param[in] nr number of elements in the r direction !> @param[in] nz number of elements in the z direction !> @param[in] init_type MPI type of the initial data !> @param[inout] mpi_coltype final type usable in communications !--------------------------------------------------------------------------- SUBROUTINE init_coltypempi(nr, nz, block_size, stride, init_type, mpi_coltype) INTEGER, INTENT(IN) :: nr INTEGER, INTENT(IN) :: nz INTEGER, INTENT(IN) :: block_size INTEGER, INTENT(IN) :: stride INTEGER, INTENT(IN) :: init_type INTEGER, INTENT(OUT) :: mpi_coltype INTEGER :: temp_mpi_coltype, ierr INTEGER(KIND=MPI_ADDRESS_KIND):: init_type_lb, init_type_extent - !(nrank(2), nrank(1), 1, 10, MPI_DOUBLE_PRECISION, rhsoverlap_type) + !(nrank(2), nrank(1), 1, 10, db_type, rhsoverlap_type) ! Create vector type of length nx CALL MPI_TYPE_VECTOR(nr, block_size, stride*block_size*nz, init_type, temp_mpi_coltype, ierr) CALL MPI_TYPE_COMMIT(temp_mpi_coltype, ierr) ! Get the size in bytes of the initial type CALL MPI_TYPE_GET_EXTENT(init_type, init_type_lb, init_type_extent, ierr) ! Resize temp_mpi_coltype such that the next item to read is at j+1 - CALL MPI_TYPE_CREATE_RESIZED(temp_mpi_coltype, init_type_lb, INT(stride*block_size,kind=MPI_ADDRESS_KIND)*init_type_extent, mpi_coltype, ierr) + CALL MPI_TYPE_CREATE_RESIZED(temp_mpi_coltype, init_type_lb, INT(stride*block_size,kind=MPI_ADDRESS_KIND)*init_type_extent ,& + & mpi_coltype, ierr) CALL MPI_TYPE_COMMIT(mpi_coltype, ierr) END SUBROUTINE init_coltypempi !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Initialize the basicdata MPI type used for inter process communications and publish it to !> the process in the communicator ! !--------------------------------------------------------------------------- SUBROUTINE init_basicdatampi() INTEGER :: nblock = 38 INTEGER:: blocklength(38) INTEGER(kind=MPI_ADDRESS_KIND):: displs(38), displ0 INTEGER, DIMENSION(38) :: types TYPE(BASICDATA) :: bdata INTEGER:: ierr CALL mpi_get_address(bdata%nlres, displs(1), ierr) CALL mpi_get_address(bdata%nlsave, displs(2), ierr) CALL mpi_get_address(bdata%newres, displs(3), ierr) CALL mpi_get_address(bdata%nlxg, displs(4), ierr) CALL mpi_get_address(bdata%nlppform, displs(5), ierr) CALL mpi_get_address(bdata%nlPhis, displs(6), ierr) CALL mpi_get_address(bdata%nlclassical, displs(7), ierr) types(1:7)=MPI_LOGICAL CALL mpi_get_address(bdata%nplasma, displs(8), ierr) CALL mpi_get_address(bdata%nz, displs(9), ierr) CALL mpi_get_address(bdata%it0d, displs(10), ierr) CALL mpi_get_address(bdata%it2d, displs(11), ierr) CALL mpi_get_address(bdata%itparts, displs(12), ierr) CALL mpi_get_address(bdata%itgraph, displs(13), ierr) CALL mpi_get_address(bdata%distribtype, displs(14), ierr) CALL mpi_get_address(bdata%nblock, displs(15), ierr) CALL mpi_get_address(bdata%nrun, displs(16), ierr) types(8:16)=MPI_INTEGER CALL mpi_get_address(bdata%job_time, displs(17), ierr) CALL mpi_get_address(bdata%extra_time, displs(18), ierr) CALL mpi_get_address(bdata%tmax, displs(19), ierr) CALL mpi_get_address(bdata%dt, displs(20), ierr) CALL mpi_get_address(bdata%potinn, displs(21), ierr) CALL mpi_get_address(bdata%potout, displs(22), ierr) CALL mpi_get_address(bdata%B0, displs(23), ierr) CALL mpi_get_address(bdata%n0, displs(24), ierr) CALL mpi_get_address(bdata%temp, displs(25), ierr) CALL mpi_get_address(bdata%Rcurv, displs(26), ierr) CALL mpi_get_address(bdata%width, displs(27), ierr) CALL mpi_get_address(bdata%H0, displs(28), ierr) CALL mpi_get_address(bdata%P0, displs(29), ierr) - types(17:29)=MPI_DOUBLE_PRECISION + types(17:29)=db_type blocklength(1:29) = 1 CALL mpi_get_address(bdata%femorder, displs(30), ierr) CALL mpi_get_address(bdata%ngauss, displs(31), ierr) CALL mpi_get_address(bdata%nnr, displs(32), ierr) types(30:32)=MPI_INTEGER blocklength(30:31) = 2 blocklength(32) = 3 CALL mpi_get_address(bdata%lz, displs(33), ierr) blocklength(33) = 2 CALL mpi_get_address(bdata%radii, displs(34), ierr) CALL mpi_get_address(bdata%plasmadim, displs(35), ierr) - types(33:35)=MPI_DOUBLE_PRECISION + types(33:35)=db_type blocklength(34:35) = 4 CALL mpi_get_address(bdata%resfile, displs(36), ierr) types(36)=MPI_CHARACTER blocklength(36) = 64 CALL mpi_get_address(bdata%partperiodic, displs(37), ierr) types(37)=MPI_LOGICAL blocklength(37) = 1 CALL mpi_get_address(bdata%samplefactor, displs(38), ierr) types(38)=MPI_INTEGER blocklength(38) = 1 CALL mpi_get_address(bdata, displ0, ierr) displs=displs-displ0 CALL MPI_Type_create_struct(nblock, blocklength, displs, types, basicdata_type, ierr) CALL MPI_Type_commit(basicdata_type,ierr) END SUBROUTINE init_basicdatampi END MODULE mpihelper diff --git a/src/resume.f90 b/src/resume.f90 index f6437d6..6e35eb8 100644 --- a/src/resume.f90 +++ b/src/resume.f90 @@ -1,97 +1,98 @@ SUBROUTINE resume ! ! Resume from previous run ! Use beam Use basic Use fields Use sort Use maxwellsource IMPLICIT NONE ! ! Local vars and arrays INTEGER, DIMENSION(:), ALLOCATABLE:: displs,sendcounts INTEGER:: i !________________________________________________________________________________ WRITE(*,'(a)') ' Resume from previous run' !________________________________________________________________________________ ! ! Open and read initial conditions from restart file CALL chkrst(0) IF (newres) THEN cstep=0 time=0 END IF IF(mpisize .gt. 1) THEN ALLOCATE(sendcounts(0:mpisize-1), displs(0:mpisize-1)) IF(mpirank.eq.0) THEN CALL bound(parts) CALL localisation(parts) CALL quicksortparts(parts,1,parts%Nptot) CALL distribpartsonprocs - CALL MPI_Bcast(size(parts%Z),1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(size(parts%Z,1),1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) WRITE(*,*) "nplocs_all: ", Nplocs_all - CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER,MPI_IN_PLACE, 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(Nplocs_all, mpisize, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) DO i=0,mpisize-1 displs(i)=i sendcounts(i)=2 END DO CALL MPI_Scatterv(Zbounds,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,2,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) displs(0)=0 sendcounts(0)=Nplocs_all(0) DO i=1,mpisize-1 displs(i)=displs(i-1)+Nplocs_all(i-1) sendcounts(i)=Nplocs_all(i) END DO CALL MPI_BARRIER(MPI_COMM_WORLD, ierr) - CALL MPI_Scatterv(parts%R,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%Z,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%THET,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%UZ,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%UR,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%R,sendcounts, displs,db_type,MPI_IN_PLACE,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%Z,sendcounts, displs,db_type,MPI_IN_PLACE,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%THET,sendcounts, displs,db_type,MPI_IN_PLACE,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%UZ,sendcounts, displs,db_type,MPI_IN_PLACE,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%UR,sendcounts, displs,db_type,MPI_IN_PLACE,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,db_type,MPI_IN_PLACE,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Rindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Zindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%partindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) ELSE CALL MPI_Bcast(nplasma,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) ! nplasma used temporarily to store required size for parts CALL creat_parts(parts,nplasma) - CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER, Nplocs_all(mpirank), 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(Nplocs_all, mpisize, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) sendcounts=2 CALL MPI_Scatterv(Zbounds,sendcounts,displs,MPI_INTEGER,Zbounds(mpirank),2,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) parts%Nptot=Nplocs_all(mpirank) WRITE(*,*) mpirank, "received Nplocs_all : ", Nplocs_all WRITE(*,*) mpirank, "received Zbounds : ", Zbounds CALL MPI_BARRIER(MPI_COMM_WORLD, ierr) - CALL MPI_Scatterv(parts%R,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%R,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%Z,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%Z,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%THET,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%THET,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%UZ,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UZ,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%UR,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UR,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UTHET,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%R,sendcounts, displs,db_type,parts%R,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%Z,sendcounts, displs,db_type,parts%Z,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%THET,sendcounts, displs,db_type,parts%THET,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%UZ,sendcounts, displs,db_type,parts%UZ,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%UR,sendcounts, displs,db_type,parts%UR,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,db_type,parts%UTHET,Nplocs_all(mpirank),db_type,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Rindex,sendcounts, displs,MPI_INTEGER,parts%Rindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Zindex,sendcounts, displs,MPI_INTEGER,parts%Zindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%partindex,sendcounts, displs,MPI_INTEGER,parts%partindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + parts%Nptot=Nplocs_all(mpirank) END IF - CALL MPI_Bcast(nplasma,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Bcast(V,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Bcast(qsim,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) - CALL MPI_Bcast(msim,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(nplasma,1,db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(V,1,db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(qsim,1,db_type,0,MPI_COMM_WORLD,ierr) + CALL MPI_Bcast(msim,1,db_type,0,MPI_COMM_WORLD,ierr) END IF ! Compute Self consistent electric field CALL bound(parts) CALL localisation(parts) ! Init Electric and Magnetic Fields CALL fields_init IF (nlmaxwellsource) CALL maxwellsource_init(lu_in, time) CALL rhscon(parts) CALL poisson CALL EForcescomp(parts) IF (newres) THEN CALL diagnostics END IF ! END SUBROUTINE resume diff --git a/src/tesend.f90 b/src/tesend.f90 index 24328ad..2a9352c 100644 --- a/src/tesend.f90 +++ b/src/tesend.f90 @@ -1,74 +1,74 @@ SUBROUTINE tesend ! ! Test for run completion ! USE basic ! IMPLICIT NONE ! ! Local vars and arrays CHARACTER(len=*), PARAMETER :: stop_file = 'mystop' LOGICAL :: mlexist REAL(kind=db) :: eltime, step_time, tremain INTEGER :: rem_steps !________________________________________________________________________________ !!$ WRITE(*,'(a/)') '=== Test for run completion ===' !________________________________________________________________________________ ! 1. Some processors had set nlend ! IF( nlend ) THEN WRITE(*,'(/a)') 'NLEND set to .TRUE.!' RETURN END IF !________________________________________________________________________________ ! 2. NRUN modified through "stop file" ! INQUIRE(file=stop_file, exist=mlexist) IF( mlexist ) THEN IF(mpirank .eq. 0) THEN OPEN(lu_stop, file=stop_file) READ(lu_stop,*) rem_steps ! Modify remaining steps "on the fly" CALL MPI_BARRIER(MPI_COMM_WORLD, ierr) - CALL MPI_Bcast(rem_steps, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Bcast(rem_steps, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) WRITE(*,'(/a,i4,a)') 'Stop file found: will exit in', rem_steps, ' steps' CLOSE(lu_stop, status='delete') ELSE CALL MPI_BARRIER(MPI_COMM_WORLD, ierr) - CALL MPI_Bcast(rem_steps, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Bcast(rem_steps, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) END IF nrun = step + rem_steps END IF !________________________________________________________________________________ ! 3. Test on NRUN ! nlend = step .GE. nrun IF ( nlend ) THEN WRITE(*,'(/a)') 'NRUN steps done' RETURN END IF !________________________________________________________________________________ ! 4. Test on TMAX ! nlend = time .GE. tmax IF ( nlend ) THEN WRITE(*,'(/a)') 'TMAX reached' RETURN END IF !________________________________________________________________________________ ! 5. Test on time allocated to job ! CALL timera(-1, '', eltime) ! Current elapsed time tremain = job_time - eltime ! CALL timera(1, 'Main loop', step_time) step_time = 1.2 * step_time / step ! Averaged time per step + 20% ! nlend = tremain .LT. (step_time+extra_time) IF( nlend ) THEN WRITE(*,'(/a,f8.3)') 'Allocated Job time exhausted:, remaining time =', tremain END IF RETURN ! nlend = .FALSE. !________________________________________________________________________________ END SUBROUTINE tesend