(mpi) bug fix/disable test_wind with mpi, remove ifdefs (#55)

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

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,28 +332,28 @@ 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
 !
 !--sanity checks on magnetic field
 !
  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

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
 
 !----------------------------------------------------------------

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

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'