if atm is requested on an interface then enforce it for both states

src/GRHD/callbacks.jl

@@ -71,6 +71,11 @@ function compute_atmosphere_mask(env, P, mesh::Mesh1d{SpectralElement})
     broadcast_volume!(cons2prim_doublecartoon, P.equation, env.mesh)
   elseif formulation(P) === :cartoon
     broadcast_volume!(cons2prim_cartoon, P.equation, env.mesh)
+    broadcast_volume!(determine_atmosphere, P.equation, env.mesh)
+    @unpack c2p_reset_atm = get_static_variables(env.mesh)
+    @unpack bdry_c2p_reset_atm = get_bdry_variables(env.mesh)
+    dg1d.interpolate_face_data!(mesh, c2p_reset_atm, bdry_c2p_reset_atm)
+    broadcast_faces!(equalize_atmosphere_spherical1d, P.equation, env.mesh)
   else
     TODO()
   end

src/GRHD/rhs.jl

@@ -1711,13 +1711,15 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
   @unpack flr_D, flr_Sr, flr_τ,
           ldg_D, ldg_Sr, ldg_τ,
           max_v, vr, p, smoothed_mu,
-          src_D, src_Sr, src_τ                = get_static_variables(cache)
+          src_D, src_Sr, src_τ,
+          c2p_reset_atm                       = get_static_variables(cache)
   @unpack nflr_D, nflr_Sr, nflr_τ,
           bdry_D, bdry_Sr, bdry_τ,
           bdry_max_v, bdry_vr, bdry_p,
           bdry_ldg_D, bdry_ldg_Sr, bdry_ldg_τ,
           bdry_rhs_D, bdry_rhs_Sr, bdry_rhs_τ,
-          bdry_smoothed_mu                    = get_bdry_variables(cache)
+          bdry_smoothed_mu,
+          bdry_c2p_reset_atm,                 = get_bdry_variables(cache)
 
   if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
     broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
@@ -1732,6 +1734,10 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
   dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
   dg1d.interpolate_face_data!(mesh, vr,    bdry_vr)
   dg1d.interpolate_face_data!(mesh, p,     bdry_p)
+  dg1d.interpolate_face_data!(mesh, c2p_reset_atm, bdry_c2p_reset_atm)
+
+  broadcast_volume!(determine_atmosphere, P.equation, mesh)
+  broadcast_faces!(equalize_atmosphere_spherical1d, P.equation, mesh)
 
   broadcast_volume!(flux_source_spherical1d, equation, mesh)
   impose_symmetry_sources!(P, mesh)
@@ -1967,7 +1973,6 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::AbstractArtificial
   end
 
   if !P.prms.atm_evolve
-    broadcast_volume!(determine_atmosphere, P.equation, mesh)
     broadcast_volume!(stop_atmosphere_evolution_doublecartoon, P.equation, mesh)
   end
 

src/GRHD/setup.jl

@@ -479,7 +479,8 @@ end
 function register_analysis!(mesh::Mesh1d, P)
   register_variables!(mesh,
       static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
-                              :c2p_freeze_atm, :c2p_init_admissible, :v)
+                              :c2p_freeze_atm, :c2p_init_admissible, :v),
+      bdry_variablenames = (:bdry_c2p_reset_atm,)
   )
   if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
     register_variables!(mesh, global_variablenames = (:Mtot,))
@@ -489,7 +490,8 @@ function register_analysis!(mesh::Mesh2d, P)
   register_variables!(mesh,
       static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
                               :c2p_freeze_atm, :c2p_init_admissible, :v,
-                              :buf_c2p_reset_atm,)
+                              :buf_c2p_reset_atm,),
+      bdry_variablenames = (:bdry_c2p_reset_atm,)
   )
   if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
     TODO()

src/GRHD/spherical1d.jl

@@ -439,3 +439,35 @@ end
   end
   @returns rhs_D, rhs_Sr, rhs_τ
 end
+
+
+@with_signature function equalize_atmosphere_spherical1d(eq::Equation)
+  @accepts D, Sr, τ, ρ, vr, v, ϵ, p, ρhW2, c2p_reset_atm, r
+  @accepts [bdry] bdry_c2p_reset_atm
+  c2p_reset_atm = Float64(c2p_reset_atm > 0 || bdry_c2p_reset_atm > 0)
+  bdry_c2p_reset_atm = c2p_reset_atm
+  if c2p_reset_atm > 0
+    @unpack atm_density, cold_eos = eq
+    ρatm = atm_density
+    patm = cold_eos(Pressure, Density, ρatm)
+    ϵatm = cold_eos(InternalEnergy, Density, ρatm)
+    ρ  = ρatm
+    vr = 0.0
+    v  = 0.0
+    W  = 1.0
+    p  = patm
+    ϵ  = ϵatm
+    c2p_reset_atm = 1.0
+    h  = 1.0+ϵ+p/ρ
+    W2 = 1.0/(1.0-v^2)
+    W  = sqrt(W2)
+    ρhW2 = ρ*h*W2
+    # recomputing conservatives in case we limited something
+    D  = W*ρ
+    Sr = 0.0
+    τ  = ρhW2 - p - D
+  end
+  # @show r, c2p_reset_atm, bdry_c2p_reset_atm
+  @returns D, Sr, τ, ρ, vr, v, ϵ, p, ρhW2, c2p_reset_atm
+  @returns [bdry] bdry_c2p_reset_atm
+end