diff --git a/src/.depend b/src/.depend index 956458c..87d0d46 100644 --- a/src/.depend +++ b/src/.depend @@ -1,33 +1,20 @@ auxval.o : auxval.f90 fields_mod.o basic_mod.o constants.o -auxval__genmod.o : auxval__genmod.f90 basic_mod.o : basic_mod.f90 mpihelper_mod.o constants.o beam_mod.o : beam_mod.f90 basic_mod.o mpihelper_mod.o constants.o chkrst.o : chkrst.f90 fields_mod.o beam_mod.o basic_mod.o -chkrst__genmod.o : chkrst__genmod.f90 constants.o : constants.f90 -cp2bk__genmod.o : cp2bk__genmod.f90 diagnose.o : diagnose.f90 xg_mod.o beam_mod.o basic_mod.o -diagnose__genmod.o : diagnose__genmod.f90 endrun.o : endrun.f90 fields_mod.o beam_mod.o basic_mod.o -endrun__genmod.o : endrun__genmod.f90 energies.o : energies.f90 basic_mod.o fields_mod.o : fields_mod.f90 mpihelper_mod.o beam_mod.o basic_mod.o constants.o inital.o : inital.f90 mpihelper_mod.o fields_mod.o beam_mod.o basic_mod.o -inital__genmod.o : inital__genmod.f90 main.o : main.f90 basic_mod.o mpihelper_mod.o : mpihelper_mod.f90 constants.o mv2bk.o : mv2bk.f90 -mv2bk__genmod.o : mv2bk__genmod.f90 newrun.o : newrun.f90 -newrun__genmod.o : newrun__genmod.f90 restart.o : restart.f90 -restart__genmod.o : restart__genmod.f90 resume.o : resume.f90 fields_mod.o beam_mod.o -resume__genmod.o : resume__genmod.f90 start.o : start.f90 basic_mod.o -start__genmod.o : start__genmod.f90 stepon.o : stepon.f90 beam_mod.o fields_mod.o constants.o basic_mod.o -stepon__genmod.o : stepon__genmod.f90 tesend.o : tesend.f90 basic_mod.o -tesend__genmod.o : tesend__genmod.f90 xg_mod.o : xg_mod.f90 fields_mod.o beam_mod.o basic_mod.o constants.o diff --git a/src/Makefile b/src/Makefile index bef4f93..2b73c1b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,64 +1,78 @@ .DEFAULT_GOAL := all ifeq ($(PLATFORM),) $(error Please specify the env variable PLATFORM (mac, intel)) else $(info *** Using $(PLATFORM).mk ***) include $(PLATFORM).mk endif include .depend PROG = espic2d SRCS = main.f90 basic_mod.f90 newrun.f90 restart.f90 \ auxval.f90 inital.f90 resume.f90 start.f90 diagnose.f90 \ stepon.f90 tesend.f90 endrun.f90 chkrst.f90 mv2bk.f90 \ constants.f90 xg_mod.f90 fields_mod.f90 beam_mod.f90 \ mpihelper_mod.f90 SRCS_C = extra.c +MKDIR_P = mkdir -p + +OUT_DIR = release + F90FLAGS += -I$(BSPLINES)/include -I$(FUTILS)/include \ - -I$(MUMPS)/include + -I$(MUMPS)/include -I../src LDFLAGS += -L$(BSPLINES)/lib -L$(FUTILS)/lib -L${HDF5}/lib -L${HDF5}/lib \ -L$(MUMPS)/lib -L$(PARMETIS)/lib -L/usr/local/xgrafix_1.2/src-double LIBS += -lbsplines -lpppack -lfutils -lhdf5_fortran -lhdf5 -lz $(MUMPSLIBS) -lpputils2 \ -lXGF -lXGC -lX11 -OBJS = ${SRCS:.f90=.o} ${SRCS_F90:.F90=.o} ${SRCS_C:.c=.o} +OBJS =${SRCS:.f90=.o} ${SRCS_F90:.F90=.o} ${SRCS_C:.c=.o} + +debug: F90FLAGS = $(DEBUGFLAGS) -I$(BSPLINES)/include -I$(FUTILS)/include \ + -I$(MUMPS)/include +debug: OUT_DIR=debug +debug: $(OUT_DIR) all + +profile: F90FLAGS+=$(PROFILEFLAGS) +profile: LDFLAGS+= $(PROFILEFLAGS) +profile: OUT_DIR=profile +profile: $(OUT_DIR) all -all: $(PROG) +.PHONY: directories + +all: directories $(PROG) $(PROG): $(OBJS) - $(F90) $(LDFLAGS) -o $@ $(OBJS) $(LIBS) + $(F90) $(LDFLAGS) -o $@ $(addprefix ./$(OUT_DIR)/,$(OBJS)) $(LIBS) + mv *.mod ./$(OUT_DIR)/ tags: etags *.f90 clean: rm -f $(OBJS) *.mod distclean: clean rm -f $(PROG) *~ a.out *.o TAGS -debug: F90FLAGS = $(DEBUGFLAGS) -I$(BSPLINES)/include -I$(FUTILS)/include \ - -I$(MUMPS)/include -debug: all +directories: ${OUT_DIR} -profile: F90FLAGS+=$(PROFILEFLAGS) -profile: LDFLAGS+= $(PROFILEFLAGS) -profile: all +${OUT_DIR}: + ${MKDIR_P} ${OUT_DIR} .SUFFIXES: $(SUFFIXES) .f90 .c .f90.o: - $(F90) $(F90FLAGS) -c $< + $(F90) $(F90FLAGS) -c -o ./$(OUT_DIR)/$@ $< .c.o: - $(CC) $(CCFLAGS) -c $< + $(CC) $(CCFLAGS) -c -o ./$(OUT_DIR)/$@ $< depend .depend: makedepf90 *.[fF]90 > .depend diff --git a/src/auxval.f90 b/src/auxval.f90 index e57dee7..6abb70a 100644 --- a/src/auxval.f90 +++ b/src/auxval.f90 @@ -1,100 +1,101 @@ 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=qsim/msim*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) ) + ALLOCATE( Zbounds(0:mpisize), Nplocs_all(0:mpisize-1)) Zbounds(0) = Zgrid(0) Zbounds(mpisize) = Zgrid(nz) leftproc=mpirank-1 IF(mpirank .eq. 0) leftproc = mpisize-1 rightproc = mpirank+1 IF(mpirank .eq. mpisize-1) rightproc = 0 WRITE(*,*) "mpirank, left, right", mpirank, 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 ALLOCATE(partdensity(nz*nr)) !________________________________________________________________________________ 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 f8945fa..afb8487 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,297 +1,303 @@ 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=.TRUE. ! Show graphical interface Xgrafix INTEGER :: nrun=1 ! Number of time steps to run DOUBLE PRECISION :: job_time=3600.0 ! Time allocated to this job in seconds DOUBLE PRECISION :: tmax=100000.0 ! Maximum simulation time DOUBLE PRECISION :: extra_time=60.0 ! Extra time allocated DOUBLE PRECISION :: dt=1 ! Time step DOUBLE PRECISION :: 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, it2d, itparts ! Number of iterations between 0d, 2d, particles, values writes to hdf5 INTEGER :: ittext=1 ! Number of iterations between text outputs in the console 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, leftproc ! Rank of next and 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 ! FID for resfile INTEGER :: fidrst ! FID 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 classicaly LOGICAL :: nlperiod(2)=(/.false.,.false./) ! Set periodic splines on or off INTEGER :: nplasma ! Number of superparticles INTEGER :: nblock ! Number of intervals in Z for stable distribution initialisation DOUBLE PRECISION :: potinn,potout ! Potential at inner and outer metalic walls DOUBLE PRECISION :: B0 ! Max magnitude of magnetic field DOUBLE PRECISION, allocatable :: Bz(:), Br(:) ! Magnetic field components DOUBLE PRECISION, allocatable :: Athet(:) ! Magnetic potential TYPE(spline2d) :: splrz ! Spline at r and z DOUBLE PRECISION, allocatable :: Ez(:), Er(:) ! Electric field components DOUBLE PRECISION, allocatable :: pot(:) ! Electro static potential DOUBLE PRECISION :: radii(4) ! Inner and outer radius of cylinder and radii of fine mesh region DOUBLE PRECISION :: 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 DOUBLE PRECISION :: H0=0, P0=0 ! Initial value of Hamiltonian and canonical angular momentum ! for distribution 2 DOUBLE PRECISION :: lz(2) ! Length of cylinder in z direction DOUBLE PRECISION :: n0 ! Average density of physical plasma DOUBLE PRECISION, ALLOCATABLE:: partdensity(:) ! Number of particles per gridcell computed every it2d INTEGER :: nz, nr, nnr(3) ! Number of intervals in z and r DOUBLE PRECISION :: dz, dr(3) ! Cell size in z and r for each region DOUBLE PRECISION, allocatable :: zgrid(:), rgrid(:) ! Nodes positions DOUBLE PRECISION :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm ! Normalization constants DOUBLE PRECISION :: qsim ! Charge of superparticles DOUBLE PRECISION :: 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 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhs ! right hand side of the poisson equation solver DOUBLE PRECISION :: omegac ! Cyclotronic frequency DOUBLE PRECISION :: omegap ! Plasma frequency DOUBLE PRECISION :: V ! Normalised volume DOUBLE PRECISION :: temp ! Initial temperature of plasma DOUBLE PRECISION :: Rcurv ! Magnetic field curvature coefficient DOUBLE PRECISION :: Width ! Distance between two magnetic mirrors DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Zbounds ! Index of bounds for local processus in Z direction TYPE(BASICDATA) :: bdata CONTAINS ! !================================================================================ SUBROUTINE basic_data ! ! Define basic data ! use mpihelper IMPLICIT NONE ! ! Local vars and arrays - CHARACTER(len=256) :: line + 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, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0 !________________________________________________________________________________ ! 1. Process Standard Input File ! IF(mpirank .eq. 0) THEN - DO - READ(*,'(a)', END=110) line - WRITE(lu_in, '(a)') TRIM(line) - END DO -110 CONTINUE - REWIND(lu_in) + 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)) + !DO + ! READ(*,'(a)', END=110) line + ! WRITE(lu_in, '(a)') TRIM(line) + !END DO +!110 CONTINUE + !REWIND(lu_in) !________________________________________________________________________________ ! 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) bdata=BASICDATA( nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical, & & nplasma, nz, it0d, it2d, itparts, itgraph, distribtype, nblock, & & nrun, job_time, extra_time, tmax, dt, potinn, potout, B0, n0, temp, & & Rcurv, width, H0, P0, femorder, ngauss, nnr, lz, radii, plasmadim, resfile ) - END IF CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation CALL MPI_Bcast(bdata, 1, basicdata_type, 0, MPI_COMM_WORLD, ierr) CALL readbdata ! 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 DOUBLE PRECISION, OPTIONAL, INTENT(out) :: eltime ! INTEGER, PARAMETER :: ncmax=128 INTEGER, SAVE :: icall=0, nc=0 DOUBLE PRECISION, SAVE :: startt0=0.0 CHARACTER(len=16), SAVE :: which(ncmax) INTEGER :: lstr, found, i DOUBLE PRECISION :: seconds DOUBLE PRECISION, 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 END SUBROUTINE readbdata !================================================================================ END MODULE basic diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index db48056..5da6530 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,1148 +1,1136 @@ MODULE beam ! USE constants USE mpi USE mpihelper USE basic, ONLY: mpirank, mpisize IMPLICIT NONE TYPE particles INTEGER :: Nptot INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex, Zindex ! Indices of the particles on the grid DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: & & R,Z, & & WZ, WR, & & pot, Er, Ez DOUBLE PRECISION, DIMENSION(:),POINTER:: & & UR, URold, & & UTHET, UTHETold, & & UZ, UZold, & & Gamma, Gammaold END TYPE particles ! TYPE(particles) :: parts ! Storage for all the particles SAVE :: parts TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffer to send and receive particle from left and right processes ! Diagnostics (scalars) DOUBLE PRECISION :: ekin ! Total kinetic energz (J) DOUBLE PRECISION :: epot ! Total potential energy (J) DOUBLE PRECISION :: etot, etot0 ! Current and Initial total energy (J) INTEGER, DIMENSION(3), SAVE :: ireducerequest=MPI_REQUEST_NULL ! INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !F.B. get local numbers of particles ! CONTAINS 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%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%Zindex(nparts)) ALLOCATE(p%pot(nparts)) ALLOCATE(p%Er(nparts)) ALLOCATE(p%Ez(nparts)) ALLOCATE(p%GAMMAold(nparts)) parts%URold=0 parts%UZold=0 parts%UTHETold=0 parts%rindex=0 parts%zindex=0 parts%WR=0 parts%WZ=0 parts%UR=0 parts%UZ=0 parts%UTHET=0 parts%Z=0 parts%R=0 parts%Gamma=1 parts%Er=0 parts%Ez=0 parts%pot=0 parts%gammaold=1 END IF END SUBROUTINE creat_parts !______________________________________________________________________________ SUBROUTINE load_parts ! Loads the particles at the beginning of the simulation and create the parts variable if necessary - - USE constants, ONLY: me, kb, elchar - USE basic, ONLY: nplasma, rnorm, plasmadim, temp, distribtype, H0, P0, Rcurv, ierr, B0, width, qsim, msim, vnorm, & - & omegac, zgrid, nz, rnorm, tnorm, bnorm, V, n0, nblock, nlclassical + USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical USE mpi - DOUBLE PRECISION :: vth, v2 - INTEGER :: i, j, n, blockstart, blockend, sumblocksize, addedpart, remainparts - INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize - DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET, ra, rb - DOUBLE PRECISION :: r0, athetpos, Vperp, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL + + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) - Allocate(ra(nblock),rb(nblock)) ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V and uniform in R - ! Initial distribution in z with normalisation - CALL loduni(1,parts%Z(:)) - parts%Z(:)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(:))/rnorm - ! Initial distribution in r with normalisation - CALL loduni(2,parts%R(:)) - parts%R=(plasmadim(3)+parts%R(:)*(plasmadim(4)-plasmadim(3)))/rnorm - - ! Initial velocities distribution - vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity - CALL lodgaus(3,VZ(:)) - CALL lodgaus(5,VR(:)) - CALL lodgaus(7,VTHET(:)) - VZ=VZ*vth - VR=VR*vth - VTHET=VTHET*vth - + CALL loaduniformRZ(VR,VZ,VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 - r0=(plasmadim(4)+plasmadim(3))/2 - halfLz=(zgrid(nz)+zgrid(0))/2 - - 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) - MirrorRatio=(Rcurv-1)/(Rcurv+1) - deltar2=MirrorRatio*cos(2*pi*Le/width) - rb(n)=r0*(1+sqrt(deltar2))/(1-deltar2) - ra(n)=r0*(1-sqrt(deltar2))/(1-deltar2) - END DO - VOL=SUM(2*pi*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) - qsim=VOL*n0*elchar/nplasma - msim=abs(qsim)/elchar*me - V=VOL/rnorm**3 - - !H0=me*omegac**2*r0**2*0.5 - gamma=1/sqrt(1-omegac**2*r0**2/vlight**2) - !P0=H0/omegac - - blockstart=1 - blockend=0 - ALLOCATE(blocksize(nblock)) - DO n=1,nblock - blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock - END DO - remainparts=parts%Nptot-SUM(blocksize) - addedpart=1 - n=nblock/2 - j=1 - DO WHILE(remainparts .GT. 0) - blocksize(n)=blocksize(n)+addedpart - remainparts=remainparts-addedpart - n=n+j - j=-1*(j+SIGN(1,j)) - END DO - DO n=1,nblock - blockstart=blockend+1 - blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot) - ! Initial distribution in z with normalisation between magnetic mirrors - CALL loduni(1,parts%Z(blockstart:blockend)) - parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm - !CALL lodinvr(3,parts%R(blockstart:blockend),ra(n),rb(n)) - CALL loduni(3,parts%R(blockstart:blockend)) - parts%R(blockstart:blockend)=ra(n)+(rb(n)-ra(n))*parts%R(blockstart:blockend) - parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm - END DO - - - - ! 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=abs(P0/rg/me) - Acomp=-SIGN(elchar/me*Athetpos,qsim) - VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp) - - - END DO - VTHET=VTHET/vnorm - VZ=0._db - VR=0._db - - + CALL loadDavidson(VR,VZ,VTHET) + CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2 + CALL loadDavidson(VR,VZ,VTHET) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT IF(nlclassical) THEN parts%Gamma=1 ELSE parts%Gamma=sqrt(1/(1-VR**2+VZ**2+VTHET**2)) END IF parts%UR=parts%Gamma*VR parts%UZ=parts%Gamma*VZ parts%UTHET=parts%Gamma*VTHET DEALLOCATE(VZ, VR, VTHET) CALL localisation IF(mpisize .gt. 1) THEN - ALLOCATE(Nplocs_all(0:mpisize-1)) CALL distribpartsonprocs END IF END SUBROUTINE load_parts !________________________________________________________________________________ SUBROUTINE bound USE basic, ONLY: zgrid, nz, Zbounds, cstep, mpirank INTEGER :: i, rsendnbparts=0, lsendnbparts=0 INTEGER, DIMENSION(parts%Nptot) :: sendhole sendhole=0 lsendnbparts=0 rsendnbparts=0 IF (parts%Nptot .NE. 0) THEN ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot IF(cstep .ne. 0) THEN IF (parts%Z(i) .ge. Zbounds(mpirank+1)) THEN !$OMP CRITICAL (nbparts) rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i !$OMP END CRITICAL (nbparts) ELSE IF (parts%Z(i) .lt. Zbounds(mpirank)) THEN !$OMP CRITICAL (nbparts) lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i !$OMP END CRITICAL (nbparts) END IF END IF DO WHILE (parts%Z(i) .GT. zgrid(nz)) parts%Z(i) = parts%Z(i) - zgrid(nz) + zgrid(0) END DO DO WHILE (parts%Z(i) .LT. zgrid(0)) parts%Z(i) = parts%Z(i) + zgrid(nz) - zgrid(0) END DO END DO !$OMP END PARALLEL DO CALL particlescommunication(lsendnbparts, rsendnbparts, sendhole) END IF END subroutine bound !________________________________________________________________________________ SUBROUTINE localisation - USE basic, ONLY: zgrid, rgrid, step, nz, dz, nr, nnr, dr, rnorm + USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm INTEGER :: j,i, nblostparts=0 INTEGER, DIMENSION(parts%Nptot) :: losthole i=1 losthole=0 IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot parts%zindex(i)=(parts%Z(i)-zgrid(0))/dz IF (parts%R(i) .GT. rgrid(0) .AND. parts%R(i) .LT. rgrid(nnr(1))) THEN parts%rindex(i)=(parts%R(i)-rgrid(0))/dr(1) ELSE IF(parts%R(i) .GT. rgrid(nnr(1)) .AND. parts%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)))/dr(2)+nnr(1) ELSE IF(parts%R(i) .GT. rgrid(nnr(1)+nnr(2)) .AND. parts%R(i) .LT. rgrid(nr)) THEN parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)+nnr(2)))/dr(3)+nnr(1)+nnr(2) ELSE WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in r:", parts%R(i)*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm !nlend=.TRUE. WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Particle removed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" !$OMP CRITICAL (lostparts) nblostparts=nblostparts+1 losthole(nblostparts)=i !$OMP END CRITICAL (lostparts) END IF END DO !$OMP END PARALLEL DO IF(nblostparts.gt.0) THEN DO j=1,nblostparts CALL move_part(parts%Nptot, losthole(j)) parts%Nptot=parts%Nptot-1 END DO END IF !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot parts%WZ(i)=(parts%Z(i)-zgrid(parts%zindex(i)))/(zgrid(parts%zindex(i)+1)-zgrid(parts%zindex(i))); parts%WR(i)=(parts%R(i)-rgrid(parts%rindex(i)))/(rgrid(parts%rindex(i)+1)-rgrid(parts%rindex(i))); END DO !$OMP END PARALLEL DO END IF END SUBROUTINE localisation !________________________________________________________________________________________________________ SUBROUTINE comp_velocity ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : nlclassical ! Store old Velocities CALL swappointer(parts%UZold, parts%UZ) CALL swappointer(parts%URold, parts%UR) CALL swappointer(parts%UTHETold, parts%UTHET) CALL swappointer(parts%Gammaold, parts%Gamma) IF (nlclassical) THEN CALL comp_velocity_fun(gamma_classical) ELSE CALL comp_velocity_fun(gamma_relativistic) END IF END SUBROUTINE comp_velocity !____________________________________________ SUBROUTINE comp_velocity_fun(gamma) ! ! Computes the new velocity of the particles due to Lorentz force ! - USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz, nlclassical, nplasma + USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz interface subroutine gamma(gam, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gam end subroutine end interface DOUBLE PRECISION :: tau DOUBLE PRECISION:: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2 INTEGER:: J1, J2, J3, J4 INTEGER:: i ! Normalized time increment tau=omegac/2/omegap*dt/tnorm IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,parts%Nptot J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1 J2=J1+1 J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1 J4=J3+1 ! Interpolation for magnetic field BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) & & +parts%WZ(i)*parts%WR(i)*Bz(J1) BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Br(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Br(J2) & & +parts%WZ(i)*parts%WR(i)*Br(J1) ! First half of electric pulse parts%UZ(i)=parts%UZold(i)+parts%Ez(i)*tau parts%UR(i)=parts%URold(i)+parts%ER(i)*tau CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) ! Rotation along magnetic field ZBZ=tau*BRZ/parts%Gamma(i) ZBR=tau*BRR/parts%Gamma(i) ZPZ=parts%UZ(i)-ZBR*parts%UTHETold(i) !u'_{z} ZPR=parts%UR(i)+ZBZ*parts%UTHETold(i) !u'_{r} ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZ(i)-ZBZ*parts%UR(i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR parts%UZ(i)=parts%UZ(i)-ZBR2*ZPTHET !u+_{z} parts%UR(i)=parts%UR(i)+ZBZ2*ZPTHET !u+_{r} parts%UTHET(i)=parts%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration parts%UZ(i)=parts%UZ(i)+parts%EZ(i)*tau parts%UR(i)=parts%UR(i)+parts%ER(i)*tau CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) END DO !$OMP END PARALLEL DO END IF END SUBROUTINE comp_velocity_fun !______________________________________________________________________________ SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gamma gamma=1.0 END SUBROUTINE gamma_classical !______________________________________________________________________________ SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gamma gamma=sqrt(1+UZ**2+UR**2+UTHET**2) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE push - Use basic, ONLY: dt, tnorm, nplasma + Use basic, ONLY: dt, tnorm DOUBLE PRECISION:: XP, YP, COSA, SINA, U1, U2 INTEGER :: i IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2) DO i=1,parts%Nptot ! Local Cartesian coordinates XP=parts%R(i)+dt/tnorm*parts%UR(i)/parts%Gamma(i) YP=dt/tnorm*parts%UTHET(i)/parts%Gamma(i) ! Conversion to cylindrical coordiantes parts%Z(i)=parts%Z(i)+dt/tnorm*parts%UZ(i)/parts%Gamma(i) parts%R(i)=sqrt(XP**2+YP**2) ! Velocity in rotated reference frame IF (parts%R(i) .EQ. 0) THEN COSA=1 SINA=0 ELSE COSA=XP/parts%R(i) SINA=YP/parts%R(i) END IF U1=COSA*parts%UR(i)+SINA*parts%UTHET(i) U2=-SINA*parts%UR(i)+COSA*parts%UTHET(i) parts%UR(i)=U1 parts%UTHET(i)=U2 END DO !$OMP END PARALLEL DO END IF END SUBROUTINE push !_______________________________________________________________________________ !_______________________________________________________________________________ SUBROUTINE diagnostics ! ! Compute energies ! !!$ CALL energies USE constants, ONLY: vlight USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, partdensity, ierr, step, it2d, nlend, nr, itparts, nplasma DOUBLE PRECISION :: gammamid INTEGER :: i INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE) INTEGER, DIMENSION(:), ALLOCATABLE :: displs CALL MPI_WAIT(ireducerequest(3), stat(2*MPI_STATUS_SIZE+1:3*MPI_STATUS_SIZE), ierr) partdensity=0 ekin=0 epot=0 etot=0 !$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity) DEFAULT(SHARED) PRIVATE(i,gammamid) DO i=1,parts%Nptot ! epot=epot+qsim*parts%pot(i)*phinorm IF(.not. nlclassical) THEN gammamid=1/sqrt(1-(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))/(parts%Gammaold(i)*parts%Gamma(i))) ekin=ekin+msim*vlight**2*(gammamid-1) ELSE ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i))) END IF IF(modulo(step,it2d).eq. 0 .or. nlend) THEN partdensity(parts%zindex(i)+1+parts%rindex(i)*nz)=partdensity(parts%zindex(i)+1+parts%rindex(i)*nz)+1 END IF END DO !$OMP END PARALLEL DO IF(mpirank .eq.0 ) THEN CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(1), ierr) CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(2), ierr) ELSE CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(1), ierr) CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(2), ierr) END IF ! Stores the particle density on the grid IF(modulo(step,it2d).eq. 0 .or. nlend) THEN IF(mpirank .eq.0 ) THEN CALL MPI_IREDUCE(MPI_IN_PLACE, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) ELSE CALL MPI_IREDUCE(partdensity, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) END IF END IF IF(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. 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 CALL MPI_WAIT(ireducerequest(1), stat(1:MPI_STATUS_SIZE), ierr) CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr) etot=epot+ekin ! ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ.1) THEN etot0 = etot END IF IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend) .and. step .ne. 0) THEN IF(mpirank .ne. 0) THEN ALLOCATE(displs(0:mpisize-1)) displs=0 ! Send Particles informations to root process CALL MPI_Gatherv(parts%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) ELSE ALLOCATE(displs(0:mpisize-1)) - displs(0)=1 + displs(0)=0 Do i=1,mpisize-1 displs(i)=displs(i-1)+Nplocs_all(i-1) END DO ! Receive particle information from all processes CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) END IF END IF end subroutine diagnostics !_______________________________________________________________________________ SUBROUTINE adapt_vinit ! ! Computes the velocity at time -dt/2 from velocities computed at time 0 ! USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm - DOUBLE PRECISION :: tau, EPR, EPZ, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & + DOUBLE PRECISION :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, ZBZ2, ZBR2, Vperp, v2 INTEGER :: J1, J2, J3, J4, i DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET - IF(distribtype .EQ. 2) THEN + IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3) THEN ALLOCATE(VR(parts%Nptot),VZ(parts%Nptot),VTHET(parts%Nptot)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=parts%UTHET/parts%Gamma*vnorm DO i=1,parts%Nptot Vperp=sqrt(MAX(2*H0/me+2*elchar/me*parts%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN parts%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 parts%Gamma(i)=sqrt(1/(1-v2)) END IF parts%UR(i)=parts%Gamma(i)*VR(i)/vnorm parts%UZ(i)=parts%Gamma(i)*VZ(i)/vnorm parts%UTHET(i)=parts%Gamma(i)*VTHET(i)/vnorm END DO DEALLOCATE(VR,VZ,VTHET) END IF ! Normalised time increment tau=-omegac/2/omegap*dt/tnorm ! Store old Velocities CALL swappointer(parts%UZold, parts%UZ) CALL swappointer(parts%URold, parts%UR) CALL swappointer(parts%UTHETold, parts%UTHET) CALL swappointer(parts%Gammaold, parts%Gamma) IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,parts%Nptot J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1 J2=J1+1 J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1 J4=J3+1 - - EPZ=parts%Ez(i) - EPR=parts%Er(i) ! Interpolation for magnetic field BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) & & +parts%WZ(i)*parts%WR(i)*Bz(J1) BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Br(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Br(J2) & & +parts%WZ(i)*parts%WR(i)*Br(J1) ! Half inverse Rotation along magnetic field ZBZ=tau*BRZ/parts%Gammaold(i) ZBR=tau*BRR/parts%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(parts%UZold(i)-ZBR*parts%UTHETold(i))/SQR !u-_{z} ZPR=(parts%URold(i)+ZBZ*parts%UTHETold(i))/SQR !u-_{r} ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZold(i)-ZBZ*parts%URold(i))/SQR !u-_{theta} parts%UZ(i)=ZPZ parts%UR(i)=ZPR parts%UTHET(i)=ZPTHET ! half of decceleration - parts%UZ(i)=parts%UZ(i)+EPZ*tau - parts%UR(i)=parts%UR(i)+EPR*tau + parts%UZ(i)=parts%UZ(i)+parts%Ez(i)*tau + parts%UR(i)=parts%UR(i)+parts%Er(i)*tau IF(.not. nlclassical) THEN parts%Gamma(i)=sqrt(1+parts%UZ(i)**2+parts%UR(i)**2+parts%UTHET(i)**2) END IF END DO !$OMP END PARALLEL DO END IF END SUBROUTINE adapt_vinit !_______________________________________________________________________________ SUBROUTINE distribpartsonprocs ! Computes the start and end indices for the Z boundaries on local processus ! Computes the particle indices from initial particle loading vector, that stay in current process USE basic, ONLY: nz, nplasma, ierr, Zbounds - INTEGER:: k, sig + INTEGER:: k DOUBLE PRECISION:: idealnbpartsperproc INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k) INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector INTEGER:: partindexlast ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector INTEGER, DIMENSION(0:nz):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1) INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction INTEGER:: Zperproc - DOUBLE PRECISION:: allowedratio=0.2 partspercol=0 idealnbpartsperproc = FLOOR(REAL(nplasma)/REAL(mpisize)) Zmin=MINVAL(parts%Zindex) Zmax=MAXVAL(parts%Zindex) Zperproc=(Zmax-Zmin)/mpisize partindexstart=1 totparts=0 partindexlast=parts%Nptot DO k=0,mpisize-2 Nplocs_all(k)=idealnbpartsperproc END DO NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize) DO k=1, mpisize-1 Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc)) END DO ! Store local end of particle vector partindexlast=SUM(Nplocs_all(0:mpirank)) ! Store local start of particle vector IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1 IF(mpirank .ne. 0) THEN DO k= partindexstart, partindexlast CALL move_part(k,k-partindexstart+1) END DO END IF parts%Nptot=Nplocs_all(mpirank) WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot IF(INT(SUM(Nplocs_all)) .ne. nplasma) THEN WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " nplasma:", nplasma CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF END SUBROUTINE distribpartsonprocs !_______________________________________________________________________________ SUBROUTINE particlescommunication(lsendnbparts, rsendnbparts, sendholes) USE mpihelper, ONLY: particle_type USE basic, ONLY: rightproc, leftproc, ierr INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0 INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE) INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL 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 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 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 IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN CALL AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes) END IF CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr) CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr) 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 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 CALL Addincomingparts(rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes) 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 !_______________________________________________________________________________ SUBROUTINE Addincomingparts(rrecvnbparts, lrecvnbparts, sendnbparts, sendholes) ! USE mpihelper INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) INTEGER k,partpos, partdiff ! computes if we received less particles than we sent partdiff=sendnbparts-rrecvnbparts-lrecvnbparts IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN partpos=abs(sendholes(k)) ELSE parts%Nptot=parts%Nptot+1 partpos=parts%Nptot END IF CALL Insertincomingpart(rrecvpartbuff, k, partpos) END DO END IF IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN partpos=abs(sendholes(k+rrecvnbparts)) ELSE parts%Nptot=parts%Nptot+1 partpos=parts%Nptot END IF CALL Insertincomingpart(lrecvpartbuff, k, partpos) END DO END IF IF(partdiff .gt. 0) THEN DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts CALL move_part(parts%Nptot,abs(sendholes(k))) parts%Nptot = parts%Nptot-1 END DO END IF ! END SUBROUTINE Addincomingparts !_______________________________________________________________________________ SUBROUTINE Insertincomingpart(buffer, bufferindex, partsindex) USE mpihelper INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer parts%Rindex(partsindex) = buffer(bufferindex)%Rindex parts%Zindex(partsindex) = buffer(bufferindex)%Zindex parts%R(partsindex) = buffer(bufferindex)%R parts%Z(partsindex) = buffer(bufferindex)%Z parts%UZ(partsindex) = buffer(bufferindex)%UZ parts%UR(partsindex) = buffer(bufferindex)%UR parts%UTHET(partsindex) = buffer(bufferindex)%UTHET parts%Gamma(partsindex) = buffer(bufferindex)%Gamma ! ! END SUBROUTINE Insertincomingpart !_______________________________________________________________________________ SUBROUTINE Insertsentpart(buffer, bufferindex, partsindex) USE mpihelper INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer buffer(bufferindex)%Rindex = parts%Rindex(partsindex) buffer(bufferindex)%Zindex = parts%Zindex(partsindex) buffer(bufferindex)%R = parts%R(partsindex) buffer(bufferindex)%Z = parts%Z(partsindex) buffer(bufferindex)%UZ = parts%UZ(partsindex) buffer(bufferindex)%UR = parts%UR(partsindex) buffer(bufferindex)%UTHET = parts%UTHET(partsindex) buffer(bufferindex)%Gamma = parts%Gamma(partsindex) ! ! END SUBROUTINE Insertsentpart !_______________________________________________________________________________ SUBROUTINE AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes) ! USE mpihelper INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff) IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff) 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))) Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(sendholes(k)) IF(sendholes(k) .GT. 0) THEN rsendpos=rsendpos+1 CALL Insertsentpart(rsendpartbuff, rsendpos, partpos) ELSE IF(sendholes(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(lsendpartbuff, lsendpos, partpos) END IF END DO ! ! END SUBROUTINE AddPartSendBuffers !_______________________________________________________________________________ SUBROUTINE exchange_parts(index1, index2) INTEGER, INTENT(IN) :: index1, index2 - DOUBLE PRECISION:: R,Z,UR,URold, UZ, UZold, UTHET, UTHETold, Gamma, Gammaold, WZ, WR, Er, Ez, pot + DOUBLE PRECISION:: R,Z,UR, UZ, UTHET, Gamma INTEGER :: Rindex, Zindex !! Exchange particle at index1 with particle at index2 ! Store part at index1 in temporary value Gamma = parts%Gamma(index1) - Gammaold = parts%Gammaold(index1) R = parts%R(index1) Z = parts%Z(index1) UR = parts%UR(index1) - URold = parts%URold(index1) UTHET = parts%UTHET(index1) - UTHETold = parts%UTHETold(index1) UZ = parts%UZ(index1) - UZold = parts%UZold(index1) - WZ = parts%WZ(index1) - WR = parts%WR(index1) - Er = parts%Er(index1) - Ez = parts%Ez(index1) - pot = parts%pot(index1) Rindex = parts%Rindex(index1) Zindex = parts%Zindex(index1) ! Move part at index2 in part at index 1 parts%Gamma(index1) = parts%Gamma(index2) - parts%Gammaold(index1) = parts%Gammaold(index2) parts%R(index1) = parts%R(index2) parts%Z(index1) = parts%Z(index2) parts%UR(index1) = parts%UR(index2) - parts%URold(index1) = parts%URold(index2) parts%UTHET(index1) = parts%UTHET(index2) - parts%UTHETold(index1) = parts%UTHETold(index2) parts%UZ(index1) = parts%UZ(index2) - parts%UZold(index1) = parts%UZold(index2) - parts%WZ(index1) = parts%WZ(index2) - parts%WR(index1) = parts%WR(index2) - parts%Er(index1) = parts%Er(index2) - parts%Ez(index1) = parts%Ez(index2) - parts%pot(index1) = parts%pot(index2) parts%Rindex(index1) = parts%Rindex(index2) parts%Zindex(index1) = parts%Zindex(index2) ! Move temporary values from part(index1) to part(index2) parts%Gamma(index2) = Gamma - parts%Gammaold(index2) = Gammaold parts%R(index2) = R parts%Z(index2) = Z parts%UR(index2) = UR - parts%URold(index2) = URold parts%UTHET(index2) = UTHET - parts%UTHETold(index2) = UTHETold parts%UZ(index2) = UZ - parts%UZold(index2) = UZold - parts%WZ(index2) = WZ - parts%WR(index2) = WR - parts%Er(index2) = Er - parts%Ez(index2) = Ez - parts%pot(index2) = pot parts%Rindex(index2) = Rindex parts%Zindex(index2) = Zindex END SUBROUTINE exchange_parts !_______________________________________________________________________________ SUBROUTINE move_part(sourceindex, destindex) !! Move particle with index sourceindex to particle with index destindex !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex ! Move part at sourceindex in part at destindex parts%Gamma(destindex) = parts%Gamma(sourceindex) parts%R(destindex) = parts%R(sourceindex) parts%Z(destindex) = parts%Z(sourceindex) parts%UR(destindex) = parts%UR(sourceindex) parts%UTHET(destindex) = parts%UTHET(sourceindex) parts%UZ(destindex) = parts%UZ(sourceindex) parts%Rindex(destindex) = parts%Rindex(sourceindex) parts%Zindex(destindex) = parts%Zindex(sourceindex) END SUBROUTINE move_part !_______________________________________________________________________________ !________________________________________________________________________________ ! SUBROUTINE loduni(nbase, y) ! ! Load a uniform distribution using the Hammersley's sequence. ! (nbase = 0 => Random sampling) ! INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER :: n, i ! n = SIZE(y) ! SELECT CASE (nbase) CASE(0) CALL RANDOM_NUMBER(y) CASE(1) DO i=1,n y(i) = (i-0.5_db)/n END DO CASE(2:) DO i=1,n y(i) = rev(nbase,i) END DO CASE default WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase END SELECT ! END SUBROUTINE loduni !________________________________________________________________________________ SUBROUTINE lodgaus(nbase, y) ! ! Sample y from the Gauss distributrion ! DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase ! INTEGER :: n, i, iflag DOUBLE PRECISION :: r(SIZE(y)) DOUBLE PRECISION :: t ! DOUBLE PRECISION :: c0, c1, c2 DOUBLE PRECISION :: d1, d2, d3 DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/ DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/ ! n = SIZE(y) CALL loduni(nbase, r) r = 1.0E-5_db + 0.99998_db*r ! DO i=1,n iflag = 1 IF (r(i) .GT. 0.5_db) THEN r(i) = 1.0_db - r(i) iflag = -1 END IF t = SQRT(LOG(1.0_db/(r(i)*r(i)))) y(i) = t - (c0+c1*t+c2*t**2) / (1.0_db+d1*t+d2*t**2+d3*t**3) y(i) = y(i) * iflag END DO y = y - SUM(y)/REAL(n,db) END SUBROUTINE lodgaus !________________________________________________________________________________ SUBROUTINE lodinvr(nbase, y, ra, rb) ! ! Sample y from the distribution f=1/(r)H(r-ra)H(rb-r) for where H is the heavyside function ! DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra ! DOUBLE PRECISION :: r(SIZE(y)) ! CALL loduni(nbase, r) r=r*log(rb/ra) y=ra*exp(r) END SUBROUTINE lodinvr !________________________________________________________________________________ REAL(db) FUNCTION rev(nbase,i) ! ! Return an element of the Hammersley's sequence ! INTEGER, INTENT(IN) :: nbase ! Base of the Hammersley's sequence (IN) INTEGER, INTENT(IN) :: i ! Index of the sequence (IN) ! ! Local vars INTEGER :: j1, j2 REAL(db) :: xs, xsi !----------------------------------------------------------------------- xs = 0._db xsi = 1.0_db j2 = i DO xsi = xsi/nbase j1 = j2/nbase xs = xs + (j2-nbase*j1)*xsi j2 = j1 IF( j2.LE.0 ) EXIT END DO rev = xs END FUNCTION rev !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) DOUBLE PRECISION :: bessi1,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 !________________________________________________________________________________ SUBROUTINE destroy_parts(p) TYPE(particles) :: p p%nptot=0 IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z) IF(ALLOCATED(p%R)) DEALLOCATE(p%R) IF(ALLOCATED(p%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) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE clean_beam ! CALL destroy_parts(parts) ! ! END SUBROUTINE clean_beam !________________________________________________________________________________ RECURSIVE SUBROUTINE quicksortparts(leftlimit, rightlimit) ! Sorts the particle according to their Z position using quicksort algorithm INTEGER,INTENT(IN):: leftlimit, rightlimit DOUBLE PRECISION:: pivot - INTEGER::i, cnt + INTEGER::i, cnt, mid IF(leftlimit .ge. rightlimit) RETURN + mid=(leftlimit+rightlimit)/2 + IF(parts%Z(mid).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,mid) + IF(parts%Z(rightlimit).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,rightlimit) + IF(parts%Z(mid).lt.parts%Z(rightlimit)) CALL exchange_parts(rightlimit,mid) + ! Store the pivot point for comparison pivot=parts%Z(rightlimit) cnt=leftlimit ! Move all parts with Z smaller than pivot to the left of pivot DO i=leftlimit, rightlimit IF(parts%Z(i) .le. pivot) THEN CALL exchange_parts(i,cnt) cnt=cnt+1 END IF END DO + WRITE(*,*)"sorted parts: ", leftlimit,"to", rightlimit CALL quicksortparts(leftlimit,cnt-2) CALL quicksortparts(cnt,rightlimit) END SUBROUTINE quicksortparts !________________________________________________________________________________ SUBROUTINE swappointer( pointer1, pointer2) DOUBLE PRECISION, DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2 DOUBLE PRECISION, DIMENSION(:), POINTER:: temppointer temppointer=>pointer1 pointer1=>pointer2 pointer2=>temppointer END SUBROUTINE swappointer !_______________________________________________________________________________ + SUBROUTINE loaduniformRZ(VR,VZ,VTHET) + USE basic, ONLY: plasmadim, rnorm, temp, vnorm + USE constants, ONLY: me, kb, elchar + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET + DOUBLE PRECISION:: vth + ! Initial distribution in z with normalisation + CALL loduni(1,parts%Z(:)) + parts%Z(:)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(:))/rnorm + ! Initial distribution in r with normalisation + CALL loduni(2,parts%R(:)) + parts%R=(plasmadim(3)+parts%R(:)*(plasmadim(4)-plasmadim(3)))/rnorm + + ! Initial velocities distribution + vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity + CALL lodgaus(3,VZ(:)) + CALL lodgaus(5,VR(:)) + CALL lodgaus(7,VTHET(:)) + VZ=VZ*vth + VR=VR*vth + VTHET=VTHET*vth + END SUBROUTINE loaduniformRZ +!_______________________________________________________________________________ + SUBROUTINE loadDavidson(VR,VZ,VTHET) + 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 + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE::ra, rb + DOUBLE PRECISION :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL + 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(2*H0/(me*omegac**2)) + halfLz=(zgrid(nz)+zgrid(0))/2 + + 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) + MirrorRatio=(Rcurv-1)/(Rcurv+1) + deltar2=MirrorRatio*cos(2*pi*Le/width) + rb(n)=r0*(1+sqrt(deltar2))/(1-deltar2) + ra(n)=r0*(1-sqrt(deltar2))/(1-deltar2) + END DO + VOL=SUM(2*pi*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) + qsim=VOL*n0*elchar/nplasma + msim=abs(qsim)/elchar*me + V=VOL/rnorm**3 + + !H0=me*omegac**2*r0**2*0.5 + gamma=1/sqrt(1-omegac**2*r0**2/vlight**2) + !P0=H0/omegac + + blockstart=1 + blockend=0 + ALLOCATE(blocksize(nblock)) + DO n=1,nblock + blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock + END DO + remainparts=parts%Nptot-SUM(blocksize) + addedpart=1 + n=nblock/2 + j=1 + DO WHILE(remainparts .GT. 0) + blocksize(n)=blocksize(n)+addedpart + remainparts=remainparts-addedpart + n=n+j + j=-1*(j+SIGN(1,j)) + END DO + DO n=1,nblock + blockstart=blockend+1 + blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot) + ! Initial distribution in z with normalisation between magnetic mirrors + CALL loduni(1,parts%Z(blockstart:blockend)) + parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm + IF(distribtype.eq. 2) THEN + CALL loduni(3,parts%R(blockstart:blockend)) + parts%R(blockstart:blockend)=ra(n)+(rb(n)-ra(n))*parts%R(blockstart:blockend) + ELSE IF(distribtype.eq.3) THEN + CALL lodinvr(3,parts%R(blockstart:blockend),ra(n),rb(n)) + END IF + parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm + END DO + + + + ! 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=abs(P0/rg/me) + Acomp=-SIGN(elchar/me*Athetpos,qsim) + VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp) + END DO + VTHET=VTHET/vnorm + VZ=0._db + VR=0._db + END SUBROUTINE loadDavidson END MODULE beam diff --git a/src/chkrst.f90 b/src/chkrst.f90 index 301437e..fcd4a69 100644 --- a/src/chkrst.f90 +++ b/src/chkrst.f90 @@ -1,121 +1,122 @@ SUBROUTINE chkrst(flag) ! ! Process checkpoint/restart file ! USE basic USE futils USE beam USE fields IMPLICIT NONE INTEGER, INTENT(in) :: flag INTEGER, DIMENSION(:), ALLOCATABLE:: displs INTEGER :: i ! 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) CALL getatt(fidrst, '/Basic', 'cstep', cstep) CALL getatt(fidrst, '/Basic', 'time', time) CALL getatt(fidrst, '/Basic', 'msim', msim) CALL getatt(fidrst, '/Basic', 'qsim', qsim) + CALL getatt(fidrst, '/Basic', 'V', V) CALL getatt(fidrst, '/var0d', 'etot0', etot0) CALL getatt(fidrst, '/var0d', 'epot', epot) CALL getatt(fidrst, '/var0d', 'ekin', ekin) CALL getatt(fidrst, '/var0d', 'etot', etot) - ! Init Electric and Magnetic Fields - CALL fields_init CALL getatt(fidrst, '/Parts', 'nplasma', nplasma) CALL creat_parts(parts,nplasma) CALL getarr(fidrst, '/Parts/Z', parts%Z) CALL getarr(fidrst, '/Parts/R', parts%R) 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 closef(fidrst) WRITE(*,'(3x,a)') & & "Reading from restart file "//TRIM(rstfile)//" completed!" !________________________________________________________________________________ ! 2. Create and write to restart file (DP reals) ! CASE(1) ALLOCATE(displs(0:mpisize-1)) displs(0)=1 DO i=1, mpisize-1 displs(i)=displs(i-1)+Nplocs_all(i-1)-1 END DO IF(mpisize .gt. 1) THEN IF(mpirank .ne. 0) THEN ! Send Particles informations to root process CALL MPI_Gatherv(parts%Zindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & & MPI_INTEGER,0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(parts%Rindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) RETURN ELSE ! Receive particle information from all processes CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & & MPI_INTEGER,0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) END IF + parts%Nptot=sum(Nplocs_all) END IF 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', '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', parts%Nptot) + CALL attach(fidrst, '/Parts', 'nplasma', nplasma) CALL putarr(fidrst, '/Parts/Z', parts%Z) CALL putarr(fidrst, '/Parts/R', parts%R) 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/Zindex', parts%Zindex) CALL putarr(fidrst, '/Parts/Rindex', parts%Rindex) ! ! Fields ! CALL closef(fidrst) WRITE(*,'(3x,a)') & & "Writing to restart file "//TRIM(rstfile)//" completed!" ! END SELECT ! END SUBROUTINE chkrst diff --git a/src/diagnose.f90 b/src/diagnose.f90 index 70b2eb0..b4dd903 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,214 +1,214 @@ 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 mpi IMPLICIT NONE ! INTEGER, INTENT(in) :: kstep ! ! Local vars and arrays INTEGER, PARAMETER :: BUFSIZE = 20 CHARACTER(len=128) :: str, fname INTEGER :: partdensitystatus(MPI_STATUS_SIZE) ! 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) 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/input.", "input") 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%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%pot), "/data/part/pot", "pot") 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(partdensity), "/data/fields/partdensity", "partdensity") IF(mpisize .gt. 1) CALL creatd(fidres, 1, SHAPE(Nplocs_all), "/data/part/Nplocs_all", "Nplocs_all") CALL putarr(fidres, "/data/fields/Br", Br) CALL putarr(fidres, "/data/fields/Bz", Bz) ! ! 1.2 Restart run ! ELSE CALL cp2bk(resfile) ! backup previous result file CALL openf(resfile, fidres) 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 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, status='delete') + 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(modulo(step,ittext).eq. 0 .or. nlend) THEN WRITE(*,'(a,1x,i6.6,a1,i6.6,20x,a,1pe10.3)') & & '*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time WRITE(*,'(a,4(1pe12.4))') 'Epot, Ekin, Etot, Eerr', epot, ekin, etot, etot-etot0 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 htable_endstep(hbuf0) IF(mpisize .gt. 1) CALL append(fidres, "/data/part/Nplocs_all", REAL(Nplocs_all,kind=db)) END IF ! ! 2.2 1d profiles IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL append(fidres, "/data/fields/pot", pot) CALL append(fidres, "/data/fields/Er", Er) CALL append(fidres, "/data/fields/Ez", Ez) CALL MPI_WAIT(ireducerequest(3), partdensitystatus, ierr) CALL append(fidres, "/data/fields/partdensity", partdensity) 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/R", parts%R) CALL append(fidres, "/data/part/Z", parts%Z) CALL append(fidres, "/data/part/UZ", parts%UZ) CALL append(fidres, "/data/part/UR", parts%UR) CALL append(fidres, "/data/part/UTHET", parts%UTHET) CALL append(fidres, "/data/part/pot", parts%pot) !WRITE(*,'(a,6(1pe12.3))') 'R-Z bounds of space = ', rgrid(0), rgrid(nnr(1)), rgrid(nnr(2)+nnr(1)), rgrid(nr), zgrid(0), zgrid(nz) !WRITE(*,'(a,4(1pe12.3))') 'MinMax part. R, Z =', & ! & MINVAL(parts%R), MAXVAL(parts%R), MINVAL(parts%Z), MAXVAL(parts%Z) ! END IF ! ! 2.4 3d profiles ! ! ! 2.5 Xgrafix ! IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0 .AND. mpisize .eq. 1) 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/fields_mod.f90 b/src/fields_mod.f90 index 6625737..f06f7c5 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,324 +1,320 @@ MODULE fields ! ! Module containing all variables and routines used to obtain |E(r,z) and |B(r,z) ! USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, splrz, nlperiod, phinorm, nlPhis, nrank, rhs, mpirank, mpisize USE beam, ONLY: parts USE bsplines USE mumps_bsplines USE mpi IMPLICIT NONE DOUBLE PRECISION, allocatable, SAVE :: matcoef(:,:), sol(:), vec1(:), vec2(:) DOUBLE PRECISION, ALLOCATABLE, SAVE :: fun(:,:), fun2(:,:) TYPE(mumps_mat), SAVE :: femat ! FEM matrix CONTAINS !________________________________________________________________________________ SUBROUTINE fields_init USE basic, ONLY: pot, nlperiod, nrank, rhs USE bsplines USE mumps_bsplines - INTEGER :: nrz(2), i + INTEGER :: nrz(2) ! Set up 2d spline splrz 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))) ALLOCATE(sol(nrank(1)*nrank(2))) ALLOCATE(fun(1:femorder(1)+1,0:1),fun2(1:femorder(2)+1,0:1))!Arrays keeping values of b-splines at gauss node ! Calculate magnetic field components in grid points (Davidson analitical formula employed) 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))) CALL magnet ! Calculate and factorise FEM matrix (depended only on mesh) CALL fematrix CALL factor(femat) END SUBROUTINE fields_init !________________________________________________________________________________ SUBROUTINE rhscon USE bsplines USE mumps_bsplines USE mpi - USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, Zbounds, nplasma + USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma USE beam, ONLY: parts USE mpihelper INTEGER :: zleft, rleft, irow, jcol, it, jw, mu, i - INTEGER :: process, zoverlap, ileftstart - INTEGER, DIMENSION(2):: rhsoverlapstatus(2*MPI_STATUS_SIZE) - INTEGER, DIMENSION(:), ALLOCATABLE:: rhsstatus - ALLOCATE(rhsstatus((mpisize-1)*MPI_STATUS_SIZE)) rhs=0 !Initialise rhs ! Assemble rhs IF (parts%Nptot .ne. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, i, irow, jcol, mu) PRIVATE(fun, fun2) REDUCTION(+:rhs) DO i=1,parts%Nptot CALL locintv(splrz%sp1, parts%Z(i), zleft) CALL locintv(splrz%sp2, parts%R(i), rleft) CALL basfun(parts%Z(i), splrz%sp1, fun, zleft+1) CALL basfun(parts%R(i), splrz%sp2, fun2, rleft+1) DO jw=1,(femorder(2)+1) DO it=1,(femorder(1)+1) irow=zleft+it jcol=rleft+jw mu=irow+(jcol-1)*(nrank(1)) rhs(mu)=rhs(mu)+fun(it,0)*fun2(jw,0)/2/pi*omegap/omegac*V/nplasma END DO END DO END DO !$OMP END PARALLEL DO END IF IF(mpisize .gt. 1) THEN IF (mpirank.eq.0) THEN CALL MPI_REDUCE(MPI_IN_PLACE, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, & & MPI_SUM, 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_REDUCE(rhs, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, & & MPI_SUM, 0, MPI_COMM_WORLD, ierr) END IF END IF ! Dirichlet BC on RHS IF(rgrid(0).ne.0.0_db) rhs(1:nrank(1))=potinn rhs((nrank(2)-1)*nrank(1)+1:nrank(2)*nrank(1))=potout END SUBROUTINE rhscon !========================================================================! SUBROUTINE poisson USE bsplines USE mumps_bsplines USE futils INTEGER:: partend, i, ierr, jder(2) jder(1)=0 jder(2)=0 partend=parts%Nptot IF(mpirank.eq.0) CALL bsolve(femat,rhs,sol) CALL MPI_Bcast(sol(1), nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) matcoef= reshape(sol, (/nrank(1),nrank(2)/)) CALL gridval(splrz, vec1, vec2, pot, jder, matcoef) CALL getgrad(splrz, parts%Z(1:partend), parts%R(1:partend), parts%pot(1:partend), parts%EZ(1:partend), parts%ER(1:partend)) IF (nlphis) THEN Do i=1,partend parts%EZ(i)=-parts%Ez(i) parts%ER(i)=-parts%ER(i) END DO ELSE Do i=1,partend parts%EZ(i)=0 parts%ER(i)=0 END DO END IF IF( mpirank .eq. 0) THEN 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 fematrix USE bsplines USE mumps_bsplines DOUBLE PRECISION, ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:) DOUBLE PRECISION, ALLOCATABLE :: f(:,:), aux(:) DOUBLE PRECISION, ALLOCATABLE :: coefs(:) DOUBLE PRECISION :: 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(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) ALLOCATE(fun(1:femorder(1)+1,0:0),fun2(1:femorder(2)+1,0:0))!Arrays keeping values of b-splines at gauss node END SUBROUTINE fematrix SUBROUTINE fe_dirichlet DOUBLE PRECISION, 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) DOUBLE PRECISION, INTENT(in) :: x(2) INTEGER, INTENT(out) :: idt(:), idw(:) DOUBLE PRECISION, 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 !________________________________________________________________________________ SUBROUTINE magnet USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm USE constants, ONLY: Pi DOUBLE PRECISION :: 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) DOUBLE PRECISION :: bessi0,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/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) DOUBLE PRECISION :: bessi1,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 !________________________________________________________________________________ SUBROUTINE clean_fields Use bsplines DEALLOCATE(matcoef) DEALLOCATE(pot) DEALLOCATE(rhs) DEALLOCATE(sol) DEALLOCATE(fun,fun2) DEALLOCATE(Br,Bz) DEALLOCATE(Er,Ez) DEALLOCATE(vec1,vec2) Call DESTROY_SP(splrz) END SUBROUTINE clean_fields END MODULE fields diff --git a/src/main.f90 b/src/main.f90 index 6228327..a1fa92d 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,84 +1,84 @@ 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) CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, mpisize, 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 stepon CALL tesend + CALL stepon CALL diagnose(step) 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/mv2bk.f90 b/src/mv2bk.f90 index 0720a3d..86207af 100644 --- a/src/mv2bk.f90 +++ b/src/mv2bk.f90 @@ -1,38 +1,36 @@ SUBROUTINE mv2bk(file) ! ! Move file to a backup file if it exists ! IMPLICIT NONE CHARACTER(len=*), INTENT(in) :: file ! LOGICAL :: ex CHARACTER(len=32) :: bkfile - CHARACTER(len=128) :: cmd ! INQUIRE(file=file, exist=ex) IF( ex ) THEN bkfile = TRIM(file) // '_old' CALL move_file(file, LEN_TRIM(file), bkfile, LEN_TRIM(bkfile)) WRITE(*,'(a)') " Moving existing " // TRIM(file) // " to " & & // TRIM(bkfile) END IF END SUBROUTINE mv2bk ! SUBROUTINE cp2bk(file) ! ! Copy file to a backup file if it exists ! IMPLICIT NONE CHARACTER(len=*), INTENT(in) :: file LOGICAL :: ex CHARACTER(len=32) :: bkfile - CHARACTER(len=128) :: cmd ! INQUIRE(file=file, exist=ex) IF( ex ) THEN bkfile = TRIM(file) // '_old' CALL copy_file(file, LEN_TRIM(file), bkfile, LEN_TRIM(bkfile)) WRITE(*,'(a)') " Copying existing " // TRIM(file) // " to " & & // TRIM(bkfile) END IF END SUBROUTINE cp2bk diff --git a/src/resume.f90 b/src/resume.f90 index 6735716..8b8d64a 100644 --- a/src/resume.f90 +++ b/src/resume.f90 @@ -1,23 +1,80 @@ SUBROUTINE resume ! ! Resume from previous run ! Use beam + Use basic Use fields 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) - + CALL chkrst(0) + IF(mpisize .gt. 1) THEN + ALLOCATE(sendcounts(0:mpisize-1), displs(0:mpisize-1)) + IF(mpirank.eq.0) THEN + CALL quicksortparts(1,nplasma) + CALL distribpartsonprocs + CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER,MPI_IN_PLACE, 1,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_DOUBLE_PRECISION,MPI_IN_PLACE,2,MPI_DOUBLE_PRECISION,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 + WRITE(*,*) mpirank, "I am the root process" + 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%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%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) + WRITE(*,*) mpirank, "I am the root process" + ELSE + WRITE(*,*) mpirank, "I am not the root process" + CALL creat_parts(parts,nplasma) + CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER, Nplocs_all(mpirank), 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + sendcounts=2 + CALL MPI_Scatterv(Zbounds,sendcounts,displs,MPI_DOUBLE_PRECISION,Zbounds(mpirank),2,MPI_DOUBLE_PRECISION,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%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%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) + WRITE(*,*) mpirank, "I am not the root process" + END IF + 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) + END IF ! Compute Self consistent electric field CALL bound CALL localisation + ! Init Electric and Magnetic Fields + CALL fields_init CALL rhscon CALL poisson ! END SUBROUTINE resume