MPI fixes

src/main/dtype_kdtree.F90

@@ -84,7 +84,6 @@ subroutine get_mpitype_of_kdnode(dtype)
  integer, parameter              :: ndata = 20
 
  integer, intent(out)            :: dtype
- integer                         :: dtype_old
  integer                         :: nblock, blens(ndata), mpitypes(ndata)
  integer(kind=MPI_ADDRESS_KIND)  :: disp(ndata)
 

src/main/evolve.F90

@@ -174,7 +174,7 @@ subroutine evol(infile,logfile,evfile,dumpfile)
  endif
 #else
  use_global_dt = .true.
- nskip   = npart
+ nskip   = int(ntot)
  nactive = npart
  istepfrac = 0 ! dummy values
  nbinmax   = 0
@@ -251,6 +251,10 @@ subroutine evol(infile,logfile,evfile,dumpfile)
 
     !--print summary of timestep bins
     if (iverbose >= 2) call write_binsummary(npart,nbinmax,dtmax,timeperbin,iphase,ibin,xyzh)
+#else
+    !--If not using individual timestepping, set nskip to the total number of particles
+    !  across all nodes
+    nskip = int(ntot)
 #endif
 
     if (gravity .and. icreate_sinks > 0 .and. ipart_rhomax /= 0) then
@@ -307,7 +311,7 @@ subroutine evol(infile,logfile,evfile,dumpfile)
     if (abs(tcheck-time) > 1.e-4) call warning('evolve','time out of sync',var='error',val=abs(tcheck-time))
 
     if (id==master .and. (iverbose >= 1 .or. inbin <= 3)) &
-       call print_dtlog_ind(iprint,istepfrac,2**nbinmaxprev,time,dt,nactivetot,tcpu2-tcpu1,npart)
+       call print_dtlog_ind(iprint,istepfrac,2**nbinmaxprev,time,dt,nactivetot,tcpu2-tcpu1,ntot)
 
     !--if total number of bins has changed, adjust istepfrac and dt accordingly
     !  (ie., decrease or increase the timestep)
@@ -333,7 +337,7 @@ subroutine evol(infile,logfile,evfile,dumpfile)
 !
 !--write log every step (NB: must print after dt has been set in order to identify timestep constraint)
 !
-    if (id==master) call print_dtlog(iprint,time,dt,dtforce,dtcourant,dterr,dtmax,dtrad,dtprint,dtinject,npart)
+    if (id==master) call print_dtlog(iprint,time,dt,dtforce,dtcourant,dterr,dtmax,dtrad,dtprint,dtinject,ntot)
 #endif
 
 !   check that MPI threads are synchronised in time

src/main/initial.F90

@@ -131,7 +131,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  use readwrite_infile, only:read_infile,write_infile
  use readwrite_dumps,  only:read_dump,write_fulldump
  use part,             only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,&
-                            npartoftype,maxtypes,ndusttypes,alphaind,ntot,ndim, &
+                            npartoftype,maxtypes,ndusttypes,alphaind,ntot,ndim,update_npartoftypetot,&
                             maxphase,iphase,isetphase,iamtype, &
                             nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,igas,idust,massoftype,&
                             epot_sinksink,get_ntypes,isdead_or_accreted,dustfrac,ddustevol,&
@@ -225,7 +225,6 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
  character(len=*), intent(out) :: logfile,evfile,dumpfile
  logical,          intent(in), optional :: noread
  integer         :: ierr,i,j,nerr,nwarn,ialphaloc,merge_n,merge_ij(maxptmass)
- integer(kind=8) :: npartoftypetot(maxtypes)
  real            :: poti,dtf,hfactfile,fextv(3)
  real            :: hi,pmassi,rhoi1
  real            :: dtsinkgas,dtsinksink,fonrmax,dtphi2,dtnew_first,dtinject
@@ -305,7 +304,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
 !--get total number of particles (on all processors)
 !
  ntot           = reduceall_mpi('+',npart)
- npartoftypetot = reduce_mpi('+',npartoftype)
+ call update_npartoftypetot
  if (id==master) write(iprint,"(a,i12)") ' npart total   = ',ntot
  if (npart > 0) then
     if (id==master .and. maxalpha==maxp)  write(iprint,*) 'mean alpha  initial: ',sum(alphaind(1,1:npart))/real(npart)
@@ -429,7 +428,6 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     call initialise_externalforces(iexternalforce,ierr)
     if (ierr /= 0) call fatal('initial','error in external force settings/initialisation')
     call get_grforce_all(npart,xyzh,metrics,metricderivs,vxyzu,dens,fext,dtextforce)
-    write(iprint,*) 'dt(extforce)  = ',dtextforce
  endif
 #else
  if (iexternalforce > 0) then
@@ -451,10 +449,14 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
        endif
     enddo
     !$omp end parallel do
-    write(iprint,*) 'dt(extforce)  = ',dtextforce
  endif
 #endif
 
+ if (iexternalforce > 0) then
+    dtextforce = reduceall_mpi('min',dtextforce)
+    if (id==master) write(iprint,*) 'dt(extforce)  = ',dtextforce
+ endif
+
 !
 !-- Set external force to zero on boundary particles
 !
@@ -643,6 +645,14 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     endif
  enddo
  !$omp end parallel do
+
+ xmin = reduceall_mpi('min',xmin)
+ ymin = reduceall_mpi('min',ymin)
+ zmin = reduceall_mpi('min',zmin)
+ xmax = reduceall_mpi('max',xmax)
+ ymax = reduceall_mpi('max',ymax)
+ zmax = reduceall_mpi('max',zmax)
+
  dx = abs(xmax - xmin)
  dy = abs(ymax - ymin)
  dz = abs(zmax - zmin)

src/main/mpi_balance.F90

@@ -96,7 +96,7 @@ end subroutine balance_init
 !----------------------------------------------------------------
 subroutine balancedomains(npart)
  use io,     only:id,master,iverbose,fatal
- use part,   only:shuffle_part,count_dead_particles,ibelong
+ use part,   only:shuffle_part,count_dead_particles,ibelong,update_npartoftypetot
  use timing, only:getused,printused
  use mpiutils, only:barrier_mpi
  implicit none
@@ -144,6 +144,10 @@ subroutine balancedomains(npart)
     print*,id,'ntot_start',ntot_start
     call fatal('balance','number of particles before and after balance not equal')
  endif
+
+ !  Update particle types
+ call update_npartoftypetot
+
  if (id==master .and. iverbose >= 3) call printused(tstart)
 
  return

src/main/mpi_dens.F90

@@ -67,13 +67,13 @@ module mpidens
 #ifdef MPI
 subroutine get_mpitype_of_celldens(dtype)
  use mpi
+ use io, only:error
 
  integer, parameter              :: ndata = 20
 
  integer, intent(out)            :: dtype
- integer                         :: dtype_old
  integer                         :: nblock, blens(ndata), mpitypes(ndata)
- integer(kind=4)                 :: disp(ndata)
+ integer(kind=MPI_ADDRESS_KIND)  :: disp(ndata)
 
  type(celldens)                 :: cell
  integer(kind=MPI_ADDRESS_KIND)  :: addr,start,lb,extent
@@ -196,20 +196,14 @@ subroutine get_mpitype_of_celldens(dtype)
  call MPI_GET_ADDRESS(cell%pad,addr,mpierr)
  disp(nblock) = addr - start
 
- call MPI_TYPE_STRUCT(nblock,blens(1:nblock),disp(1:nblock),mpitypes(1:nblock),dtype,mpierr)
+ call MPI_TYPE_CREATE_STRUCT(nblock,blens(1:nblock),disp(1:nblock),mpitypes(1:nblock),dtype,mpierr)
  call MPI_TYPE_COMMIT(dtype,mpierr)
 
  ! check extent okay
  call MPI_TYPE_GET_EXTENT(dtype,lb,extent,mpierr)
  if (extent /= sizeof(cell)) then
-    dtype_old = dtype
-    lb = 0
-    extent = sizeof(cell)
-    call MPI_TYPE_CREATE_RESIZED(dtype_old,lb,extent,dtype,mpierr)
-    call MPI_TYPE_COMMIT(dtype,mpierr)
-    call MPI_TYPE_FREE(dtype_old,mpierr)
+    call error('mpi_dens','MPI_TYPE_GET_EXTENT has calculated the extent incorrectly')
  endif
-
 end subroutine get_mpitype_of_celldens
 #endif
 

src/main/mpi_force.F90

@@ -68,13 +68,13 @@ module mpiforce
 #ifdef MPI
 subroutine get_mpitype_of_cellforce(dtype)
  use mpi
+ use io, only:error
 
  integer, parameter              :: ndata = 20
 
  integer, intent(out)            :: dtype
- integer                         :: dtype_old
  integer                         :: nblock, blens(ndata), mpitypes(ndata)
- integer(kind=4)                 :: disp(ndata)
+ integer(kind=MPI_ADDRESS_KIND)  :: disp(ndata)
 
  type(cellforce)                 :: cell
  integer(kind=MPI_ADDRESS_KIND)  :: addr,start,lb,extent
@@ -203,18 +203,13 @@ subroutine get_mpitype_of_cellforce(dtype)
  call MPI_GET_ADDRESS(cell%pad,addr,mpierr)
  disp(nblock) = addr - start
 
- call MPI_TYPE_STRUCT(nblock,blens(1:nblock),disp(1:nblock),mpitypes(1:nblock),dtype,mpierr)
+ call MPI_TYPE_CREATE_STRUCT(nblock,blens(1:nblock),disp(1:nblock),mpitypes(1:nblock),dtype,mpierr)
  call MPI_TYPE_COMMIT(dtype,mpierr)
 
  ! check extent okay
  call MPI_TYPE_GET_EXTENT(dtype,lb,extent,mpierr)
  if (extent /= sizeof(cell)) then
-    dtype_old = dtype
-    lb = 0
-    extent = sizeof(cell)
-    call MPI_TYPE_CREATE_RESIZED(dtype_old,lb,extent,dtype,mpierr)
-    call MPI_TYPE_COMMIT(dtype,mpierr)
-    call MPI_TYPE_FREE(dtype_old,mpierr)
+    call error('mpi_force','MPI_TYPE_GET_EXTENT has calculated the extent incorrectly')
  endif
 
 end subroutine get_mpitype_of_cellforce

src/main/part.F90

@@ -361,8 +361,9 @@ module part
  integer(kind=8) :: ntot
  integer         :: ideadhead = 0
 
- integer :: npartoftype(maxtypes)
- real    :: massoftype(maxtypes)
+ integer         :: npartoftype(maxtypes)
+ integer(kind=8) :: npartoftypetot(maxtypes)
+ real            :: massoftype(maxtypes)
 
  integer :: ndustsmall,ndustlarge,ndusttypes
 !
@@ -564,6 +565,7 @@ subroutine init_part
  npart = 0
  nptmass = 0
  npartoftype(:) = 0
+ npartoftypetot(:) = 0
  massoftype(:)  = 0.
 !--initialise point mass arrays to zero
  xyzmh_ptmass = 0.
@@ -1803,4 +1805,11 @@ real function Omega_k(i)
 
 end function Omega_k
 
+subroutine update_npartoftypetot
+ use mpiutils, only:reduceall_mpi
+
+ npartoftypetot = reduceall_mpi('+',npartoftype)
+
+end subroutine update_npartoftypetot
+
 end module part

src/main/readwrite_dumps_common.F90

@@ -118,7 +118,7 @@ end subroutine get_options_from_fileid
 !  and perform basic sanity checks
 !+
 !---------------------------------------------------------------
-subroutine check_arrays(i1,i2,npartoftype,npartread,nptmass,nsinkproperties,massoftype,&
+subroutine check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkproperties,massoftype,&
                         alphafile,tfile,phantomdump,got_iphase,got_xyzh,got_vxyzu,got_alpha, &
                         got_krome_mols,got_krome_gamma,got_krome_mu,got_krome_T,got_x,got_z,got_mu, &
                         got_abund,got_dustfrac,got_sink_data,got_sink_vels,got_Bxyz,got_psi,got_dustprop,got_pxyzu,got_VrelVf, &
@@ -132,7 +132,7 @@ subroutine check_arrays(i1,i2,npartoftype,npartread,nptmass,nsinkproperties,mass
  use io,   only:warning,id,master
  use options,    only:alpha,use_dustfrac,use_var_comp
  use sphNGutils, only:itype_from_sphNG_iphase,isphNG_accreted
- integer,         intent(in)    :: i1,i2,npartoftype(:),npartread,nptmass,nsinkproperties
+ integer,         intent(in)    :: i1,i2,noffset,npartoftype(:),npartread,nptmass,nsinkproperties
  real,            intent(in)    :: massoftype(:),alphafile,tfile
  logical,         intent(in)    :: phantomdump,got_iphase,got_xyzh(:),got_vxyzu(:),got_alpha,got_dustprop(:)
  logical,         intent(in)    :: got_VrelVf,got_dustgasprop(:),got_x,got_z,got_mu
@@ -372,7 +372,7 @@ subroutine check_arrays(i1,i2,npartoftype,npartread,nptmass,nsinkproperties,mass
 !
  if (.not.got_iorig) then
     do i=i1,i2
-       iorig(i) = i
+       iorig(i) = i + noffset
     enddo
     norig = i2
     if (id==master .and. i1==1) write(*,*) 'WARNING: Particle IDs not in dump; resetting IDs'

src/main/readwrite_dumps_fortran.F90

@@ -211,7 +211,8 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG)
                  lightcurve,store_temperature,use_dustgrowth,store_dust_temperature,gr
  use eos,   only:ieos,eos_is_non_ideal,eos_outputs_mu
  use io,    only:idump,iprint,real4,id,master,error,warning,nprocs
- use part,  only:xyzh,xyzh_label,vxyzu,vxyzu_label,Bevol,Bevol_label,Bxyz,Bxyz_label,npart,npartoftype,maxtypes, &
+ use part,  only:xyzh,xyzh_label,vxyzu,vxyzu_label,Bevol,Bevol_label,Bxyz,Bxyz_label,npart,maxtypes, &
+                 npartoftype,npartoftypetot,update_npartoftypetot, &
                  alphaind,rhoh,divBsymm,maxphase,iphase,iamtype_int1,iamtype_int11, &
                  nptmass,nsinkproperties,xyzmh_ptmass,xyzmh_ptmass_label,vxyz_ptmass,vxyz_ptmass_label,&
                  maxptmass,get_pmass,h2chemistry,nabundances,abundance,abundance_label,mhd,&
@@ -256,7 +257,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG)
  integer            :: ipass,k,l
  integer            :: ierr,ierrs(29)
  integer            :: nblocks,nblockarrays,narraylengths
- integer(kind=8)    :: nparttot,npartoftypetot(maxtypes)
+ integer(kind=8)    :: nparttot
  logical            :: sphNGdump,write_itype,use_gas
  character(len=lenid)  :: fileid
  type(dump_h)          :: hdr
@@ -267,17 +268,17 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG)
 !--allow non-MPI calls to create MPI dump files
 #ifdef MPI
  nparttot = reduceall_mpi('+',npart)
- npartoftypetot = reduceall_mpi('+',npartoftype)
+ call update_npartoftypetot
 #else
  if (present(ntotal)) then
     nparttot = ntotal
-    npartoftypetot = npartoftype
+    call update_npartoftypetot
     if (all(npartoftypetot==0)) then
        npartoftypetot(1) = ntotal
     endif
  else
     nparttot = npart
-    npartoftypetot = npartoftype
+    call update_npartoftypetot
  endif
 #endif
  nblocks = nprocs
@@ -503,7 +504,8 @@ end subroutine write_fulldump_fortran
 subroutine write_smalldump_fortran(t,dumpfile)
  use dim,        only:maxp,maxtypes,use_dust,lightcurve,use_dustgrowth
  use io,         only:idump,iprint,real4,id,master,error,warning,nprocs
- use part,       only:xyzh,xyzh_label,npart,npartoftype,Bxyz,Bxyz_label,&
+ use part,       only:xyzh,xyzh_label,npart,Bxyz,Bxyz_label,&
+                      npartoftypetot,update_npartoftypetot,&
                       maxphase,iphase,h2chemistry,nabundances,&
                       nptmass,nsinkproperties,xyzmh_ptmass,xyzmh_ptmass_label,&
                       abundance,abundance_label,mhd,dustfrac,iamtype_int11,&
@@ -521,14 +523,14 @@ subroutine write_smalldump_fortran(t,dumpfile)
  integer         :: nums(ndatatypes,4)
  integer         :: ierr,ipass,k
  integer         :: nblocks,nblockarrays,narraylengths
- integer(kind=8) :: nparttot,npartoftypetot(maxtypes)
+ integer(kind=8) :: nparttot
  logical         :: write_itype
  type(dump_h)    :: hdr
 !
 !--collect global information from MPI threads
 !
  nparttot = reduceall_mpi('+',npart)
- npartoftypetot = reduceall_mpi('+',npartoftype)
+ call update_npartoftypetot
  nblocks = nprocs
 
  narraylengths = 2
@@ -777,7 +779,7 @@ subroutine read_dump_fortran(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ie
 #ifdef INJECT_PARTICLES
        call allocate_memory(maxp_hard)
 #else
-       call allocate_memory(int(min(nprocs,3)*nparttot/nprocs))
+       call allocate_memory(int(min(nprocs,4)*nparttot/nprocs))
 #endif
     endif
 !
@@ -1273,7 +1275,7 @@ subroutine read_phantom_arrays(i1,i2,noffset,narraylengths,nums,npartread,nparto
  !
  ! check for errors
  !
- call check_arrays(i1,i2,npartoftype,npartread,nptmass,nsinkproperties,massoftype,&
+ call check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkproperties,massoftype,&
                    alphafile,tfile,phantomdump,got_iphase,got_xyzh,got_vxyzu,got_alpha, &
                    got_krome_mols,got_krome_gamma,got_krome_mu,got_krome_T,got_x,got_z,got_mu, &
                    got_abund,got_dustfrac,got_sink_data,got_sink_vels,got_Bxyz,got_psi,got_dustprop,got_pxyzu,got_VrelVf, &
@@ -1352,7 +1354,8 @@ subroutine unfill_header(hdr,phantomdump,got_tags,nparttot, &
  use dim,        only:maxdustlarge,use_dust
  use io,         only:master ! check this
  use eos,        only:isink
- use part,       only:maxtypes,igas,idust,ndustsmall,ndustlarge,ndusttypes
+ use part,       only:maxtypes,igas,idust,ndustsmall,ndustlarge,ndusttypes,&
+                      npartoftypetot
  use units,      only:udist,umass,utime,set_units_extra,set_units
  use dump_utils, only:extract,dump_h
  use fileutils,  only:make_tags_unique
@@ -1365,7 +1368,7 @@ subroutine unfill_header(hdr,phantomdump,got_tags,nparttot, &
  integer,         intent(out) :: ierr
  integer         :: nparttoti,npartoftypetoti(maxtypes),ntypesinfile,nptinfile
  integer         :: ierr1,ierrs(3),i,counter
- integer(kind=8) :: npartoftypetot(maxtypes),ntypesinfile8
+ integer(kind=8) :: ntypesinfile8
  character(len=10) :: dust_label(maxdustlarge)
 
  ierr = 0