diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl index 51c57d9beda2293f3e67bbfa7f035477bbb623b2..c0265dbb940fd512191e1acae8c39feeb2507bcd 100644 --- a/src/GRHD/callbacks.jl +++ b/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 diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl index 132b6119749fa040af19616fcccccbd6fede40a0..cfb3ca24578a087d53583513df98a30678184f30 100644 --- a/src/GRHD/rhs.jl +++ b/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 diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl index 6e85e21af62693f7a56d5caf4671bbee2d30f04b..7d7235051c2924c7d550c0043e4cdb93da60ce6b 100644 --- a/src/GRHD/setup.jl +++ b/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() diff --git a/src/GRHD/spherical1d.jl b/src/GRHD/spherical1d.jl index 5a6646c4909fcd3ba3180c4c9bc1508e50743df4..ce0bfd81c01a6f52e861f86be485f74602d8a8b1 100644 --- a/src/GRHD/spherical1d.jl +++ b/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