Step Restarts in AstroBear

In AstroBEAR there are a variety of conditions that can arise that require the code to resolve the equations with a smaller time step. When the hydrodynamic time step is too large, the code must restart the entire root step in order to keep all of the levels of the AMR hierarchy synchronized. There are several places in the code where a restart might be requested. This is done by setting the variable lRequestRestart to true - and printing out an error message detailing the cause of this request. Then at the next convenient global communication, every processor checks to see if any other processors have requested a restart - and if so, then lRestartStep is set to true. This then triggers the exiting of the current AMR step, and the restoration of the AMR hierarchy, (and Particle positions etc..) followed by a reduction in time step.

Places where lRequestRestart is reduced

• distribution_control.f90::DistributeChildrens - (called at the start of each level step after determining where to place new children)
• elliptic_control.f90::EllipticAdvance - (prevents invoking Hypre with potentially bad data)
• time_step.f90::TestBadCFL - (called after each root step)

Places where lRequestRestart is set

• time_step.f90::PrintAdvance - (called after each level step - checks maximum hydro speed and elliptic speed on that level, as well as the max particle speed.
• elliptic_declarations.f90::CheckErr - any generic hypre error will trigger a restart request
• elliptic_declarations.f90:!CheckIters - If lRestartOnHypreMaxIters is true, this will trigger a restart if the number of iterations for the solver exceeds the maximum set
• riemann_solvers.f90::CalcFlux - if any of the fluxes have Nan's. Also prints out left and right states and coordinates
• sweep_scheme.f90::Reconstruct - If calc_eigens fails (only for interp orders 2 and 3) - although on closer inspection - this doesn't actually do anything… Probably code to check eigen system was removed for speed?
• EOS.f90::Protect_Info
  • If there are any Nan's in q
  • If density protections and lRestartOnDensityProtections is true
  • If pressure protections and lRestartOnPressureProtections is true
• source_control.f90::ExplicitSrc -If source step results in a Nan in q
• source_control.f90::SourceCell - if rk45 exceeds max iterations or if the timestep becomes insignificant

Places that call Protect_Info

• poisson.f90::SetBoxValues - protection to make sure that good values for density are passed to poisson solver
• radtransfer.f90::RadSetVectorValues - protection to make sure that good values are sent to radiation diffusion solver
• module_control.f90::BeforeStep, AfterStep, CheckData called from scrambler.f90 at beginning of simulations

In addition, there are several simple protections that occur during the hydro update at the beginning, after reconstruction and after the CTU update

Places that protect without triggering restarts in hydro code

• init_prims
• reconstruct
• CTU

  IF (q(1) < MinDensity) THEN
         q(1)=MinDensity
         q(m_low:m_high)=0d0
         IF (iE /= 0) CALL PrimSetTemp(q, MinTemp)
  ELSE IF (iE /= 0) THEN
         iE_old=q(iE)
         CALL PrimSetTemp(q, MinTemp)
         q(iE)=MAX(q(iE),iE_old)
  END IF

How max speeds are calculated

Elliptic max speed depends on the solver

• Poisson

  elliptic_maxspeed(level)=max(elliptic_maxspeed(level),sqrt(nDim*&
                  maxval(abs(Info%q(ir(1,1):ir(1,2),ir(2,1):ir(2,2),ir(3,1):ir(3,2),iPhiGas)-&
                  Info%q(il(1,1):il(1,2),il(2,1):il(2,2),il(3,1):il(3,2),iPhiGas)))))
  

Hydro max speed

• MUSCL uses maximum speed from the riemann solvers - which is usually just the max of the left and right most waves

    max(ABS(SL),ABS(SR))
  

• sweep uses the GetMaxSpeed routine to find the fastest eigen speed on the grid after an update

     SUBROUTINE GetMaxSpeed(Info)
        TYPE(InfoDef)  :: Info
        REAL(KIND=qPrec), DIMENSION(:,:,:,:), POINTER :: q
        REAL(KIND=qPREC) :: myMaxSpeed, fast_speed(3), A2, B2, B, dx, A
        INTEGER :: i,j,k,level
        level=Info%level
        dx=levels(level)%dx
        q=>Info%q
        myMaxSpeed=0d0
        DO i=1, Info%mX(1)
           DO j=1, Info%mX(2)
              DO k=1, Info%mX(3)
                 A=SoundSpeed(q(i,j,k,:))
                 A2=q(i,j,k,1)*A**2
                 IF (iE_rad /= 0) THEN
                    A2=A2+4d0/9d0*q(i,j,k,iE_rad)*(1d0-exp(-q(i,j,k,iKappa)*dx))
                    A=sqrt(A2/q(i,j,k,1))
                 END IF
                 IF (lMHD) THEN
                    B2 = SUM(q(i,j,k,iBx:iBz)**2)
                    B = A2+B2
                    fast_speed(1:nDim) = sqrt(half*(B+sqrt(B**2-4d0*A2*q(i,j,k,iB(1:nDim))**2))/q(i,j,k,1))
                    mymaxspeed= maxval(abs(q(i,j,k,imom(1:nDim)))/q(i,j,k,1)+fast_speed(1:nDim))
                 ELSE
  !                  mymaxspeed=maxval(abs(q(i,j,k,imom(1:nDim)))+A2)/q(i,j,k,1)                                                                                
                    mymaxspeed=sqrt(sum(q(i,j,k,imom(1:nDim))**2))/q(i,j,k,1)+A
                 END IF
                 IF (mymaxspeed > maxspeed(level)) THEN
                    maxspeed(level)=mymaxspeed
                    maxspeed_location(level,:)=CellPos(Info, i, j, k)
                 END IF
              END DO
           END DO
        END DO
     END SUBROUTINE GetMaxSpeed
  

Particle max speed

This just finds the fastest particle speed

ParticleMaxSpeed=max(ParticleMaxSpeed, sqrt(sum(vel(i,1:nDim)**2)))