Merge pull request #537 from HughyG/divb

removing hdivB timestep constraint as did not work

src/main/checkconserved.f90

@@ -95,7 +95,6 @@ subroutine check_conservation_error(val,ref,tol,label,decrease)
  character(len=*), intent(in) :: label
  logical, intent(in), optional :: decrease
  real :: err
- character(len=20) :: string
 
  if (abs(ref) > 1.e-3) then
     err = (val - ref)/abs(ref)
@@ -113,12 +112,7 @@ subroutine check_conservation_error(val,ref,tol,label,decrease)
        call error('evolve',trim(label)//' is not being conserved due to corotating frame',var='err',val=err)
     else
        call error('evolve','Large error in '//trim(label)//' conservation ',var='err',val=err)
-       call get_environment_variable('I_WILL_NOT_PUBLISH_CRAP',string)
-       if (.not. (trim(string)=='yes')) then
-          print "(2(/,a))",' You can ignore this error and continue by setting the ',&
-                           ' environment variable I_WILL_NOT_PUBLISH_CRAP=yes to continue'
-          call fatal('evolve',' Conservation errors too large to continue simulation')
-       endif
+       call do_not_publish_crap('evolve','Conservation errors too large to continue simulation')
     endif
  else
     if (iverbose >= 2) print "(a,es10.3)",trim(label)//' error is ',err
@@ -133,24 +127,31 @@ end subroutine check_conservation_error
 !  so is related to the checks performed here
 !+
 !----------------------------------------------------------------
-subroutine check_magnetic_stability(hdivBB_xa)
- use options, only:hdivbbmax_max
+subroutine check_magnetic_stability(hdivBonB_ave,hdivBonB_max)
  use io,      only:fatal
- real, intent(in) :: hdivBB_xa(:)
+ real, intent(in) :: hdivBonB_ave,hdivBonB_max
 
- if (hdivbbmax_max < 1.1) then
-    ! In this regime, we assume the user has not modified this value,
-    ! either by choice or by being unaware of this.  This warning will
-    ! appear in this case.
-    if (hdivBB_xa(1) > 100 .or. hdivBB_xa(2) > 0.1) then
-       ! Tricco, Price & Bate (2016) suggest the average should remain lower than 0.01,
-       ! but we will increase it here due to the nature of the exiting the code
-       ! The suggestion of 512 was empirically determined in Dobbs & Wurster (2021)
-       call fatal('evolve','h|divb|/b is too large; recommend hdivbbmax_max = 512; set >1.2 to suppress this message.')
-    endif
+ if (hdivBonB_max > 100 .or. hdivBonB_ave > 0.1) then
+   ! Tricco, Price & Bate (2016) suggest the average should remain lower than 0.01,
+   ! but we will increase it here due to the nature of the exiting the code
+   ! The suggestion of 512 was empirically determined in Dobbs & Wurster (2021)
+   call do_not_publish_crap('evolve','h|divb|/b is too large; recommend to increase the overcleanfac')
  endif
 
 end subroutine check_magnetic_stability
 
+subroutine do_not_publish_crap(subr,msg)
+ use io, only:fatal
+ character(len=*), intent(in) :: subr,msg
+ character(len=20) :: string
+
+ call get_environment_variable('I_WILL_NOT_PUBLISH_CRAP',string)
+ if (.not. (trim(string)=='yes')) then
+    print "(2(/,a))",' You can ignore this error and continue by setting the ',&
+                     ' environment variable I_WILL_NOT_PUBLISH_CRAP=yes to continue'
+    call fatal(subr,msg)
+ endif
+
+end subroutine do_not_publish_crap
 !----------------------------------------------------------------
 end module checkconserved

src/main/energies.F90

@@ -26,7 +26,8 @@ module energies
  implicit none
 
  logical,         public    :: gas_only,track_mass,track_lum
- real,            public    :: ekin,etherm,emag,epot,etot,totmom,angtot,mtot,xyzcom(3),hdivBB_xa(2)
+ real,            public    :: ekin,etherm,emag,epot,etot,totmom,angtot,mtot,xyzcom(3)
+ real,            public    :: hdivBonB_ave,hdivBonB_max
  real,            public    :: vrms,rmsmach,accretedmass,mdust(maxdusttypes),mgas
  real,            public    :: xmom,ymom,zmom
  real,            public    :: totlum
@@ -730,8 +731,8 @@ subroutine compute_energies(t)
  endif
 
  if (mhd) then
-    hdivBB_xa(1) = ev_data(iev_max,iev_hdivB)
-    hdivBB_xa(2) = ev_data(iev_ave,iev_hdivB)
+    hdivBonB_max = ev_data(iev_max,iev_hdivB)
+    hdivBonB_ave = ev_data(iev_ave,iev_hdivB)
  endif
 
  if (maxp==maxp_hard) then

src/main/evolve.F90

@@ -37,7 +37,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
                             dtmax_ifactor,dtmax_ifactorWT,dtmax_dratio,check_dtmax_for_decrease,&
                             idtmax_n,idtmax_frac,idtmax_n_next,idtmax_frac_next
  use evwrite,          only:write_evfile,write_evlog
- use energies,         only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0,hdivBB_xa
+ use energies,         only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0,hdivBonB_ave,hdivBonB_max
  use checkconserved,   only:etot_in,angtot_in,totmom_in,mdust_in,&
                             init_conservation_checks,check_conservation_error,&
                             check_magnetic_stability
@@ -396,7 +396,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
              call check_conservation_error(mdust(j),mdust_in(j),1.e-1,'dust mass',decrease=.true.)
           enddo
        endif
-       if (mhd) call check_magnetic_stability(hdivBB_xa)
+       if (mhd) call check_magnetic_stability(hdivBonB_ave,hdivBonB_max)
        if (id==master) then
           if (np_e_eq_0  > 0) call warning('evolve','N gas particles with energy = 0',var='N',ival=int(np_e_eq_0,kind=4))
           if (np_cs_eq_0 > 0) call warning('evolve','N gas particles with sound speed = 0',var='N',ival=int(np_cs_eq_0,kind=4))

src/main/force.F90

@@ -150,27 +150,26 @@ module forces
        idBevolyi      = 10, &
        idBevolzi      = 11, &
        idivBdiffi     = 12, &
-       ihdivBBmax     = 13, &
  !--dust array indexing
-       ifdragxi       = 14, &
-       ifdragyi       = 15, &
-       ifdragzi       = 16, &
-       iddustevoli    = 17, &
-       iddustevoliend = 17 +   (maxdustsmall-1), &
-       idudtdusti     = 18 +   (maxdustsmall-1), &
-       idudtdustiend  = 18 + 2*(maxdustsmall-1), &
-       ideltavxi      = 19 + 2*(maxdustsmall-1), &
-       ideltavxiend   = 19 + 3*(maxdustsmall-1), &
-       ideltavyi      = 20 + 3*(maxdustsmall-1), &
-       ideltavyiend   = 20 + 4*(maxdustsmall-1), &
-       ideltavzi      = 21 + 4*(maxdustsmall-1), &
-       ideltavziend   = 21 + 5*(maxdustsmall-1), &
-       idvix          = 22 + 5*(maxdustsmall-1), &
-       idviy          = 23 + 5*(maxdustsmall-1), &
-       idviz          = 24 + 5*(maxdustsmall-1), &
-       idensgasi      = 25 + 5*(maxdustsmall-1), &
-       icsi           = 26 + 5*(maxdustsmall-1), &
-       idradi         = 26 + 5*(maxdustsmall-1) + 1
+       ifdragxi       = 13, &
+       ifdragyi       = 14, &
+       ifdragzi       = 15, &
+       iddustevoli    = 16, &
+       iddustevoliend = 16 +   (maxdustsmall-1), &
+       idudtdusti     = 17 +   (maxdustsmall-1), &
+       idudtdustiend  = 17 + 2*(maxdustsmall-1), &
+       ideltavxi      = 18 + 2*(maxdustsmall-1), &
+       ideltavxiend   = 18 + 3*(maxdustsmall-1), &
+       ideltavyi      = 19 + 3*(maxdustsmall-1), &
+       ideltavyiend   = 19 + 4*(maxdustsmall-1), &
+       ideltavzi      = 20 + 4*(maxdustsmall-1), &
+       ideltavziend   = 20 + 5*(maxdustsmall-1), &
+       idvix          = 21 + 5*(maxdustsmall-1), &
+       idviy          = 22 + 5*(maxdustsmall-1), &
+       idviz          = 23 + 5*(maxdustsmall-1), &
+       idensgasi      = 24 + 5*(maxdustsmall-1), &
+       icsi           = 25 + 5*(maxdustsmall-1), &
+       idradi         = 25 + 5*(maxdustsmall-1) + 1
 
  private
 
@@ -1621,7 +1620,6 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g
              else
                 Bj1 = 0.0
              endif
-             fsum(ihdivBBmax) = max( hj*abs(divcurlB(1,j))*Bj1, fsum(ihdivBBmax))
              !
              ! non-ideal MHD terms
              !
@@ -2522,7 +2520,7 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv
  use dim,            only:mhd,mhd_nonideal,lightcurve,use_dust,maxdvdx,use_dustgrowth,gr,use_krome,&
                           store_dust_temperature,do_nucleation,update_muGamma,h2chemistry
  use eos,            only:ieos,iopacity_type
- use options,        only:alpha,ipdv_heating,ishock_heating,psidecayfac,overcleanfac,hdivbbmax_max, &
+ use options,        only:alpha,ipdv_heating,ishock_heating,psidecayfac,overcleanfac, &
                           use_dustfrac,damp,icooling,implicit_radiation
  use part,           only:rhoanddhdrho,iboundary,igas,maxphase,maxvxyzu,nptmass,xyzmh_ptmass,eos_vars, &
                           massoftype,get_partinfo,tstop,strain_from_dvdx,ithick,iradP,sinks_have_heating,&
@@ -2598,7 +2596,7 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv
  real    :: Bxyzi(3),curlBi(3),dvdxi(9),straini(6)
  real    :: xi,yi,zi,B2i,f2i,divBsymmi,betai,frac_divB,divBi,vcleani
  real    :: pri,spsoundi,drhodti,divvi,shearvisc,fac,pdv_work
- real    :: psii,dtau,hdivbbmax
+ real    :: psii,dtau
  real    :: eni,dudtnonideal
  real    :: dustfraci(maxdusttypes),dustfracisum
  real    :: tstopi(maxdusttypes),tseff,dtdustdenom
@@ -2965,21 +2963,7 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv
              ! new cleaning evolving d/dt (psi/c_h)
              dBevol(4,i) = -vcleani*fsum(idivBdiffi)*rho1i - psii*dtau - 0.5*psii*divvi
 
-             ! timestep from cleaning
-             !   1. the factor of 10 in hdivbbmax is empirical from checking how much
-             !      spurious B-fields are decreased in colliding flows
-             !   2. if overcleaning is on (i.e. hdivbbmax > 1.0), then factor of 2 is
-             !      from empirical tests to ensure that overcleaning with individual
-             !      timesteps is stable
-             if (B2i > 0.) then
-                hdivbbmax = hi*abs(divBi)/sqrt(B2i)
-             else
-                hdivbbmax = 0.0
-             endif
-             hdivbbmax = max( overcleanfac, 10.*hdivbbmax, 10.*fsum(ihdivBBmax) )
-             hdivbbmax = min( hdivbbmax, hdivbbmax_max )
-             if (hdivbbmax > 1.0) hdivbbmax = 2.0*hdivbbmax
-             dtclean   = C_cour*hi/(hdivbbmax * vwavei + tiny(0.))
+             dtclean   = C_cour*hi/(vcleani + tiny(0.))
           endif
        endif
 

src/main/options.f90

@@ -36,7 +36,7 @@ module options
 
  real, public :: alpha,alphau,beta
  real, public :: alphamax
- real, public :: alphaB, psidecayfac, overcleanfac, hdivbbmax_max
+ real, public :: alphaB, psidecayfac, overcleanfac
  integer, public :: ishock_heating,ipdv_heating,icooling,iresistive_heating
  integer, public :: ireconav
 
@@ -134,8 +134,6 @@ subroutine set_default_options
  alphaB            = 1.0
  psidecayfac       = 1.0     ! psi decay factor (MHD only)
  overcleanfac      = 1.0     ! factor to increase signal velocity for (only) time steps and psi cleaning
- hdivbbmax_max     = 1.0     ! if > overcleanfac, then use B/(h*|div B|) as a coefficient for dtclean;
- !                           ! this is the max value allowed; test suggest =512 for magnetised colliding flows
  beta              = 2.0     ! beta viscosity term
  avdecayconst      = 0.1     ! decay time constant for viscosity switches
 

src/main/readwrite_infile.F90

@@ -34,7 +34,6 @@ module readwrite_infile
 !   - dtwallmax          : *maximum wall time between dumps (hhh:mm, 000:00=ignore)*
 !   - dumpfile           : *dump file to start from*
 !   - flux_limiter       : *limit radiation flux*
-!   - hdivbbmax_max      : *max factor to decrease cleaning timestep propto B/(h|divB|)*
 !   - hfact              : *h in units of particle spacing [h = hfact(m/rho)^(1/3)]*
 !   - ien_type           : *energy variable (0=auto, 1=entropy, 2=energy, 3=entropy_s)*
 !   - implicit_radiation : *use implicit integration (Whitehouse, Bate & Monaghan 2005)*
@@ -76,7 +75,7 @@ module readwrite_infile
  use options,   only:nfulldump,nmaxdumps,twallmax,iexternalforce,tolh, &
                      alpha,alphau,alphaB,beta,avdecayconst,damp,rkill, &
                      ipdv_heating,ishock_heating,iresistive_heating,ireconav, &
-                     icooling,psidecayfac,overcleanfac,hdivbbmax_max,alphamax,calc_erot,rhofinal_cgs, &
+                     icooling,psidecayfac,overcleanfac,alphamax,calc_erot,rhofinal_cgs, &
                      use_mcfost,use_Voronoi_limits_file,Voronoi_limits_file,use_mcfost_stellar_parameters,&
                      exchange_radiation_energy,limit_radiation_flux,iopacity_type,mcfost_computes_Lacc,&
                      mcfost_uses_PdV,implicit_radiation,mcfost_keep_part,ISM, mcfost_dust_subl
@@ -202,7 +201,6 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint)
     call write_inopt(alphaB,'alphaB','shock resistivity parameter',iwritein)
     call write_inopt(psidecayfac,'psidecayfac','div B diffusion parameter',iwritein)
     call write_inopt(overcleanfac,'overcleanfac','factor to increase cleaning speed (decreases time step)',iwritein)
-    call write_inopt(hdivbbmax_max,'hdivbbmax_max','max factor to decrease cleaning timestep propto B/(h|divB|)',iwritein)
  endif
  call write_inopt(beta,'beta','beta viscosity',iwritein)
  if (maxalpha==maxp .and. maxp > 0) then
@@ -481,8 +479,6 @@ subroutine read_infile(infile,logfile,evfile,dumpfile)
        read(valstring,*,iostat=ierr) psidecayfac
     case('overcleanfac')
        read(valstring,*,iostat=ierr) overcleanfac
-    case('hdivbbmax_max')
-       read(valstring,*,iostat=ierr) hdivbbmax_max
     case('beta')
        read(valstring,*,iostat=ierr) beta
     case('ireconav')
@@ -684,10 +680,6 @@ subroutine read_infile(infile,logfile,evfile,dumpfile)
        if (psidecayfac < 0.) call fatal(label,'stupid value for psidecayfac')
        if (psidecayfac > 2.) call warn(label,'psidecayfac set outside recommended range (0.1-2.0)')
        if (overcleanfac < 1.0) call warn(label,'overcleanfac less than 1')
-       if (hdivbbmax_max < overcleanfac) then
-          call warn(label,'Resetting hdivbbmax_max = overcleanfac')
-          hdivbbmax_max = overcleanfac
-       endif
     endif
     if (beta < 0.)     call fatal(label,'beta < 0')
     if (beta > 4.)     call warn(label,'very high beta viscosity set')

src/setup/setup_sphereinbox.f90

@@ -110,7 +110,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact_
  use dust,         only:ilimitdustflux
  use timestep,     only:dtmax,tmax,dtmax_dratio,dtmax_min
  use centreofmass, only:reset_centreofmass
- use options,      only:nfulldump,rhofinal_cgs,hdivbbmax_max,use_dustfrac,icooling
+ use options,      only:nfulldump,rhofinal_cgs,use_dustfrac,icooling
  use kernel,       only:hfact_default
  use mpidomain,    only:i_belong
  use ptmass,       only:icreate_sinks,r_crit,h_acc,h_soft_sinksink
@@ -664,7 +664,6 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact_
     icreate_sinks = icreate_sinks_setup
     r_crit        = r_crit_setup
     h_acc         = h_acc_setup
-    hdivbbmax_max = 1. !512.
     ! reset defaults based upon options
     if (density_contrast > 1.) dtmax_dratio = 1.258
     if (density_contrast < 1.+epsilon(density_contrast) .and. maxvxyzu>=4) then