diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index 64da389..aef01ca 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,339 +1,340 @@ 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=128) :: 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=100 !< Number of iterations between 2d values writes to hdf5 INTEGER :: itparts=1000 !< Number of iterations between particles values writes to hdf5 INTEGER :: ittext=10 !< Number of iterations between text outputs in the console INTEGER :: itrestart=10000 !< Number of iterations between save of restart.h5 file INTEGER :: ittracer=100 !< Number of iterations between save of traced particles position and velocity INTEGER :: itcelldiag=100000 !< Number of iterations between save of celldiag diagnostic INTEGER :: nbcelldiag=0 !< Number of celldiagnostics 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 INTEGER :: lu_partfile = 120 !< particle loading file, see beam::loadpartfile ! ! HDF5 file CHARACTER(len=256) :: resfile = "results.h5" !< Main result file CHARACTER(len=256) :: rstfile = "restart.h5" !< Restart file CHARACTER(len=256) :: magnetfile = "" !< H5 file containing the magnetic field definition CHARACTER(len=256) :: partfile(10)="" !< Particle loading 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 :: nlfreezephi= .FALSE. !< Freeze the Poisson solver to the field obtained at (re-)start 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 :: nbspecies = 1 !< Number of particles species also counting tracing particles 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), SAVE :: splrz !< Spline at r and z for total electric field TYPE(spline2d), SAVE :: splrz_ext !< Spline at r and z for external electric field 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), 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 REAL(kind=db) :: partmass=me !< Mass of physical particle 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) :: temp !< Initial temperature of plasma REAL(kind=db) :: Rcurv !< Magnetic field curvature coefficient REAL(kind=db) :: Width !< Distance between two magnetic mirrors REAL(kind=db) :: weights_scale=1.0 !< Scale factor for the particle weights on restart (only for newres=.true.) INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction REAL(kind=db):: invdz, invdr(3) 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) :: inputfilename INTEGER, EXTERNAL :: OMP_GET_MAX_THREADS ! 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, partfile, partmass, nbspecies, & - & ittracer, itcelldiag, nbcelldiag, magnetfile, weights_scale + & ittracer, itcelldiag, nbcelldiag, magnetfile, weights_scale, nlfreezephi !________________________________________________________________________________ ! 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) #if _DEBUG==1 WRITE(*,*) "Compiled in debug mode" #endif 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 WRITE(*,'(a,i4.2,a,i4.2,a)')"Running on ",mpisize," tasks with", omp_get_max_threads() ," openMP threads" 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(funcstr) CHARACTER(len=*), INTENT(in) :: funcstr INTEGER :: j, ind loc = 0 DO j=1,nc ind = INDEX(which(j), funcstr(1:lstr)) IF( ind .GT. 0 .AND. LEN_TRIM(which(j)) .EQ. lstr ) THEN loc = j 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 f34f9b3..a02596f 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,2098 +1,2107 @@ !------------------------------------------------------------------------------ ! 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 :: Nploc !< Local number of simulated particles INTEGER :: Nptot !< Total number of simulated particles INTEGER :: Newindex !< Stores the higher partindex for the creation of new particles REAL(kind=db) :: m !< Particle mass REAL(kind=db) :: q !< Particle charge REAL(kind=db) :: weight !< Number of particles represented by one macro-particle REAL(kind=db) :: qmRatio !< Charge over mass ratio REAL(kind=db) :: H0 REAL(kind=db) :: P0 REAL(kind=db) :: temperature LOGICAL :: Davidson=.false. LOGICAL :: is_test= .false. INTEGER, DIMENSION(4) :: nblost !< number of particles lost in (z_min,z_max,r_min,r_max) since last gather INTEGER :: nbadded !< number of particles added by source since last gather 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 :: BZ !< axial radial relative distances to the left grid line REAL(kind=db), DIMENSION(:), ALLOCATABLE :: BR !< radial relative distances to the bottom grid line REAL(kind=db), DIMENSION(:), ALLOCATABLE :: pot !< Electric potential REAL(kind=db), DIMENSION(:), ALLOCATABLE :: potxt !< External 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 Real(kind=db), Dimension(:,:),ALLOCATABLE:: geomweight !< geometric weight at the particle position LOGICAL:: collected !< Stores if the particles data have been collected to MPI root process during this timestep INTEGER, DIMENSION(:), ALLOCATABLE:: addedlist END TYPE particles ! !TYPE(particles) :: parts !< Storage for all the particles !SAVE :: parts TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist ! Diagnostics (scalars) REAL(kind=db) :: ekin=0 !< Total kinetic energy (J) REAL(kind=db) :: epot=0 !< Total potential energy (J) REAL(kind=db) :: etot=0 !< Current total energy (J) REAL(kind=db) :: etot0=0 !< Initial total energy (J) REAL(kind=db) :: loc_etot0=0 !< theoretical local total energy (J) REAL(kind=db) :: Energies(4) !< (1) kinetic energy, (2) potential energy, (3) total energy and (4) gained/lossed energy due to gain or loss of particles (J) ! INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process ! 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%Nploc = nparts p%Nptot = nparts ALLOCATE(p%Z(nparts)) ALLOCATE(p%R(nparts)) ALLOCATE(p%THET(nparts)) ALLOCATE(p%BZ(nparts)) ALLOCATE(p%BR(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%potxt(nparts)) ALLOCATE(p%Er(nparts)) ALLOCATE(p%Ez(nparts)) ALLOCATE(p%GAMMAold(nparts)) Allocate(p%geomweight(nparts,0:2)) p%newindex=0 p%nblost=0 p%nbadded=0 p%partindex=-1 p%URold=0 p%UZold=0 p%UTHETold=0 p%rindex=0 p%zindex=0 p%BR=0 p%BZ=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%potxt=0 p%gammaold=1 p%collected=.false. p%Davidson=.false. p%is_test=.false. p%m=me p%q=-elchar p%qmRatio=p%q/p%m p%weight=1.0_db p%H0=0 p%P0=0 p%temperature=0 p%geomweight=0 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, nbspecies, Zbounds USE mpi INTEGER:: i, j REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V, uniform in Z and 1/R in R CALL loaduniformRZ(partslist(1), VR, VZ, VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(3) !Stable distribution from Davidson 4.95 p.119 but with constant distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodinvr) CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodgausr) CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz @brief Checks for each particle if the z position is outside of the local/global simulation space. !> Depending on the boundary conditions, the leaving particles are sent to the correct neighbouring MPI process !> or deleted. ! !> @param[in] p particles structure ! !> @author Guillaume Le Bars EPFL/SPC !--------------------------------------------------------------------------- SUBROUTINE bound(p) - USE basic, ONLY: zgrid, nz, Zbounds, mpirank, step, leftproc, rightproc + USE basic, ONLY: zgrid, nz, Zbounds, mpirank, step, leftproc, rightproc, partperiodic USE IFPORT IMPLICIT NONE type(particles), INTENT(INOUT):: p INTEGER :: i, rsendnbparts, lsendnbparts, nblostparts INTEGER :: receivednbparts, partdiff, lostindex, sentindex INTEGER, DIMENSION(p%Nploc) :: sendhole INTEGER, DIMENSION(p%Nploc) :: losthole LOGICAL:: leftcomm, rightcomm, sent INTEGER, ALLOCATABLE:: partstoremove(:) rsendnbparts=0 lsendnbparts=0 nblostparts=0 losthole=0 sendhole=0 receivednbparts=0 ! We communicate with the left processus leftcomm = leftproc .ne. -1 ! We communicate with the right processus rightcomm = rightproc .ne. -1 IF (p%Nploc .gt. 0) THEN ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) DO i=1,p%Nploc ! If the particle is to the right of the local simulation space, it is sent to the right MPI process IF (p%Z(i) .ge. zgrid(Zbounds(mpirank+1))) THEN + IF(partperiodic) THEN + DO WHILE (p%Z(i) .GT. zgrid(nz)) + p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0) + END DO + END IF !$OMP CRITICAL (nbparts) IF(rightcomm) THEN rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i - DO WHILE (p%Z(i) .GT. zgrid(nz)) - p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0) - END DO - ELSE + ELSE if(.not. partperiodic) THEN nblostparts=nblostparts+1 losthole(nblostparts)=i p%nblost(2)=p%nblost(2)+1 END IF !$OMP END CRITICAL (nbparts) ! If the particle is to the left of the local simulation space, it is sent to the left MPI process ELSE IF (p%Z(i) .lt. zgrid(Zbounds(mpirank))) THEN + IF(partperiodic) THEN + DO WHILE (p%Z(i) .LT. zgrid(0)) + p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0) + END DO + END IF !$OMP CRITICAL (nbparts) IF(leftcomm) THEN ! We send the particle to the left process lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i - DO WHILE (p%Z(i) .LT. zgrid(0)) - p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0) - END DO - ELSE + ELSE if(.not. partperiodic) THEN ! we destroy the particle nblostparts=nblostparts+1 losthole(nblostparts)=i p%nblost(1)=p%nblost(1)+1 END IF !$OMP END CRITICAL (nbparts) END IF END DO !$OMP END PARALLEL DO END IF IF(mpisize .gt. 1) THEN ! We send the particles leaving the local simulation space to the closest neighbour CALL particlescommunication(p, lsendnbparts, rsendnbparts, sendhole, receivednbparts, (/leftproc,rightproc/)) END IF ! If the boundary conditions are not periodic, we delete the corresponding particles IF(nblostparts .gt. 0 .and. step .ne. 0) THEN DO i=1,nblostparts CALL delete_part(p, losthole(i), .false. ) END DO !WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank END IF ! computes if we received less particles than we sent partdiff=max(lsendnbparts+rsendnbparts-receivednbparts,0) IF(nblostparts + partdiff .gt. 0) THEN ALLOCATE(partstoremove(nblostparts+partdiff)) partstoremove(1:partdiff)=abs(sendhole(receivednbparts+1:receivednbparts+partdiff)) partstoremove(partdiff+1:partdiff+nblostparts)=abs(losthole(1:nblostparts)) call qsort(partstoremove,size(partstoremove),sizeof(partstoremove(1)),compare_int) ! If we received less particles than we sent, or lost particles we fill the remaining holes with the particles from the end of the ! parts arrays DO i=1,nblostparts+partdiff CALL move_part(p, p%Nploc, partstoremove(i)) p%partindex(p%Nploc)=-1 p%Nploc = p%Nploc-1 END DO END IF CONTAINS INTEGER(2) function compare_int(a,b) INTEGER :: a, b, c c=b-a if(c.ne.0) c=1 compare_int=sign(c,b-a) end function 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. !> @param[in] p particles structure !--------------------------------------------------------------------------- SUBROUTINE localisation(p) USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, mpirank Use geometry, ONLY: r_a, geom_weight USE IFPORT type(particles), INTENT(INOUT):: p INTEGER :: i, nblostparts, iend,nbunch INTEGER, DIMENSION(p%Nploc) :: losthole INTEGER, DIMENSION(2):: nblost losthole=0 nblostparts=0 nblost=0 nbunch=64 IF (p%Nploc .gt. 0) THEN !$OMP PARALLEL DEFAULT(SHARED), private(i,iend) !$OMP DO DO i=1,p%Nploc,nbunch ! Avoid segmentation fault by accessing non relevant data iend=min(i+nbunch-1,p%Nploc) call geom_weight(p%Z(i:iend), p%r(i:iend), p%geomweight(i:iend,:)) END DO !$OMP END DO !$OMP DO reduction(+:nblost,nblostparts) DO i=1,p%Nploc if(p%geomweight(i,0).le.0 .or. p%R(i) .ge. rgrid(nr) .or. p%R(i) .le. rgrid(0)) then ! If the particle is outside of the simulation space in the r direction, it is deleted. !!$OMP CRITICAL (lostparts) nblostparts=nblostparts+1 losthole(i)=i !!$OMP END CRITICAL (lostparts) if(p%R(i) .le. r_a) then nblost(1)=nblost(1)+1 else nblost(2)=nblost(2)+1 end if else call p_calc_rzindex(p,i) end if END DO !$OMP END DO !$OMP END PARALLEL IF(nblostparts.gt.0) THEN p%nblost(3:4)=nblost+p%nblost(3:4) call qsort(losthole,p%Nploc,sizeof(losthole(1)),compare_int) DO i=1,nblostparts CALL delete_part(p,losthole(i)) END DO END IF END IF CONTAINS INTEGER(2) function compare_int(a,b) INTEGER :: a, b, c c=b-a if(c.ne.0) c=1 compare_int=sign(c,b-a) end function END SUBROUTINE localisation subroutine p_calc_rzindex(p,i) use basic, only: rgrid,zgrid,invdz,invdr, nnr, nr integer::i type(particles)::p 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) End if p%zindex(i)=(p%Z(i)-zgrid(0))*invdz end subroutine p_calc_rzindex SUBROUTINE comp_mag_p(p) USE basic, ONLY: zgrid, rgrid, dz, BZ, BR, nz, invdz type(particles), INTENT(INOUT):: p INTEGER :: i, nblostparts Real(kind=db):: WZ,WR INTEGER:: j1,j2,j3,j4 !$OMP PARALLEL DO SIMD DEFAULT(SHARED) Private(J1,J2,J3,J4,WZ,WR) DO i=1,p%Nploc WZ=(p%Z(i)-zgrid(p%zindex(i)))*invdz; WR=(p%R(i)-rgrid(p%rindex(i)))/(rgrid(p%rindex(i)+1)-rgrid(p%rindex(i))); J1=(p%rindex(i))*(nz+1) + p%zindex(i)+1 J2=(p%rindex(i))*(nz+1) + p%zindex(i)+2 J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1 J4=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+2 ! Interpolation for magnetic field p%BZ(i)=(1-WZ)*(1-WR)*Bz(J4) & & +WZ*(1-WR)*Bz(J3) & & +(1-WZ)*WR*Bz(J2) & & +WZ*WR*Bz(J1) p%BR(i)=(1-WZ)*(1-WR)*Br(J4) & & +WZ*(1-WR)*Br(J3) & & +(1-WZ)*WR*Br(J2) & & +WZ*WR*Br(J1) END DO !$OMP END PARALLEL DO SIMD end subroutine comp_mag_p !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief General routine to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity(p) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : nlclassical type(particles), INTENT(INOUT):: p ! Store old Velocities CALL swappointer(p%UZold, p%UZ) CALL swappointer(p%URold, p%UR) CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (nlclassical) THEN CALL comp_velocity_fun(p, gamma_classical) ELSE CALL comp_velocity_fun(p, gamma_relativistic) END IF END SUBROUTINE comp_velocity !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine called by comp_velocity to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using the routine computing gamma as an input. This avoid the nlclassical flag check on each particle. ! !> @param[in] gamma the function used to compute the value of the lorentz factor \f$\gamma\f$ !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity_fun(p, gamma) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : bnorm, dt, BZ, BR, nz interface REAL(kind=db) FUNCTION gamma(UZ, UR, UTHET) USE constants REAL(kind=db), INTENT(IN):: UR,UZ,UTHET end FUNCTION 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 tau=p%qmRatio*bnorm*0.5*dt IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,p%Nploc !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 p%Gamma(i)=gamma(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) ZBZ=tau*p%BZ(i)/p%Gamma(i) ZBR=tau*p%BR(i)/p%Gamma(i) ZPZ=p%UZ(i)-ZBR*p%UTHETold(i) !u'_{z} ZPR=p%UR(i)+ZBZ*p%UTHETold(i) !u'_{r} ZPTHET=p%UTHETold(i)+(ZBR*p%UZ(i)-ZBZ*p%UR(i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR p%UZ(i)=p%UZ(i)-ZBR2*ZPTHET !u+_{z} p%UR(i)=p%UR(i)+ZBZ2*ZPTHET !u+_{r} p%UTHET(i)=p%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration p%UZ(i)=p%UZ(i)+p%EZ(i)*tau p%UR(i)=p%UR(i)+p%ER(i)*tau ! Final computation of the Lorentz factor p%Gamma(i)=gamma(p%UZ(i), p%UR(i), p%UTHET(i)) END DO !$OMP END PARALLEL DO SIMD END IF p%collected=.false. END SUBROUTINE comp_velocity_fun !--------------------------------------------------------------------------- !> @author !> 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 REAL(kind=db) FUNCTION gamma_classical(UZ, UR, UTHET) #if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_classical) #endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_classical=1.0 END FUNCTION gamma_classical !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the relativistic simulations. !> This routine computes the Lorentz factor \f$\gamma=\sqrt{1+\mathbf{\gamma\beta}^2}\f$ ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- ELEMENTAL REAL(kind=db) FUNCTION gamma_relativistic(UZ, UR, UTHET) #if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_relativistic) #endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_relativistic=sqrt(1+UZ**2+UR**2+UTHET**2) END FUNCTION gamma_relativistic !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles position at time t+1 !> This routine computes the particles position at time t+1 according to the Bunemann algorithm. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE push(p) Use basic, ONLY: dt, tnorm type(particles), INTENT(INOUT):: p REAL(kind=db):: XP, YP, COSA, SINA, U1, U2, ALPHA INTEGER :: i IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(XP, YP, COSA, SINA, U1, U2, ALPHA) DO i=1,p%Nploc ! Local Cartesian coordinates XP=p%R(i)+dt/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 p%collected=.false. END SUBROUTINE push !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes several diagnostic quantities !> This routine computes the total kinetic and electric potential energy. !> It keeps track of the reference energy and the number of particle per mpi node. ! !--------------------------------------------------------------------------- SUBROUTINE diagnostics ! ! Compute energies ! USE constants, ONLY: vlight USE basic, ONLY: phinorm, cstep, nlclassical, ierr, step, nlend,& & itparts INTEGER:: i ! 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,partslist(1)%Nploc ! Potential energy epot=epot+(partslist(1)%pot(i)+partslist(1)%potxt(i)) ! Kinetic energy IF(.not. nlclassical) THEN ekin=ekin+(partslist(1)%Gamma(i)-1) ELSE ekin=ekin+0.5*( partslist(1)%UR(i)**2 & & + partslist(1)%UZ(i)**2 & & + partslist(1)%UTHET(i)**2 ) END IF END DO !$OMP END PARALLEL DO SIMD epot=epot*phinorm*0.5*partslist(1)%q*partslist(1)%weight ekin=ekin*partslist(1)%m*partslist(1)%weight*vlight**2 ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ. 0) THEN ! Compute the local total energy loc_etot0 = epot+ekin etot0=0 END IF !etot=loc_etot0 ! Compute the total energy etot=epot+ekin Energies=(/ekin,epot,etot,loc_etot0/) ! The computed energy is sent to the root process IF(mpisize .gt.1) THEN IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) etot0=etot0+Energies(4) ekin=Energies(1) epot=Energies(2) etot=Energies(3) ELSE CALL MPI_REDUCE(Energies, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) END IF ELSE etot0=etot0+loc_etot0 END IF loc_etot0=0 ! Send the local number of particles on each node to the root process IF(mpisize .gt. 1) THEN Nplocs_all(mpirank)=partslist(1)%Nploc IF(mpirank .eq.0 ) THEN CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) partslist(1)%Nptot=sum(Nplocs_all) ELSE CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) END IF ELSE partslist(1)%Nptot=partslist(1)%Nploc 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(partslist(1)) 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, Nploc INTEGER:: i IF(p%collected) RETURN ! exit subroutine if particles have already been collected during this time step ALLOCATE(Nploc(0:mpisize-1)) ALLOCATE(displs(0:mpisize-1)) displs=0 Nploc(mpirank)=p%Nploc CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, Nploc, 1, MPI_INTEGER,& & MPI_COMM_WORLD, ierr) p%Nptot=sum(Nploc) IF(p%Nptot .eq. 0 ) THEN p%collected=.true. RETURN END IF IF(mpirank .ne. 0) THEN ! Send Particles informations to root process CALL MPI_Gatherv(p%Z, Nploc(mpirank), db_type, p%Z, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%R, Nploc(mpirank), db_type, p%R, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%THET, Nploc(mpirank), db_type, p%THET, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%UR, Nploc(mpirank), db_type, p%UR, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%UZ, Nploc(mpirank), db_type, p%UZ, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%UTHET, Nploc(mpirank), db_type, p%UTHET, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%pot, Nploc(mpirank), db_type, p%pot, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%Rindex, Nploc(mpirank), MPI_INTEGER, p%Rindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%Zindex, Nploc(mpirank), MPI_INTEGER, p%Zindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%partindex, Nploc(mpirank), MPI_INTEGER, p%partindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) ELSE Do i=1,mpisize-1 displs(i)=displs(i-1)+Nploc(i-1) END DO IF(p%Nptot .gt. size(p%R,1)) THEN CALL change_parts_allocation(p,p%Nptot-size(P%R,1)) END IF ! Receive particle information from all processes CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%Z, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%R, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%THET, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%UR, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%UZ, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%UTHET, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%pot, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), MPI_INTEGER, p%Rindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), MPI_INTEGER, p%Zindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), MPI_INTEGER, p%partindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) p%partindex(sum(Nploc)+1:)=-1 END IF p%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 : bnorm, dt, BZ, BR, nlclassical, phinorm, nz, distribtype, vnorm type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, Vperp, v2 INTEGER :: J1, J2, J3, J4, i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ! In case Davidson distribution is used the longitudinal and radial velocities are adapted to take into account the ! electric potential. IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4 .or. p%Davidson) THEN ALLOCATE(VR(p%Nploc),VZ(p%Nploc),VTHET(p%Nploc)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=p%UTHET/p%Gamma*vnorm DO i=1,p%Nploc Vperp=sqrt(MAX(2*p%H0/p%m-2*p%qmRatio*p%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN p%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 p%Gamma(i)=sqrt(1/(1-v2/vnorm**2)) END IF p%UR(i)=p%Gamma(i)*VR(i)/vnorm p%UZ(i)=p%Gamma(i)*VZ(i)/vnorm p%UTHET(i)=p%Gamma(i)*VTHET(i)/vnorm END DO DEALLOCATE(VR,VZ,VTHET) END IF RETURN ! Normalised time increment !tau=-omegac/2/omegap*dt/tnorm tau=p%qmRatio*bnorm*0.5*dt ! Store old Velocities CALL swappointer(p%UZold, p%UZ) CALL swappointer(p%URold, p%UR) CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR) DO i=1,p%Nploc ! 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%BZ(i))*(1-p%BR(i))*Bz(J4) & !& +p%BZ(i)*(1-p%BR(i))*Bz(J3) & !& +(1-p%BZ(i))*p%BR(i)*Bz(J2) & !& +p%BZ(i)*p%BR(i)*Bz(J1) !BRR=(1-p%BZ(i))*(1-p%BR(i))*Br(J4) & !& +p%BZ(i)*(1-p%BR(i))*Br(J3) & !& +(1-p%BZ(i))*p%BR(i)*Br(J2) & !& +p%BZ(i)*p%BR(i)*Br(J1) ! Half inverse Rotation along magnetic field ZBZ=tau*p%BZ(i)/p%Gammaold(i) ZBR=tau*p%BR(i)/p%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(p%UZold(i)-ZBR*p%UTHETold(i))/SQR !u-_{z} ZPR=(p%URold(i)+ZBZ*p%UTHETold(i))/SQR !u-_{r} ZPTHET=p%UTHETold(i)+(ZBR*p%UZold(i)-ZBZ*p%URold(i))/SQR !u-_{theta} p%UZ(i)=ZPZ p%UR(i)=ZPR p%UTHET(i)=ZPTHET ! half of decceleration p%UZ(i)=p%UZ(i)+p%Ez(i)*tau p%UR(i)=p%UR(i)+p%Er(i)*tau IF(.not. nlclassical) THEN p%Gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2) END IF END DO !$OMP END PARALLEL DO SIMD END IF END SUBROUTINE adapt_vinit !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief In the case of MPI parallelism, computes the indices of the particle assigned to the current process. !--------------------------------------------------------------------------- SUBROUTINE distribpartsonprocs(p) ! 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, Zbounds TYPE(particles), INTENT(INOUT):: p INTEGER:: k, i, nbparts REAL(kind=db):: idealnbpartsperproc 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(p%Nploc)/REAL(mpisize)) Zmin=MINVAL(p%Zindex(1:p%Nploc)) Zmax=MAXVAL(p%Zindex(1:p%Nploc)) Zperproc=(Zmax-Zmin)/mpisize IF(Zmax .eq. 0) Zmax=nz DO k=1,p%Nploc zindex=p%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 DO k=0,mpisize-1 Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1)) END DO WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " nptot", Nplocs_all(mpirank) END SUBROUTINE distribpartsonprocs SUBROUTINE keep_mpi_self_parts(p,Zbounds) TYPE(particles),INTENT(INOUT):: p INTEGER,INTENT(in)::Zbounds(0:) INTEGER :: i, partstart, old_sum,ierr partstart=1 p%Nploc=0 Do i=1,p%Nptot IF(p%Zindex(i).ge.Zbounds(mpirank).and.p%Zindex(i).lt.Zbounds(mpirank+1)) THEN p%Nploc=p%Nploc+1 CALL move_part(p,i,p%Nploc) END IF END DO old_sum=p%Nptot CALL MPI_REDUCE(p%Nploc, p%Nptot,1,MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr) IF(p%Nptot .ne. old_sum) WRITE(*,*) "Error in particle distribution kept: ", p%Nptot, "/",old_sum END SUBROUTINE keep_mpi_self_parts !_______________________________________________________________________________ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Manage the particle communication between neighbours. !> This routine is responsible to receive the incoming particles from the MPI neighbours and to send its outgoing !> particles to these neighbours ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative !--------------------------------------------------------------------------- SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, sendholes, receivednbparts, procs) USE mpihelper, ONLY: particle_type #ifdef _DEBUG USE basic, ONLY: step #endif type(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(out) :: receivednbparts INTEGER, INTENT(in) :: sendholes(:) INTEGER, INTENT(in) :: procs(2) INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0 INTEGER, ASYNCHRONOUS :: sendrequest(2), recvrequest(2) INTEGER, ASYNCHRONOUS :: sendstatus(MPI_STATUS_SIZE,2), recvstatus(MPI_STATUS_SIZE,2) TYPE(particle), ALLOCATABLE :: rrecvpartbuff(:), lrecvpartbuff(:), rsendpartbuff(:), lsendpartbuff(:) ! buffers to send and receive particle from left and right processes INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts, ierr lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL lrecvnbparts=0 rrecvnbparts=0 ! Send and receive the number of particles to exchange CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, procs(1), 0, MPI_COMM_WORLD, recvrequest(1), ierr) CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, procs(2), 0, MPI_COMM_WORLD, recvrequest(2), ierr) CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, procs(1), 0, MPI_COMM_WORLD, sendrequest(1), ierr) CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, procs(2), 0, MPI_COMM_WORLD, sendrequest(2), ierr) CALL MPI_Waitall(2,recvrequest(1:2), recvstatus(:,1:2), ierr) recvrequest=MPI_REQUEST_NULL lreceivednbparts=lrecvnbparts rreceivednbparts=rrecvnbparts ! Re/allocate enough memory to store the incoming particles ALLOCATE(rrecvpartbuff(rreceivednbparts)) ALLOCATE(lrecvpartbuff(lreceivednbparts)) ! Receive particles from left and right processes to the corresponding buffers IF ( lrecvnbparts .gt. 0) THEN CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, procs(1), 1, MPI_COMM_WORLD, recvrequest(1), ierr) END IF IF( rrecvnbparts .gt. 0) THEN CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, procs(2), 1, MPI_COMM_WORLD, recvrequest(2), ierr) END IF ALLOCATE(rsendpartbuff(rsendnbparts)) ALLOCATE(lsendpartbuff(lsendnbparts)) ! Copy the leaving particles to the corresponding send buffers IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes, lsendpartbuff, rsendpartbuff) END IF CALL MPI_Waitall(2,sendrequest(1:2), sendstatus(:,1:2), ierr) ! Send the particles to the left and right neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, procs(1), 1, MPI_COMM_WORLD, sendrequest(1), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", lsendnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, procs(2), 1, MPI_COMM_WORLD, sendrequest(2), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", rsendnbparts , " to right at step: ",step #endif END IF ! Receive the incoming parts in the receive buffers IF ( lreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(1), recvstatus(:,1), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,1) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", lreceivednbparts , " from left at step: ",step #endif END IF IF ( rreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(2), recvstatus(:,2), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,2) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", rreceivednbparts , " from right at step: ",step #endif END IF receivednbparts=rreceivednbparts+lreceivednbparts IF(p%Nploc+receivednbparts-lsendnbparts-rsendnbparts .gt. size(p%R,1)) THEN CALL change_parts_allocation(p,receivednbparts) END IF ! Copy the incoming particles from the receive buffers to the simulation parts variable CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, & & sendholes, lrecvpartbuff, rrecvpartbuff) ! Wait for the outgoing particles to be fully received by the neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(1), sendstatus(:,1), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", lsentnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(2), sendstatus(:,2), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", rsentnbparts , " to right at step: ",step #endif END IF ! ! END SUBROUTINE particlescommunication !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the receive buffers to the local simulation variable parts. !> The incoming particles will first be stored in the holes left by the outgoing particles, then they !> will be added at the end of the parts variable ! !> @param [in] rrecvnbparts number of particles received from the right neighbour (mpirank+1) !> @param [in] lrecvnbparts number of particles received from the left neighbour (mpirank-1) !> @param [in] sendnbparts total number of particles having left the local domain !> @param [in] sendholes array containing the indices of the particle having left the local domain in ascending order. !--------------------------------------------------------------------------- SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts, sendholes,lrecvpartbuff, rrecvpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(:) TYPE(particle), INTENT(IN) :: rrecvpartbuff(:), lrecvpartbuff(:) INTEGER k,partpos ! First import the particles coming from the right IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(sendholes(k)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, rrecvpartbuff, partpos, k) END DO END IF ! Then import the particles coming from the left IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(sendholes(k+rrecvnbparts)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, lrecvpartbuff, partpos, k) 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, partsindex, bufferindex) USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), 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 p%pot(partsindex) = buffer(bufferindex)%pot ! 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(:), 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) buffer(bufferindex)%pot = p%pot(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, lsendpartbuff, rsendpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(:) TYPE(particle), INTENT(OUT) :: rsendpartbuff(:), lsendpartbuff(:) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 ! Loop over the outgoing particles and fill the correct send buffer Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(sendholes(k)) IF(sendholes(k) .GT. 0) THEN rsendpos=rsendpos+1 CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos) ELSE IF(sendholes(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(p, lsendpartbuff, lsendpos, partpos) END IF END DO ! ! END SUBROUTINE AddPartSendBuffers !--------------------------------------------------------------------------- !> @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, geomweight(0:2),pot 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) pot = p%pot(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) geomweight = p%geomweight(index1,:) ! Move part at index2 in part at index 1 p%partindex(index1) = p%partindex(index2) p%Gamma(index1) = p%Gamma(index2) p%pot(index1) = p%pot(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) p%geomweight(index1,:) = p%geomweight(index2,:) ! Move temporary values from part(index1) to part(index2) p%partindex(index2) = partindex p%Gamma(index2) = Gamma p%pot(index2) = pot 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 p%geomweight(index2,:) = geomweight 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%potxt,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%BR,sizedifference) CALL change_array_size_dp(p%BZ,sizedifference) CALL change_array_size_dp2(p%geomweight,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%Nploc=MIN(p%Nploc,size(p%R)) 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_dp2(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,1) new_size=current_size+sizedifference ALLOCATE(temp(new_size,0:size(arr,2)-1)) 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_dp2 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) + SUBROUTINE add_created_part(p, buffer,nb_ins) USE bsplines USE basic, ONLY: splrz, phinorm, nlclassical USE geometry IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:) + INTEGER, OPTIONAL:: nb_ins INTEGER:: i, nptotinit, parts_size_increase, nb_added real(kind=db), ALLOCATABLE::gtildeloc(:,:) nptotinit=p%Nploc+1 - nb_added=size(buffer,1) + if(present(nb_ins)) THEN + nb_added=nb_ins + ELSE + nb_added=size(buffer,1) + end if IF(nb_added .le. 0) RETURN ! No particles to add ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nploc + nb_added .gt. size(p%Z,1)) THEN parts_size_increase=Max(floor(0.1*size(p%Z,1)),nb_added) CALL change_parts_allocation(p, parts_size_increase) END IF DO i=1,nb_added p%Nploc=p%Nploc+1 p%newindex=p%newindex+1 CALL Insertincomingpart(p, buffer, p%Nploc, i) p%partindex(p%Nploc)=p%newindex CALL geom_weight(p%Z(p%Nploc),p%R(p%Nploc),p%geomweight(p%Nploc,:)) if( .not. is_inside(p,p%Nploc) ) then p%Nploc=p%Nploc-1 p%newindex=p%newindex-1 CYCLE end if call p_calc_rzindex(p,p%Nploc) END DO nb_added=p%Nploc-nptotinit+1 if( .not. p%is_test) then IF(allocated(p%addedlist)) then call change_array_size_int(p%addedlist,2) else allocate(p%addedlist(2)) end if p%addedlist(size(p%addedlist)-1)=nptotinit p%addedlist(size(p%addedlist))=nb_added end if END SUBROUTINE add_created_part function is_inside(p,id) Use basic, ONLY: rgrid,zgrid, nr, nz IMPLICIT NONE logical :: is_inside type(particles) :: p integer :: id is_inside=.true. if(p%geomweight(id,0).le.0)then is_inside=.false. return end if if(p%R(id).ge.rgrid(nr) .or. p%R(id) .le. rgrid(0))then is_inside=.false. return end if if(p%Z(id).ge.zgrid(nz) .or. p%Z(id) .le. zgrid(0))then is_inside=.false. return end if end function is_inside SUBROUTINE calc_newparts_energy(p) USE basic, ONLY: phinorm, nlclassical type(particles)::p integer::i,n,nptotinit,nbadded, nptotend if(p%is_test) return if( allocated(p%addedlist)) then n=size(p%addedlist) !write(*,*) n, "addedlist: ", p%addedlist Do i=1,n,2 nptotinit=p%addedlist(i) nbadded=p%addedlist(i+1) p%nbadded=p%nbadded+nbadded nptotend=nptotinit+nbadded-1 loc_etot0=loc_etot0+p%q*p%weight*sum(p%pot(nptotinit:nptotend))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0+p%m*p%weight*vlight**2*sum(p%Gamma(nptotinit:nptotend)-1) ELSE loc_etot0=loc_etot0+0.5*p%m*p%weight*vlight**2*sum(p%UR(nptotinit:nptotend)**2 & & +p%UZ(nptotinit:nptotend)**2 & & +p%UTHET(nptotinit:nptotend)**2) END IF end do deallocate(p%addedlist) end if end subroutine calc_newparts_energy !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Delete particle at given index removing its energy from the system ! !> @param [in] index index of particle to be deleted !--------------------------------------------------------------------------- SUBROUTINE delete_part(p, index, replace) !! This will destroy particle at the given index USE constants, ONLY: vlight USE bsplines USE geometry USE basic, ONLY: phinorm, nlclassical,splrz TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index LOGICAL, OPTIONAL :: replace LOGICAL:: repl IF(present(replace)) THEN repl=replace ELSE repl=.true. END IF !Computes the potential at the particle position with phi_ext+phi_s IF(index .le. p%Nploc) THEN IF(.not. p%is_test) THEN loc_etot0=loc_etot0-p%q*p%weight*(p%pot(index))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0-p%m*p%weight*vlight**2*(p%Gamma(index)-1) ELSE loc_etot0=loc_etot0-0.5*p%m*p%weight*vlight**2*(p%UR(index)**2+p%UZ(index)**2+p%UTHET(index)**2) END IF END IF IF(repl) THEN ! We fill the gap CALL move_part(p, p%Nploc, index) p%partindex(p%Nploc)=-1 ! Reduce the total number of simulated parts p%Nploc=p%Nploc-1 END IF END IF END SUBROUTINE delete_part !--------------------------------------------------------------------------- !> @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 IF(sourceindex .eq. destindex) RETURN IF(sourceindex .le. 0 .or. destindex .le. 0) RETURN ! Move part at sourceindex in part at destindex Call copy_part(p,sourceindex,destindex,p) END SUBROUTINE move_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy particle with index sourceindex in particles sourcep to particle with index destindex in particles destp. !> !WARNING! This will overwrite particle at destp(destindex). ! !> @param [inout] sourcep Structure of source particles. !> @param [in] sourceindex index in parts of the particle to move. !> @param [in] destindex index in parts of the moved particle destination. !> @param [inout] destp Structure of source particles. !--------------------------------------------------------------------------- SUBROUTINE copy_part(sourcep, sourceindex, destindex, destp) !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex TYPE(particles), INTENT(IN)::sourcep TYPE(particles), INTENT(INOUT)::destp IF(sourceindex .le. 0 .or. destindex .le. 0) RETURN IF( destindex .gt. size(destp%R,1)) RETURN ! Move part at sourceindex in part at destindex destp%partindex(destindex) = sourcep%partindex(sourceindex) destp%Gamma(destindex) = sourcep%Gamma(sourceindex) destp%R(destindex) = sourcep%R(sourceindex) destp%Z(destindex) = sourcep%Z(sourceindex) destp%THET(destindex) = sourcep%THET(sourceindex) destp%UR(destindex) = sourcep%UR(sourceindex) destp%UTHET(destindex) = sourcep%UTHET(sourceindex) destp%UZ(destindex) = sourcep%UZ(sourceindex) destp%Rindex(destindex) = sourcep%Rindex(sourceindex) destp%Zindex(destindex) = sourcep%Zindex(sourceindex) destp%geomweight(destindex,:) = sourcep%geomweight(sourceindex,:) destp%pot(destindex) = sourcep%pot(sourceindex) END SUBROUTINE copy_part !________________________________________________________________________________ SUBROUTINE destroy_parts(p) TYPE(particles) :: p p%Nploc=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%BZ)) DEALLOCATE(p%BZ) IF(ALLOCATED(p%BR)) DEALLOCATE(p%BR) 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) if(allocated(p%geomweight)) Deallocate(p%geomweight) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE clean_beam ! INTEGER:: i Do i=1,size(partslist,1) CALL destroy_parts(partslist(i)) END DO ! 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(p, VR,VZ,VTHET) USE basic, ONLY: plasmadim, rnorm, temp, qsim, msim USE constants, ONLY: me, kb, elchar REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=size(VR,1) p%q=sign(elchar,qsim) p%weight=msim/me p%m=me p%qmRatio=qsim/msim ! Initial distribution in z with normalisation CALL loduni(1,p%Z(1:p%Nploc)) p%Z(1:p%Nploc)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*p%Z(1:p%Nploc))/rnorm ! Initial distribution in r with normalisation - CALL loduni(2,p%R(1:p%Nploc)) - p%R(1:p%Nploc)=(plasmadim(3)+p%R(1:p%Nploc)*(plasmadim(4)-plasmadim(3)))/rnorm + CALL lodlinr(2,p%R(1:p%Nploc),plasmadim(3),plasmadim(4)) + p%R(1:p%Nploc)=p%R(1:p%Nploc)/rnorm ! Initial velocities distribution CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) END SUBROUTINE loaduniformRZ !_______________________________________________________________________________ SUBROUTINE loadDavidson(p, VR,VZ,VTHET, lodr) USE constants, ONLY: me, kb, elchar USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, width, qsim, msim, & & omegac, zgrid, nz, rnorm, n0, nblock, temp interface subroutine lodr(nbase,y,rminus,rplus) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rplus, rminus end subroutine end interface TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb, z REAL(kind=db) :: r0, deltar2, halfLz, Mirrorratio, Le, VOL INTEGER :: j, n, blockstart, blockend, addedpart, remainparts INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=p%Nploc Allocate(ra(nblock),rb(nblock), z(0:nblock)) !r0=(plasmadim(4)+plasmadim(3))/2 r0=sqrt(4*H0/(me*omegac**2)) halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) z(0)=plasmadim(1) DO n=1,nblock ! Compute limits in radius and load radii for each part Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1) z(n)=z(0)+n*(plasmadim(2)-plasmadim(1))/nblock deltar2=1-MirrorRatio*cos(2*pi*Le/width) rb(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2)) ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/H0*deltar2)) END DO VOL=SUM(2*pi*MINVAL(ra)*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) qsim=VOL*n0*elchar/nplasma msim=abs(qsim)/elchar*me p%weight=abs(qsim)/elchar p%m=me p%q=sign(elchar,qsim) p%qmRatio=p%q/p%m blockstart=1 blockend=0 ALLOCATE(blocksize(nblock)) WRITE(*,*) "blocksize: ", size(blocksize), nblock DO n=1,nblock blocksize(n)=nplasma/VOL*2*pi*MINVAL(ra)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock END DO remainparts=p%Nploc-SUM(blocksize) addedpart=1 n=nblock/2 j=1 DO WHILE(remainparts .GT. 0) blocksize(n)=blocksize(n)+addedpart remainparts=remainparts-addedpart n=n+j j=-1*(j+SIGN(1,j)) END DO CALL loadPartSlices(p, lodr, ra, rb, z, blocksize) IF(distribtype .eq. 5) THEN CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) VZ=VZ/4 VR=VR*8 VTHET=VTHET*8 ELSE Call loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) END IF END SUBROUTINE loadDavidson SUBROUTINE loadDavidsonVelocities(p, VR,VZ,VTHET, H0, P0) USE constants, ONLY: me, kb, elchar USE basic, ONLY: rnorm, Rcurv, B0, width, vnorm, zgrid, nz TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), INTENT(IN):: H0, P0 REAL(kind=db) :: athetpos, rg, zg, halfLz, Mirrorratio, Pcomp, Acomp INTEGER :: i MirrorRatio=(Rcurv-1)/(Rcurv+1) halfLz=(zgrid(nz)+zgrid(0))/2 ! Load velocities theta velocity ! Loading of r and z velocity is done in adapt_vinit to have ! access to parts%pot DO i=1,p%Nploc ! Interpolation for Magnetic potential rg=p%R(i)*rnorm zg=(p%Z(i)-halfLz)*rnorm Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width)) Pcomp=P0/rg/p%m Acomp=-p%qmRatio*Athetpos VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/p%m)),Pcomp+Acomp) !VTHET(i)=Pcomp+Acomp END DO VTHET=VTHET/vnorm VZ=0._db VR=0._db p%Davidson=.true. p%H0=H0 p%P0=P0 END SUBROUTINE loadDavidsonvelocities SUBROUTINE loadGaussianVelocities(p, VR,VZ,VTHET, temperature) USE basic, ONLY: vnorm USE constants, ONLY: kb REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN):: temperature REAL(kind=db):: vth ! Initial velocities distribution vth=sqrt(kb*temperature/p%m)/vnorm !thermal velocity CALL lodgaus(3,VZ) CALL lodgaus(5,VR) CALL lodgaus(7,VTHET) VZ=VZ*vth VR=VR*vth VTHET=VTHET*vth p%temperature=temperature p%Davidson=.false. END SUBROUTINE loadGaussianVelocities SUBROUTINE loadPartslices(p, lodr, ra, rb, z, blocksize) USE basic, ONLY: rnorm interface subroutine lodr(nbase,y,rminus,rplus) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rplus, rminus end subroutine end interface TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN)::ra(:), rb(:), z(0:) INTEGER, DIMENSION(:), INTENT(IN) :: blocksize INTEGER :: n, blockstart, blockend, nblock nblock=size(blocksize,1) blockstart=1 blockend=0 DO n=1,nblock blockstart=blockend+1 blockend=MIN(blockstart+blocksize(n)-1,p%Nploc) ! Initial distribution in z with normalisation between magnetic mirrors CALL loduni(1,p%Z(blockstart:blockend)) p%Z(blockstart:blockend)= (z(n-1)+p%Z(blockstart:blockend)*(z(n)-z(n-1)))/rnorm CALL lodr(2,p%R(blockstart:blockend),ra(n), rb(n)) p%R(blockstart:blockend)=p%R(blockstart:blockend)/rnorm END DO END SUBROUTINE loadPartslices SUBROUTINE loadPartFile(p, partfileindex, VR, VZ, VTHET) USE basic, ONLY: partfile, nplasma, lu_partfile implicit None TYPE(particles), INTENT(INOUT):: p REAL(kind=db), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VR, VZ, VTHET INTEGER, INTENT(IN):: partfileindex INTEGER:: nblock = 0 REAL(kind=db), Dimension(:), ALLOCATABLE:: ra, rb, z INTEGER, Dimension(:), ALLOCATABLE:: npartsslice INTEGER:: velocitytype=1 !< 1) gaussian with temp 2) Davidson with H0, P0 INTEGER:: radialtype=1 !< 1) 1/R 2) uniform 3) 1/R^2 4) gauss INTEGER:: npartsalloc !< initial size of particles arrays REAL(kind=db):: mass=me REAL(kind=db):: charge=-elchar REAL(kind=db):: weight=1.0 CHARACTER(len=256) :: header=' ' !< header of csv file REAL(kind=db):: H0=3.2e-14 !< Total energy REAL(kind=db):: P0=8.66e-25 !< Canonical angular momentum REAL(kind=db):: temperature=10000 !< temperature in kelvins real(kind=db):: n0 !< density factor LOGICAL :: is_test = .false. !< Defines if particle are test particles or tracers INTEGER:: i, ierr, openerr NAMELIST /partsload/ nblock, mass, charge, weight, npartsalloc, velocitytype, & & radialtype, temperature, H0, P0, is_test, n0 OPEN(UNIT=lu_partfile,FILE=trim(partfile(partfileindex)),ACTION='READ',IOSTAT=openerr) header=' ' IF(openerr .ne. 0) THEN CLOSE(unit=lu_partfile) RETURN END IF READ(lu_partfile,partsload) IF(mpirank .eq.0) THEN WRITE(*,'(a,i2,a,a)')"reading partfile: ", partfileindex, " ", partfile(partfileindex) WRITE(*,partsload) END IF IF( nblock .ge. 1) THEN ALLOCATE(z(0:nblock),ra(nblock),rb(nblock), npartsslice(nblock)) DO WHILE(header(1:8) .ne. '//slices') READ(lu_partfile,'(a)') header END DO DO i=1,nblock READ(lu_partfile,*) z(i-1),ra(i),rb(i),npartsslice(i) END DO READ(lu_partfile,*) z(nblock) CALL creat_parts(p,max(npartsalloc,sum(npartsslice))) p%Nploc=sum(npartsslice) p%Nptot=p%Nploc IF(partfileindex .eq. 1) nplasma=p%Nploc IF( allocated(VR) ) THEN DEALLOCATE(VR,VZ,VTHET) end if if(.not. allocated(VR)) THEN ALLOCATE(VR(p%Nploc)) ALLOCATE(VZ(p%Nploc)) ALLOCATE(VTHET(p%Nploc)) END IF p%m=mass p%q=charge p%weight=weight p%qmRatio=charge/mass p%is_test=is_test p%Newindex=sum(npartsslice) SELECT CASE(radialtype) CASE(1) ! 1/R distribution in R CALL loadPartslices(p, lodunir, ra, rb, z, npartsslice) CASE(2) ! flat top distribution in R CALL loadPartslices(p, lodlinr, ra, rb, z, npartsslice) CASE(3) ! 1/R^2 distribution in R CALL loadPartslices(p, lodinvr, ra, rb, z, npartsslice) CASE(4) ! gaussian distribution in R CALL loadPartslices(p, lodgausr, ra, rb, z, npartsslice) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT SELECT CASE(velocitytype) CASE(1) ! Gaussian with temperature CALL loadGaussianVelocities(p, VR, VZ, VTHET, temperature) CASE(2) CALL loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of velocity distribution:", velocitytype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT END IF CLOSE(unit=lu_partfile) END SUBROUTINE !=============================================================================== 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, dr, dz INTEGER, INTENT(IN) ::samplefactor TYPE(particles), INTENT(INOUT):: p INTEGER:: i, j, currentindex REAL(kind=db), DIMENSION(p%Nploc) :: spreaddir ! random direction for the spread of each initial macro particle REAL(kind=db) :: dir ! Direction in which the particle is moved REAL(kind=db) :: dl ! Particle displacement used for ! Load and scale the direction angle for spreading the new particles CALL loduni(2, spreaddir) spreaddir=spreaddir*2*pi/samplefactor dl=min(dz,minval(dr))/100 DO i=1,p%Nploc DO j=1,samplefactor-1 currentindex=p%Nploc+(i-1)*(samplefactor-1)+j CALL move_part(p,i,currentindex) p%partindex(currentindex)=currentindex dir = spreaddir(i)+2*pi*j/samplefactor p%R(currentindex)=p%R(currentindex) + dl*cos(dir) p%Z(currentindex)=p%Z(currentindex) + dl*sin(dir) END DO p%partindex(i)=i p%R(i)=p%R(i) + dl*cos(spreaddir(i)) p%Z(i)=p%Z(i) + dl*sin(spreaddir(i)) END DO nplasma=nplasma*samplefactor p%weight=p%weight/samplefactor p%Nploc=p%Nploc*samplefactor END SUBROUTINE upsample END MODULE beam diff --git a/src/stepon.f90 b/src/stepon.f90 index a624d57..d718b84 100644 --- a/src/stepon.f90 +++ b/src/stepon.f90 @@ -1,62 +1,64 @@ SUBROUTINE stepon ! ! Advance one time step ! USE basic USE constants USE fields USE beam USE maxwsrce USE celldiag INTEGER:: i DO i=1,nbspecies ! Boundary conditions for plasma particles outside the plasma region CALL bound(partslist(i)) ! Localisation of particles in cells (calculation of the r and z indices) CALL localisation(partslist(i)) END DO ! Cell diag quantities IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN CALL celldiag_save(time, fidres) END IF ! The particles are injected by the source CALL maxwsrce_inject(time) DO i=1,nbspecies ! Assemble right hand side of Poisson equation CALL rhscon(partslist(i)) END DO + if (.not. nlfreezephi) THEN ! Solve Poisson equation - CALL poisson(splrz) + CALL poisson(splrz) + end if DO i=1,nbspecies ! Compute the electric field at the particle position CALL EFieldscompatparts(partslist(i)) ! Compute the magnetic field at the particle position call comp_mag_p(partslist(i)) ! Solve Newton eq. and advance velocity by delta t CALL comp_velocity(partslist(i)) ! Calculate new positions of particles at time t+delta t CALL push(partslist(i)) ! Compute the energy of added particles CALL calc_newparts_energy(partslist(i)) END DO ! Calculate main physical quantities CALL diagnostics END SUBROUTINE stepon