From 4c1244af8d46c6c510722415bdfbc613185257c2 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 14:46:28 +1100 Subject: [PATCH 01/21] fix: assign unique IDs when reading a file without IDs --- src/main/readwrite_dumps_common.F90 | 6 +++--- src/main/readwrite_dumps_fortran.F90 | 2 +- src/main/readwrite_dumps_hdf5.F90 | 1 + 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/main/readwrite_dumps_common.F90 b/src/main/readwrite_dumps_common.F90 index 46e521542..3d0632acb 100644 --- a/src/main/readwrite_dumps_common.F90 +++ b/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' diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 89be85f01..e3132d51b 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -1273,7 +1273,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, & diff --git a/src/main/readwrite_dumps_hdf5.F90 b/src/main/readwrite_dumps_hdf5.F90 index 47f132cc6..892f7d6c8 100644 --- a/src/main/readwrite_dumps_hdf5.F90 +++ b/src/main/readwrite_dumps_hdf5.F90 @@ -707,6 +707,7 @@ subroutine read_any_dump_hdf5( if (.not.smalldump) then call check_arrays(1, & npart, & + 0, & npartoftype, & npart, & nptmass, & From 69e0a92727d0da1a86b8d901585d8422caa5f556 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 15:17:15 +1100 Subject: [PATCH 02/21] reduce npartoftype(i) in printout when writing header --- src/main/writeheader.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/main/writeheader.F90 b/src/main/writeheader.F90 index 1ff19b31d..8b7fd834e 100644 --- a/src/main/writeheader.F90 +++ b/src/main/writeheader.F90 @@ -15,8 +15,8 @@ module writeheader ! :Runtime parameters: None ! ! :Dependencies: boundary, cooling, dim, dust, eos, gitinfo, growth, io, -! kernel, metric_tools, options, part, physcon, readwrite_infile, units, -! viscosity +! kernel, metric_tools, mpiutils, options, part, physcon, +! readwrite_infile, units, viscosity ! implicit none public :: write_header,write_codeinfo @@ -82,6 +82,7 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) use part, only:hfact,massoftype,mhd,& gravity,h2chemistry,periodic,npartoftype,massoftype,& labeltype,maxtypes + use mpiutils, only:reduceall_mpi use eos, only:eosinfo use cooling, only:cooling_implicit,cooling_explicit,Tfloor,ufloor use readwrite_infile, only:write_infile @@ -94,7 +95,7 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) #ifdef GR use metric_tools, only:print_metricinfo #endif - integer :: Nneigh,i + integer :: Nneigh,i,npartoftypetoti integer, intent(in) :: icall character(len=*), intent(in) :: infile,evfile,logfile,dumpfile integer(kind=8), intent(in), optional :: ntot @@ -133,9 +134,10 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) if (present(ntot)) then write(iprint,"(/,' Number of particles = ',i12)") ntot do i = 1,maxtypes - if (npartoftype(i) > 0) then + npartoftypetoti = reduceall_mpi('+', npartoftype(i)) + if (npartoftypetoti > 0) then write(iprint,"(1x,3a,i12,a,es14.6)") & - "Number & mass of ",labeltype(i)," particles: ", npartoftype(i),", ",massoftype(i) + "Number & mass of ",labeltype(i)," particles: ", npartoftypetoti,", ",massoftype(i) endif enddo write(iprint,"(a)") " " From 3431b00a0c318f3894873563b538d13779e3a548 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 15:32:09 +1100 Subject: [PATCH 03/21] store npartoftypetot as a global variable rather than reducing it on the fly everywhere --- src/main/initial.F90 | 5 ++--- src/main/mpi_balance.F90 | 6 +++++- src/main/part.F90 | 13 +++++++++++-- src/main/readwrite_dumps_fortran.F90 | 19 +++++++++++-------- src/main/readwrite_dumps_hdf5.F90 | 7 ++++--- src/main/step_leapfrog.F90 | 6 +++--- src/main/writeheader.F90 | 6 +++--- 7 files changed, 39 insertions(+), 23 deletions(-) diff --git a/src/main/initial.F90 b/src/main/initial.F90 index 114d21744..bb79d639e 100644 --- a/src/main/initial.F90 +++ b/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) diff --git a/src/main/mpi_balance.F90 b/src/main/mpi_balance.F90 index a262fe165..89a7de08c 100644 --- a/src/main/mpi_balance.F90 +++ b/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 diff --git a/src/main/part.F90 b/src/main/part.F90 index dabd0a8b0..fb786d357 100644 --- a/src/main/part.F90 +++ b/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 + end module part diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 89be85f01..8b1ba31f1 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/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, & + 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,7 +268,7 @@ 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 @@ -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 @@ -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 diff --git a/src/main/readwrite_dumps_hdf5.F90 b/src/main/readwrite_dumps_hdf5.F90 index 47f132cc6..d7184db27 100644 --- a/src/main/readwrite_dumps_hdf5.F90 +++ b/src/main/readwrite_dumps_hdf5.F90 @@ -109,7 +109,8 @@ subroutine write_dump_hdf5(t,dumpfile,fulldump,ntotal,dtind) Bextz,ndustlarge,idust,idustbound,grainsize, & graindens,h2chemistry,lightcurve,ndivcurlB, & ndivcurlv,pxyzu,dens,gamma_chem,mu_chem,T_gas_cool, & - dust_temp,rad,radprop,itemp,igasP,eos_vars,iorig + dust_temp,rad,radprop,itemp,igasP,eos_vars,iorig, & + npartoftypetot,update_npartoftypetot #ifdef NUCLEATION use part, only:nucleation #endif @@ -135,7 +136,7 @@ subroutine write_dump_hdf5(t,dumpfile,fulldump,ntotal,dtind) integer :: i integer :: ierr - integer(kind=8) :: nparttot,npartoftypetot(maxtypes) + integer(kind=8) :: nparttot logical :: ind_timesteps,const_av,prdrag,isothermal real, allocatable :: dtin(:),beta_pr(:) character(len=200) :: fileid,fstr,sstr @@ -167,7 +168,7 @@ subroutine write_dump_hdf5(t,dumpfile,fulldump,ntotal,dtind) !--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 diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 0060ff044..ce7c0c762 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -101,7 +101,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) use part, only:xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol, & rad,drad,radprop,isdead_or_accreted,rhoh,dhdrho,& iphase,iamtype,massoftype,maxphase,igas,idust,mhd,& - iamboundary,get_ntypes,npartoftype,& + iamboundary,get_ntypes,npartoftypetot,& dustfrac,dustevol,ddustevol,eos_vars,alphaind,nptmass,& dustprop,ddustprop,dustproppred,ndustsmall,pxyzu,dens,metrics,ics use cooling, only:cooling_implicit,ufloor @@ -172,7 +172,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) ! velocity predictor step, using dtsph !-------------------------------------- itype = igas - ntypes = get_ntypes(reduceall_mpi('+',npartoftype)) + ntypes = get_ntypes(npartoftypetot) pmassi = massoftype(itype) store_itype = (maxphase==maxp .and. ntypes > 1) ialphaloc = 2 @@ -422,7 +422,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) np = 0 itype = igas pmassi = massoftype(igas) - ntypes = get_ntypes(npartoftype) + ntypes = get_ntypes(npartoftypetot) store_itype = (maxphase==maxp .and. ntypes > 1) !$omp parallel default(none) & !$omp shared(xyzh,vxyzu,vpred,fxyzu,npart,hdtsph,store_itype) & diff --git a/src/main/writeheader.F90 b/src/main/writeheader.F90 index 8b7fd834e..ca4087fb2 100644 --- a/src/main/writeheader.F90 +++ b/src/main/writeheader.F90 @@ -81,6 +81,7 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) use options, only:tolh,alpha,alphau,alphaB,ieos,alphamax,use_dustfrac use part, only:hfact,massoftype,mhd,& gravity,h2chemistry,periodic,npartoftype,massoftype,& + npartoftypetot,& labeltype,maxtypes use mpiutils, only:reduceall_mpi use eos, only:eosinfo @@ -134,10 +135,9 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) if (present(ntot)) then write(iprint,"(/,' Number of particles = ',i12)") ntot do i = 1,maxtypes - npartoftypetoti = reduceall_mpi('+', npartoftype(i)) - if (npartoftypetoti > 0) then + if (npartoftypetot(i) > 0) then write(iprint,"(1x,3a,i12,a,es14.6)") & - "Number & mass of ",labeltype(i)," particles: ", npartoftypetoti,", ",massoftype(i) + "Number & mass of ",labeltype(i)," particles: ", npartoftypetot(i),", ",massoftype(i) endif enddo write(iprint,"(a)") " " From 5dad02b6f951e9330c9b1f97f7eb96680c0b4d0b Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 15:35:06 +1100 Subject: [PATCH 04/21] amend previous commit npartoftype is still required --- src/main/readwrite_dumps_fortran.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 8b1ba31f1..7e80f4a8b 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -212,7 +212,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) 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,maxtypes, & - npartoftypetot,update_npartoftypetot, & + 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,& From 82e908793dbb74ff0600d1921a839b3a3e27ed4a Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:11:50 +1100 Subject: [PATCH 05/21] use update_npartoftypetot instead of copying directly --- src/main/readwrite_dumps_fortran.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 7e80f4a8b..bfc9fb450 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -272,13 +272,13 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) #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 From 7e180b24921a4f4ed3340331674c84f1e7602571 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:15:24 +1100 Subject: [PATCH 06/21] reduce summary variables using MPI --- src/main/step_leapfrog.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 0060ff044..eb8ad6b22 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -663,6 +663,13 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) if (gr) vxyzu = vpred ! May need primitive variables elsewhere? endif enddo iterations + + ! MPI reduce summary variables + nwake = int(reduceall_mpi('+', nwake)) + nvfloorp = int(reduceall_mpi('+', nvfloorp)) + nvfloorps = int(reduceall_mpi('+', nvfloorps)) + nvfloorc = int(reduceall_mpi('+', nvfloorc)) + ! Summary statements & crash if velocity is not converged if (nwake > 0) call summary_variable('wake', iowake, 0,real(nwake) ) if (nvfloorp > 0) call summary_variable('floor',iosumflrp, 0,real(nvfloorp) ) From 898df99405f6d0a261c234f38732a4d7957c50de Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:26:00 +1100 Subject: [PATCH 07/21] replace mpi_type_struct with mpi_type_create_struct MPI_TYPE_STRUCT was deprecated since MPI 2.0 and replaced with MPI_TYPE_CREATE_STRUCT It has exactly the same interface and does exactly the same thing And then MPI_TYPE_STRUCT was removed from the standard in MPI 3.0 But MPI_TYPE_STRUCT still works, but gives the incorrect result --- src/main/mpi_dens.F90 | 14 ++++---------- src/main/mpi_force.F90 | 13 ++++--------- 2 files changed, 8 insertions(+), 19 deletions(-) diff --git a/src/main/mpi_dens.F90 b/src/main/mpi_dens.F90 index 149d012ab..1dce490d1 100644 --- a/src/main/mpi_dens.F90 +++ b/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 diff --git a/src/main/mpi_force.F90 b/src/main/mpi_force.F90 index 351666a2d..da50b516c 100644 --- a/src/main/mpi_force.F90 +++ b/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 From d3cbab0bf7b764e962af2b98850ba13734391988 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:32:12 +1100 Subject: [PATCH 08/21] reduce dtextforce across mpi tasks and only print on master --- src/main/initial.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/main/initial.F90 b/src/main/initial.F90 index 114d21744..2ab8630e2 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -429,7 +429,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 +450,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 ! From c5f469db62d2d351599f78c58f0866b7babfd370 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:33:59 +1100 Subject: [PATCH 09/21] remove unused variable declaration --- src/main/dtype_kdtree.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/main/dtype_kdtree.F90 b/src/main/dtype_kdtree.F90 index a537ea7f2..fcb75a157 100644 --- a/src/main/dtype_kdtree.F90 +++ b/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) From c61278e4cc2523d0999147750297678d3c34d0ed Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:36:17 +1100 Subject: [PATCH 10/21] initial: reduce box size across mpi tasks --- src/main/initial.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/main/initial.F90 b/src/main/initial.F90 index 2ab8630e2..fa9b88f77 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -646,6 +646,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) From c9909dd158433c46dddd07db483e8947f2d9b06e Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:39:14 +1100 Subject: [PATCH 11/21] use ntot instead of npart for nskip npart is the number of particles on the local task ntot is the total number of particles --- src/main/evolve.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index 84ce7c862..2ab01e7c3 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -174,7 +174,7 @@ subroutine evol(infile,logfile,evfile,dumpfile) endif #else use_global_dt = .true. - nskip = npart + nskip = 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 From cfdcf8520df7c56998e3238f0c8e0c9ed8999ef3 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Fri, 18 Feb 2022 16:41:07 +1100 Subject: [PATCH 12/21] use int() to convert between int types --- src/main/evolve.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index 2ab01e7c3..fb1452cfe 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -174,7 +174,7 @@ subroutine evol(infile,logfile,evfile,dumpfile) endif #else use_global_dt = .true. - nskip = ntot + nskip = int(ntot) nactive = npart istepfrac = 0 ! dummy values nbinmax = 0 From 46e31ce0aa8dbf980b6645b0059dc32ccfaec9c7 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 12:33:25 +1100 Subject: [PATCH 13/21] [indent-bot] standardised indentation --- src/main/mpi_force.F90 | 2 +- src/main/part.F90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/main/mpi_force.F90 b/src/main/mpi_force.F90 index da50b516c..316130234 100644 --- a/src/main/mpi_force.F90 +++ b/src/main/mpi_force.F90 @@ -209,7 +209,7 @@ subroutine get_mpitype_of_cellforce(dtype) ! check extent okay call MPI_TYPE_GET_EXTENT(dtype,lb,extent,mpierr) if (extent /= sizeof(cell)) then - call error('mpi_force','MPI_TYPE_GET_EXTENT has calculated the extent incorrectly') + call error('mpi_force','MPI_TYPE_GET_EXTENT has calculated the extent incorrectly') endif end subroutine get_mpitype_of_cellforce diff --git a/src/main/part.F90 b/src/main/part.F90 index fb786d357..34241bff6 100644 --- a/src/main/part.F90 +++ b/src/main/part.F90 @@ -1806,10 +1806,10 @@ real function Omega_k(i) end function Omega_k subroutine update_npartoftypetot - use mpiutils, only:reduceall_mpi + use mpiutils, only:reduceall_mpi - npartoftypetot = reduceall_mpi('+',npartoftype) + npartoftypetot = reduceall_mpi('+',npartoftype) -end subroutine +end subroutine update_npartoftypetot end module part From 8e07f1eb7d20e10581dafc92081d51e5f64e7518 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 12:39:43 +1100 Subject: [PATCH 14/21] use ntot instead of npart to print dtlog --- src/main/evolve.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index fb1452cfe..e0e766aed 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -311,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) @@ -337,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 From 19b114897d71b70987ae8cf1fb20092c8c3ecd53 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 12:43:01 +1100 Subject: [PATCH 15/21] change particle allocation ratio maximum memory allocated is 4x the uniform share, make change on hdf5 version as well --- src/main/readwrite_dumps_fortran.F90 | 2 +- src/main/readwrite_dumps_hdf5.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 4391e78e4..c7e23f8d0 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -779,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 ! diff --git a/src/main/readwrite_dumps_hdf5.F90 b/src/main/readwrite_dumps_hdf5.F90 index 0212bee10..e55500060 100644 --- a/src/main/readwrite_dumps_hdf5.F90 +++ b/src/main/readwrite_dumps_hdf5.F90 @@ -620,7 +620,7 @@ subroutine read_any_dump_hdf5( #ifdef INJECT_PARTICLES call allocate_memory(maxp_hard) #else - call allocate_memory(int(npart / nprocs) + 1) + call allocate_memory(int(min(nprocs,4)*nparttot/nprocs)) #endif if (periodic) then From 7dd759426a4846e25b68248dde757ac52d18c7c5 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 12:44:29 +1100 Subject: [PATCH 16/21] optimisations to step_leapfrog don't do ptmass reductions unless necessary --- src/main/step_leapfrog.F90 | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 9b90b908a..5f3be1286 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -1355,7 +1355,7 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time, ! ! reduction of sink-gas forces from each MPI thread ! - call reduce_in_place_mpi('+',fxyz_ptmass(:,1:nptmass)) + if (nptmass > 0) call reduce_in_place_mpi('+',fxyz_ptmass(:,1:nptmass)) !--------------------------- ! corrector during substeps @@ -1443,16 +1443,18 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time, ! ! reduction of sink particle changes across MPI ! - call reduce_in_place_mpi('+',dptmass(:,1:nptmass)) + if (nptmass > 0) then + call reduce_in_place_mpi('+',dptmass(:,1:nptmass)) - naccreted = int(reduceall_mpi('+',naccreted)) - nfail = int(reduceall_mpi('+',nfail)) + naccreted = int(reduceall_mpi('+',naccreted)) + nfail = int(reduceall_mpi('+',nfail)) - if (id==master) call update_ptmass(dptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass) + if (id==master) call update_ptmass(dptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass) - call bcast_mpi(xyzmh_ptmass(:,1:nptmass)) - call bcast_mpi(vxyz_ptmass(:,1:nptmass)) - call bcast_mpi(fxyz_ptmass(:,1:nptmass)) + call bcast_mpi(xyzmh_ptmass(:,1:nptmass)) + call bcast_mpi(vxyz_ptmass(:,1:nptmass)) + call bcast_mpi(fxyz_ptmass(:,1:nptmass)) + endif if (iverbose >= 2 .and. id==master .and. naccreted /= 0) write(iprint,"(a,es10.3,a,i4,a,i4,a)") & 'Step: at time ',timei,', ',naccreted,' particles were accreted amongst ',nptmass,' sink(s).' From 34b7ec13aff5e2503fc250de2ac1797fd38b3ed6 Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 12:46:22 +1100 Subject: [PATCH 17/21] use integer(kind=8) for np --- src/main/timestep.F90 | 8 ++++---- src/main/utils_indtimesteps.F90 | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/main/timestep.F90 b/src/main/timestep.F90 index 96ce19a15..dded777fb 100644 --- a/src/main/timestep.F90 +++ b/src/main/timestep.F90 @@ -68,10 +68,10 @@ end subroutine set_defaults_timestep !----------------------------------------------------------------- subroutine print_dtlog(iprint,time,dt,dtforce,dtcourant,dterr,dtmax,& dtrad,dtprint,dtinj,np) - integer, intent(in) :: iprint - real, intent(in) :: time,dt,dtforce,dtcourant,dterr,dtmax,dtrad - real, intent(in), optional :: dtprint,dtinj - integer, intent(in) :: np + integer, intent(in) :: iprint + real, intent(in) :: time,dt,dtforce,dtcourant,dterr,dtmax,dtrad + real, intent(in), optional :: dtprint,dtinj + integer(kind=8), intent(in) :: np character(len=20) :: str integer, save :: nplast = 0 diff --git a/src/main/utils_indtimesteps.F90 b/src/main/utils_indtimesteps.F90 index 47a6339d6..5fca6306f 100644 --- a/src/main/utils_indtimesteps.F90 +++ b/src/main/utils_indtimesteps.F90 @@ -463,7 +463,7 @@ subroutine print_dtlog_ind(iprint,ifrac,nfrac,time,dt,nactive,tcpu,np) real, intent(in) :: time,dt integer(kind=8), intent(in) :: nactive real(kind=4), intent(in) :: tcpu - integer, intent(in) :: np + integer(kind=8), intent(in) :: np character(len=120) :: string character(len=14) :: tmp integer, save :: nplast = 0 From c1d900911bdec6730c52a809f62953957ae15cce Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 13:14:03 +1100 Subject: [PATCH 18/21] use global npartoftypetot in test_growth --- src/tests/test_growth.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/tests/test_growth.F90 b/src/tests/test_growth.F90 index 7dc2b1097..c2e3747f3 100644 --- a/src/tests/test_growth.F90 +++ b/src/tests/test_growth.F90 @@ -107,7 +107,8 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) fxyzu,fext,Bevol,dBevol,dustprop,ddustprop,& dustfrac,dustevol,ddustevol,iphase,maxtypes,& VrelVf,dustgasprop,Omega_k,alphaind,iamtype,& - ndustlarge,ndustsmall,rhoh,deltav,this_is_a_test,periodic + ndustlarge,ndustsmall,rhoh,deltav,this_is_a_test,periodic, & + npartoftypetot,update_npartoftypetot use step_lf_global, only:step,init_step use deriv, only:get_derivs_global use energies, only:compute_energies @@ -130,8 +131,6 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) integer, intent(inout) :: ntests,npass logical, intent(in) :: frag,onefluid - integer(kind=8) :: npartoftypetot(maxtypes) - integer :: nx,nerror,nwarn integer :: itype,npart_previous,i,j,nsteps,modu,noutputs integer :: ncheck(4),nerr(4) @@ -242,7 +241,7 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) call set_particle_type(i,itype) enddo npartoftype(itype) = npart - npart_previous - npartoftypetot(itype) = reduceall_mpi('+',npartoftype(itype)) + call update_npartoftypetot massoftype(itype) = totmass/npartoftypetot(itype) !- setting dust particles if not one fluid From e1f66417a17ea40b968248e17d96f28cd00afa5d Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 13:30:22 +1100 Subject: [PATCH 19/21] call update_npartoftypetot in test_nonidealmhd after setting particle types --- src/tests/test_nonidealmhd.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/tests/test_nonidealmhd.F90 b/src/tests/test_nonidealmhd.F90 index 8f9c4f500..74df8cff0 100644 --- a/src/tests/test_nonidealmhd.F90 +++ b/src/tests/test_nonidealmhd.F90 @@ -377,6 +377,12 @@ subroutine test_standingshock(ntests,npass) Bevol(1:3,i) = leftstate(6:8)/leftstate(1) endif enddo + + ! + ! reduce types across MPI tasks + ! + call update_npartoftypetot + ! ! initialise runtime parameters ! From 02d78f256c539f666ac80ed2a5091d37e593456a Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 13:32:00 +1100 Subject: [PATCH 20/21] [indent-bot] standardised indentation --- src/main/readwrite_dumps_hdf5.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/readwrite_dumps_hdf5.F90 b/src/main/readwrite_dumps_hdf5.F90 index e55500060..1334ab1b5 100644 --- a/src/main/readwrite_dumps_hdf5.F90 +++ b/src/main/readwrite_dumps_hdf5.F90 @@ -620,7 +620,7 @@ subroutine read_any_dump_hdf5( #ifdef INJECT_PARTICLES call allocate_memory(maxp_hard) #else - call allocate_memory(int(min(nprocs,4)*nparttot/nprocs)) + call allocate_memory(int(min(nprocs,4)*nparttot/nprocs)) #endif if (periodic) then From 2d737cf39d3eb88f72be9aad1445eaf74f7e586a Mon Sep 17 00:00:00 2001 From: Conrad Chan Date: Mon, 21 Feb 2022 16:37:03 +1100 Subject: [PATCH 21/21] fix missing import --- src/tests/test_nonidealmhd.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/test_nonidealmhd.F90 b/src/tests/test_nonidealmhd.F90 index 74df8cff0..07c906fa4 100644 --- a/src/tests/test_nonidealmhd.F90 +++ b/src/tests/test_nonidealmhd.F90 @@ -277,7 +277,7 @@ subroutine test_standingshock(ntests,npass) use kernel, only:hfact_default,radkern use part, only:init_part,npart,xyzh,vxyzu,npartoftype,massoftype,set_particle_type,hrho,rhoh,& Bevol,fext,igas,iboundary,set_boundaries_to_active,alphaind,maxalpha,maxp,iphase,Bxyz,& - iamtype,iamboundary + iamtype,iamboundary,update_npartoftypetot use step_lf_global, only:step,init_step use deriv, only:get_derivs_global use testutils, only:checkval