diff --git a/src/main/checkoptions.F90 b/src/main/checkoptions.F90 index 6dfa483e8..2102df2aa 100644 --- a/src/main/checkoptions.F90 +++ b/src/main/checkoptions.F90 @@ -31,7 +31,7 @@ module checkoptions !------------------------------------------------------------------- subroutine check_compile_time_settings(ierr) use part, only:mhd,gravity,ngradh,h2chemistry,maxvxyzu,use_dust,gr - use dim, only:use_dustgrowth,maxtypes + use dim, only:use_dustgrowth,maxtypes,mpi,inject_parts use io, only:error,id,master,fatal,warning #ifdef GR use metric_tools, only:icoordinate,icoord_cartesian @@ -142,10 +142,14 @@ subroutine check_compile_time_settings(ierr) #endif #ifdef DUSTGROWTH - if (.not. use_dustgrowth) call error(string,'-DDUSTGROWTH but use_dustgrowth = .false.') + if (.not. use_dustgrowth) then + call error(string,'-DDUSTGROWTH but use_dustgrowth = .false.') + ierr = 16 + endif #endif - return + if (mpi .and. inject_parts) call error(string,'MPI currently not compatible with particle injection') + end subroutine check_compile_time_settings end module checkoptions diff --git a/src/main/checksetup.F90 b/src/main/checksetup.F90 index 2ea04643b..6251c415b 100644 --- a/src/main/checksetup.F90 +++ b/src/main/checksetup.F90 @@ -293,7 +293,8 @@ subroutine check_setup(nerror,nwarn,restart) ! warn about external force settings ! if (iexternalforce==iext_star .and. nptmass==0) then - print*,'WARNING: iexternalforce=1 does not conserve momentum - use a sink particle at r=0 if you care about this' + if (id==master) print "(a,/,a)",'WARNING: iexternalforce=1 does not conserve momentum:',& + ' use a sink particle at r=0 if you care about this' nwarn = nwarn + 1 endif ! @@ -331,15 +332,15 @@ subroutine check_setup(nerror,nwarn,restart) if (gravity .or. nptmass > 0) then if (.not.G_is_unity()) then if (gravity) then - print*,'ERROR: self-gravity ON but G /= 1 in code units, got G=',get_G_code() + if (id==master) print*,'ERROR: self-gravity ON but G /= 1 in code units, got G=',get_G_code() elseif (nptmass > 0) then - print*,'ERROR: sink particles used but G /= 1 in code units, got G=',get_G_code() + if (id==master) print*,'ERROR: sink particles used but G /= 1 in code units, got G=',get_G_code() endif nerror = nerror + 1 endif endif if (.not. gr .and. (gravity .or. mhd) .and. ien_type == ien_etotal) then - print*,'Cannot use total energy with self gravity or mhd' + if (id==master) print*,'Cannot use total energy with self gravity or mhd' nerror = nerror + 1 endif ! @@ -347,12 +348,12 @@ subroutine check_setup(nerror,nwarn,restart) ! if (mhd) then if (all(abs(Bxyz(:,1:npart)) < tiny(0.))) then - print*,'WARNING: MHD is ON but magnetic field is zero everywhere' + if (id==master) print*,'WARNING: MHD is ON but magnetic field is zero everywhere' nwarn = nwarn + 1 endif if (mhd_nonideal) then if (n_nden /= n_nden_phantom) then - print*,'ERROR: n_nden in nicil.f90 needs to match n_nden_phantom in config.F90; n_nden = ',n_nden + if (id==master) print*,'ERROR: n_nden in nicil.f90 needs to match n_nden_phantom in config.F90; n_nden = ',n_nden nerror = nerror + 1 endif endif diff --git a/src/main/mpi_balance.F90 b/src/main/mpi_balance.F90 index a73fbbfd2..91e604027 100644 --- a/src/main/mpi_balance.F90 +++ b/src/main/mpi_balance.F90 @@ -44,19 +44,33 @@ module mpibalance contains -subroutine allocate_balance_arrays +!---------------------------------------------------------------- +!+ +! allocate memory +!+ +!---------------------------------------------------------------- +subroutine allocate_balance_arrays() use allocutils, only:allocate_array + call allocate_array('nsent', nsent, nprocs) call allocate_array('nexpect', nexpect, nprocs) call allocate_array('nrecv', nrecv, nprocs) call allocate_array('countrequest',countrequest,nprocs) + end subroutine allocate_balance_arrays -subroutine deallocate_balance_arrays +!---------------------------------------------------------------- +!+ +! deallocate memory +!+ +!---------------------------------------------------------------- +subroutine deallocate_balance_arrays() + if (allocated(nsent )) deallocate(nsent ) if (allocated(nexpect )) deallocate(nexpect ) if (allocated(nrecv )) deallocate(nrecv ) if (allocated(countrequest)) deallocate(countrequest) + end subroutine deallocate_balance_arrays !---------------------------------------------------------------- diff --git a/src/main/partinject.F90 b/src/main/partinject.F90 index ae06241b4..5b4edb5d2 100644 --- a/src/main/partinject.F90 +++ b/src/main/partinject.F90 @@ -21,14 +21,16 @@ module partinject ! timestep_ind ! implicit none - character(len=80), parameter, public :: & ! module version - modid="$Id$" public :: add_or_update_particle, add_or_update_sink public :: update_injected_particles + + ! ! Use this flag if particles are updated rather than injected (e.g. inject_sne) ! see inject_sne for use; currently only valid for gas particles + ! logical, public :: updated_particle = .false. + private contains @@ -42,12 +44,10 @@ subroutine add_or_update_particle(itype,position,velocity,h,u,particle_number,np use part, only:maxp,iamtype,iphase,maxvxyzu,iboundary,nucleation use part, only:maxalpha,alphaind,maxgradh,gradh,fxyzu,fext,set_particle_type use part, only:mhd,Bevol,dBevol,Bxyz,divBsymm!,dust_temp - use part, only:divcurlv,divcurlB,ndivcurlv,ndivcurlB,ntot + use part, only:divcurlv,divcurlB,ndivcurlv,ndivcurlB,ntot,ibin use io, only:fatal -#ifdef IND_TIMESTEPS - use part, only:ibin + use dim, only:ind_timesteps use timestep_ind, only:nbinmax -#endif integer, intent(in) :: itype real, intent(in) :: position(3), velocity(3), h, u real, intent(in), optional :: JKmuS(:) @@ -75,32 +75,37 @@ subroutine add_or_update_particle(itype,position,velocity,h,u,particle_number,np npartoftype(itype_old) = npartoftype(itype_old) - 1 npartoftype(itype) = npartoftype(itype)+1 endif + call set_particle_type(particle_number,itype) + ! Update particle type, position, size, velocity and energy xyzh(1,particle_number) = position(1) xyzh(2,particle_number) = position(2) xyzh(3,particle_number) = position(3) if (itype /= iboundary .or. new_particle) xyzh(4,particle_number) = h + vxyzu(1,particle_number) = velocity(1) vxyzu(2,particle_number) = velocity(2) vxyzu(3,particle_number) = velocity(3) if (maxvxyzu>=4) vxyzu(4,particle_number) = u + fxyzu(:,particle_number) = 0. fext(:,particle_number) = 0. + if (mhd) then Bevol(:,particle_number) = 0. dBevol(:,particle_number) = 0. Bxyz(:,particle_number) = 0. divBsymm(particle_number) = 0. endif + if (ndivcurlv > 0) divcurlv(:,particle_number) = 0. if (ndivcurlB > 0) divcurlB(:,particle_number) = 0. if (maxalpha==maxp) alphaind(:,particle_number) = 0. if (maxgradh==maxp) gradh(:,particle_number) = 0. !if (store_dust_temperature) dust_temp(:,particle_number) = 0. -#ifdef IND_TIMESTEPS - ibin(particle_number) = nbinmax -#endif + + if (ind_timesteps) ibin(particle_number) = nbinmax if (present(jKmuS)) nucleation(:,particle_number) = JKmuS(:) end subroutine add_or_update_particle @@ -125,6 +130,7 @@ subroutine add_or_update_sink(position,velocity,radius,mass,sink_number) elseif (sink_number > nptmass+1) then call fatal('Add sink', 'Incorrect sink number (> maxptmass + 1).') endif + ! Update sink position, size, velocity and mass xyzmh_ptmass(1,sink_number) = position(1) xyzmh_ptmass(2,sink_number) = position(2) @@ -134,6 +140,7 @@ subroutine add_or_update_sink(position,velocity,radius,mass,sink_number) vxyz_ptmass(1,sink_number) = velocity(1) vxyz_ptmass(2,sink_number) = velocity(2) vxyz_ptmass(3,sink_number) = velocity(3) + end subroutine add_or_update_sink !----------------------------------------------------------------------- @@ -144,11 +151,9 @@ end subroutine add_or_update_sink !+ !----------------------------------------------------------------------- subroutine update_injected_particles(npartold,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject) -#ifdef IND_TIMESTEPS + use dim, only:ind_timesteps use timestep_ind, only:get_newbin,change_nbinmax,get_dt - use part, only:twas,ibin,ibin_old -#endif - use part, only:norig,iorig,iphase,igas,iunknown + use part, only:twas,ibin,ibin_old,norig,iorig,iphase,igas,iunknown #ifdef GR use part, only:xyzh,vxyzu,pxyzu,dens,metrics,metricderivs,fext use cons2prim, only:prim2consall @@ -162,9 +167,7 @@ subroutine update_injected_particles(npartold,npart,istepfrac,nbinmax,time,dtmax real, intent(inout) :: dt real, intent(in) :: time,dtmax,dtinject integer :: i -#ifdef IND_TIMESTEPS integer(kind=1) :: nbinmaxprev -#endif #ifdef GR real :: dtext_dum #endif @@ -184,24 +187,24 @@ subroutine update_injected_particles(npartold,npart,istepfrac,nbinmax,time,dtmax endif #endif -#ifdef IND_TIMESTEPS - ! find timestep bin associated with dtinject - nbinmaxprev = nbinmax - call get_newbin(dtinject,dtmax,nbinmax,allow_decrease=.false.) - if (nbinmax > nbinmaxprev) then ! update number of bins if needed - call change_nbinmax(nbinmax,nbinmaxprev,istepfrac,dtmax,dt) + if (ind_timesteps) then + ! find timestep bin associated with dtinject + nbinmaxprev = nbinmax + call get_newbin(dtinject,dtmax,nbinmax,allow_decrease=.false.) + if (nbinmax > nbinmaxprev) then ! update number of bins if needed + call change_nbinmax(nbinmax,nbinmaxprev,istepfrac,dtmax,dt) + endif + ! put all injected particles on shortest bin + do i=npartold+1,npart + ibin(i) = nbinmax + ibin_old(i) = nbinmax ! for particle waking to ensure that neighbouring particles are promptly woken + twas(i) = time + 0.5*get_dt(dtmax,ibin(i)) + enddo + else + ! For global timestepping, reset the timestep, since this is otherwise + ! not updated until after the call to step. + dt = min(dt,dtinject) endif - ! put all injected particles on shortest bin - do i=npartold+1,npart - ibin(i) = nbinmax - ibin_old(i) = nbinmax ! for particle waking to ensure that neighbouring particles are promptly woken - twas(i) = time + 0.5*get_dt(dtmax,ibin(i)) - enddo -#else - ! For global timestepping, reset the timestep, since this is otherwise - ! not updated until after the call to step. - dt = min(dt,dtinject) -#endif ! add particle ID do i=npartold+1,npart @@ -214,11 +217,11 @@ subroutine update_injected_particles(npartold,npart,istepfrac,nbinmax,time,dtmax do i=1,npart if (iphase(i) == iunknown) then iphase(i) = igas -#ifdef IND_TIMESTEPS - ibin(i) = nbinmax - ibin_old(i) = nbinmax - twas(i) = time + 0.5*get_dt(dtmax,ibin(i)) -#endif + if (ind_timesteps) then + ibin(i) = nbinmax + ibin_old(i) = nbinmax + twas(i) = time + 0.5*get_dt(dtmax,ibin(i)) + endif endif enddo endif diff --git a/src/tests/test_wind.f90 b/src/tests/test_wind.f90 index d9622bc73..603e38bff 100644 --- a/src/tests/test_wind.f90 +++ b/src/tests/test_wind.f90 @@ -42,7 +42,7 @@ subroutine test_wind(ntests,npass) use timestep, only:time,tmax,dt,dtmax,nsteps,dtrad,dtforce,dtcourant,dterr,print_dtlog use step_lf_global, only:step,init_step use testutils, only:checkval,update_test_scores - use dim, only:isothermal,inject_parts + use dim, only:isothermal,inject_parts,mpi use partinject, only:update_injected_particles use timestep_ind, only:nbinmax use wind, only:trvurho_1D @@ -54,7 +54,10 @@ subroutine test_wind(ntests,npass) integer :: i,ierr,nerror,istepfrac,npart_old,nfailed(9),nwarn real :: dtinject,dtlast,t,default_particle_mass,dtext,dtnew,dtprint,dtmaxold,tprint - if (inject_type /= 'wind') then + if (mpi) then + if (id==master) write(*,"(/,a,/)") '--> SKIPPING WIND TEST (currently not working with MPI)' + return + elseif (inject_type /= 'wind') then if (id==master) write(*,"(/,a,/)") '--> SKIPPING WIND TEST (need to compile with wind injection module)' return else @@ -166,7 +169,6 @@ subroutine test_wind(ntests,npass) call checkval(xyzmh_ptmass(7,1),0.,epsilon(0.),nfailed(7),'mass accreted') call checkval(npart,12180,0,nfailed(8),'number of ejected particles') call checkval(xyzmh_ptmass(15,1),1.591640703559762E-06,epsilon(0.),nfailed(9),'wind mass loss rate') - call update_test_scores(ntests,nfailed,npass) if (id==master) write(*,"(/,a)") '<-- WIND TEST COMPLETE'