diff --git a/src/GRHD/GRHD.jl b/src/GRHD/GRHD.jl
index f5c63495c8fd0233927c35c7f75f20d0e90f9d23..a3174a8b4e51f44b9b6597dced40553d094cdd95 100644
--- a/src/GRHD/GRHD.jl
+++ b/src/GRHD/GRHD.jl
@@ -60,6 +60,14 @@ dg1d.@parameters GRHD begin
atm_evolve = true
@check atm_evolve isa Bool
+ """
+ If `true` then all states on an interface are set to atmosphere if at least one state is atmosphere.
+
+ Note: So far only supported for `spherical1d` formulation.
+ """
+ atm_equalize_on_interface = false
+ @check atm_equalize_on_interface isa Bool
+
"""
Available options:
- `tov_reflective_origin`: for `id = tov` use reflection symmetry across coordinate origin;
@@ -100,12 +108,6 @@ dg1d.@parameters GRHD begin
"""
id_filename = ""
- """
- Restmass density for `atmosphere` initial data.
- """
- id_atmosphere_density = 1e-11
- @check id_atmosphere_density > 0.0
-
"""
Choice of formulation of evolution equations. Available options:
- `spherical1d`: generalized Valencia equations from Montero+2014 with choice
@@ -199,7 +201,8 @@ dg1d.@parameters GRHD begin
@check slope_limiter_tvb_M > 0.0
"""
- If `true` then a causal atmosphere is enforced after every full time step.
+ If `true` then atmosphere is enforced after every full time step in a way that it only has
+ causal effects on its surrounding.
"""
c2p_enforce_causal_atm = true
@check c2p_enforce_causal_atm isa Bool
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index 1b302c14e9d2e0372c5e4956d3ec7f0a8bca77ba..9192c61ca4275d44e2846bfc147fc9966657ca0e 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -71,16 +71,25 @@ 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)
+ if P.prms.equalize_atmosphere_spherical1d
+ 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)
+ end
else
TODO()
end
if P.prms.c2p_dynamic_atm
broadcast_volume!(determine_atmosphere, P.equation, env.mesh)
end
+ # this controls the effect of the stop_atmosphere_evolution_* methods
+ # which trigger when c2p_reset_atm > 0.0 && c2p_freeze_atm > 0.0 at a point
if enforce
@unpack c2p_freeze_atm = get_static_variables(env.mesh)
- c2p_freeze_atm .= 0.0
- elseif P.prms.c2p_dynamic_atm
+ c2p_freeze_atm .= 1.0 # always on
+ elseif enforce_causal
update_atm_domain_of_dependence!(env.mesh)
end
end
@@ -91,30 +100,64 @@ function update_atm_domain_of_dependence!(mesh::Mesh1d)
Nx = mesh.element.Npts
for (cidx,(frz,rst)) in enumerate(zip(eachcell(mesh,freeze),eachcell(mesh,reset)))
# we need to special case the reset check for when we hit an interface or not
- # reset check in bulk
+ # ========================================
+ # reset check in bulk (requires Nx ≥ 5)
# --x|x--x---x-----x---x--x|x--
# | | |
# r1 r2 r3
- # reset check on interface
- # ---x--x|x--x---- ---x--x|x--x----
- # | | | | or | | | |
- # r1 r2 r3 r4 r1 r2 r3 r4
cell = mesh.tree.cells[cidx]
# points in bulk
- for i in 2:Nx-1
- r1, r2, r3 = rst[i-1], rst[i], rst[i+1]
- frz[i] = Float64(r1==r2==r3)
+ if Nx ≥ 5
+ for i in 3:Nx-2
+ r1, r2, r3 = rst[i-1], rst[i], rst[i+1]
+ frz[i] = Float64(r1==r2==r3)
+ end
end
- # points on edge
+ # ========================================
+ # reset check close to interface
+ # --x|x--x---x-----x---x--x|x--
+ # | | | |
+ # r1 r2 r3 r4
+ # i = 2
if dg1d.has_neighbor(cell,Cart1d.Xmin)
nrst = dg1d.view_neighbor(reset, cell, Cart1d.Xmin, mesh)
- r1, r2, r3, r4 = rst[1], rst[2], nrst[Nx], nrst[Nx-1]
+ r1, r2, r3, r4 = nrst[Nx], rst[1], rst[2], rst[3]
+ frz[2] = Float64(r1==r2==r3==r4)
+ else
+ r2, r3, r4 = rst[1], rst[2], rst[3]
+ frz[2] = Float64(r2==r3==r4)
+ end
+ # i = Nx-1
+ if dg1d.has_neighbor(cell,Cart1d.Xmax)
+ nrst = dg1d.view_neighbor(reset, cell, Cart1d.Xmax, mesh)
+ r1, r2, r3, r4 = rst[Nx-2], rst[Nx-1], rst[Nx], nrst[1]
+ frz[Nx-1] = Float64(r1==r2==r3==r4)
+ else
+ r1, r2, r3 = rst[Nx-2], rst[Nx-1], rst[Nx]
+ frz[Nx-1] = Float64(r1==r2==r3)
+ end
+ # ========================================
+ # reset check on interface
+ # --x----x--x|x--x----x--
+ # | | | |
+ # r1 r2 r3 r4
+ # i = 1
+ if dg1d.has_neighbor(cell,Cart1d.Xmin)
+ nrst = dg1d.view_neighbor(reset, cell, Cart1d.Xmin, mesh)
+ r1, r2, r3, r4 = nrst[Nx-1], nrst[Nx], rst[1], rst[2]
frz[1] = Float64(r1==r2==r3==r4)
+ else
+ r3, r4 = rst[1], rst[2]
+ frz[1] = Float64(r3==r4)
end
+ # i = Nx
if dg1d.has_neighbor(cell,Cart1d.Xmax)
nrst = dg1d.view_neighbor(reset, cell, Cart1d.Xmax, mesh)
- r1, r2, r3, r4 = rst[Nx], rst[Nx-1], nrst[1], nrst[2]
+ r1, r2, r3, r4 = rst[Nx-1], rst[Nx], nrst[1], nrst[2]
frz[Nx] = Float64(r1==r2==r3==r4)
+ else
+ r1, r2 = rst[Nx-1], rst[Nx]
+ frz[Nx] = Float64(r1==r2)
end
end
freeze[1] = freeze[end] = 1.0
diff --git a/src/GRHD/cartoon.jl b/src/GRHD/cartoon.jl
index b81f5cac65d05d571d5763f5481ee8f4cff9f3a6..37bad364ba27291f3ab231e169110cb754010bc1 100644
--- a/src/GRHD/cartoon.jl
+++ b/src/GRHD/cartoon.jl
@@ -617,8 +617,8 @@ function impose_symmetry_cartoon!(P::Project{:cartoon}, mesh::Mesh)
end
@with_signature function stop_atmosphere_evolution_cartoon(eq::Equation)
- @accepts rhs_rD, rhs_rSx, rhs_rSz, rhs_rτ, c2p_reset_atm
- if c2p_reset_atm > 0
+ @accepts rhs_rD, rhs_rSx, rhs_rSz, rhs_rτ, c2p_reset_atm, c2p_freeze_atm
+ if c2p_reset_atm > 0 && c2p_freeze_atm > 0
rhs_rD = rhs_rSx = rhs_rSz = rhs_rτ = 0.0
end
@returns rhs_rD, rhs_rSx, rhs_rSz, rhs_rτ
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index f36cc28af9001b33ea35cb3e82f1ec46ceacdf1d..588d2680d1c6407353eca1e5571e8573c95301d3 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -163,7 +163,7 @@ end
γuuxx = 1.0/γxx
γuuzz = 1.0/γzz
- if (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < Ïmin)
+ if !c2p_set_atmosphere_on_failure && (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < Ïmin)
@warn "Require atmosphere but denied by mask" D Ïmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
TODO("Upsi")
elseif (c2p_freeze_atm == 0.0 && D < Ïmin) || (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
@@ -351,7 +351,7 @@ end
W = 1.0 / sqrt(1.0 - v^2)
Ï = D / W
end
- if (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && Ï < Ïmin)
+ if !c2p_set_atmosphere_on_failure && (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && Ï < Ïmin)
@warn "Require atmosphere but denied by mask" D Ïmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
TODO("Upsi")
elseif (c2p_freeze_atm == 0.0 && Ï < Ïmin) || (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
diff --git a/src/GRHD/doublecartoon.jl b/src/GRHD/doublecartoon.jl
index 551b4f5320918dd81e201d1c286cf7069d0372b7..31caea98dc4692206adc73d7febc89c802bc828d 100644
--- a/src/GRHD/doublecartoon.jl
+++ b/src/GRHD/doublecartoon.jl
@@ -570,8 +570,8 @@ end
@with_signature function stop_atmosphere_evolution_doublecartoon(eq::Equation)
- @accepts rhs_rD, rhs_rSx, rhs_rτ, c2p_reset_atm
- if c2p_reset_atm > 0
+ @accepts rhs_rD, rhs_rSx, rhs_rτ, c2p_reset_atm, c2p_freeze_atm
+ if c2p_reset_atm > 0 && c2p_freeze_atm > 0
rhs_rD = rhs_rSx = rhs_rτ = 0.0
end
@returns rhs_rD, rhs_rSx, rhs_rτ
diff --git a/src/GRHD/rescaled_spherical1d.jl b/src/GRHD/rescaled_spherical1d.jl
index f3e1c5d1c177f5648806c2e1b039c571b4b94825..166970316bbb8a4a004ad9d86eab5544b8385176 100644
--- a/src/GRHD/rescaled_spherical1d.jl
+++ b/src/GRHD/rescaled_spherical1d.jl
@@ -321,8 +321,8 @@ end
@with_signature function stop_atmosphere_evolution_rescaled_spherical1d(eq::Equation)
- @accepts rhs_rD, rhs_rSr, rhs_rτ, c2p_reset_atm
- if c2p_reset_atm > 0
+ @accepts rhs_rD, rhs_rSr, rhs_rτ, c2p_reset_atm, c2p_freeze_atm
+ if c2p_reset_atm > 0 && c2p_freeze_atm > 0
rhs_rD = rhs_rSr = rhs_rτ = 0.0
end
@returns rhs_rD, rhs_rSr, rhs_rτ
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 33984c2704e227ac24c56024599542a6db0b1494..abe5a67b080438695720fe7fd0a2285a7b52ebb3 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -108,7 +108,7 @@ function rhs_fv_central!(mesh, P::Project{:spherical1d})
bdry_D, bdry_Sr, bdry_Ï„ = get_bdry_variables(cache)
@unpack dt = get_global_variables(cache)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -145,7 +145,7 @@ function rhs_fv_central!(mesh, P::Project{:rescaled_spherical1d})
broadcast_volume!(unscale_conservatives, equation, mesh)
impose_symmetry!(P, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_rescaled_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
@@ -182,7 +182,7 @@ function rhs_fv_central!(mesh, P::Project{:doublecartoon})
bdry_rD, bdry_rSx, bdry_rτ = get_bdry_variables(cache)
broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -346,7 +346,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
D .= sD
Sr .= sSr
τ .= sτ
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -639,7 +639,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:doublecartoon}, state)
rSx .= sSx
rτ .= sτ
broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1177,7 +1177,7 @@ function rhs!(mesh::Mesh1d, P::Project{:valencia1d}, hrsc::Nothing)
broadcast_volume!(unscale_conservatives_valencia1d, equation, mesh)
impose_symmetry!(P, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_valencia1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_valencia1d, equation, mesh)
@@ -1234,7 +1234,7 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
bdry_D, bdry_Sr, bdry_Ï„,
bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -1325,7 +1325,7 @@ function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
broadcast_volume!(unscale_conservatives, equation, mesh)
impose_symmetry!(P, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_rescaled_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
@@ -1377,7 +1377,7 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::Nothing)
bdry_max_v, bdry_vx, bdry_p = get_bdry_variables(cache)
broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1454,7 +1454,7 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractReconstruction)
HRSC.reconstruct!(Sr, flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
HRSC.reconstruct!(Ï„, flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim, equation, mesh)
@@ -1550,7 +1550,7 @@ function rhs!(mesh, P::Project{:doublecartoon}, hrsc::HRSC.AbstractReconstructio
# HRSC.reconstruct!(rτ, flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1628,7 +1628,7 @@ function rhs!(mesh, P::Project{:valencia1d}, hrsc::HRSC.AbstractArtificialViscos
broadcast_volume!(unscale_conservatives_valencia1d, equation, mesh)
impose_symmetry!(P, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_valencia1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_valencia1d, equation, mesh)
@@ -1711,15 +1711,17 @@ 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
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_spherical1d, equation, mesh)
@@ -1733,6 +1735,12 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
dg1d.interpolate_face_data!(mesh, vr, bdry_vr)
dg1d.interpolate_face_data!(mesh, p, bdry_p)
+ if P.prms.atm_equalize_on_interface
+ 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)
+ end
+
broadcast_volume!(flux_source_spherical1d, equation, mesh)
impose_symmetry_sources!(P, mesh)
broadcast_faces!(llf_spherical1d, equation, mesh)
@@ -1864,7 +1872,7 @@ function rhs!(mesh::Mesh1d, P::Project{:doublecartoon}, hrsc::AbstractArtificial
bdry_rhs_rD, bdry_rhs_rSx, bdry_rhs_rτ = get_bdry_variables(cache)
broadcast_volume!(unscale_conservatives_doublecartoon, equation, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_doublecartoon_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_doublecartoon, equation, mesh)
@@ -1967,7 +1975,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
@@ -1998,7 +2005,7 @@ function rhs!(mesh::Mesh2d, P::Project{:cartoon})
bdry_max_v, bdry_vx, bdry_vz, bdry_p = get_bdry_variables(cache)
broadcast_volume!(unscale_conservatives_cartoon, equation, mesh)
- if P.prms.c2p_enforce_causal_atm
+ if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_cartoon_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_cartoon, equation, mesh)
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index 6e85e21af62693f7a56d5caf4671bbee2d30f04b..3c76c8bb3320a4cf0c4dbaa76d53108b3e526a44 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -26,6 +26,7 @@ function Project(env::Environment, prms)
:tov
end
+ atm_equalize_on_interface = prms["GRHD"]["atm_equalize_on_interface"]
av_drag = prms["GRHD"]["av_drag"]
av_sensor_abslog_D = Bool(prms["GRHD"]["av_sensor_abslog_D"])
slope_limiter_method = Symbol(prms["GRHD"]["slope_limiter_method"])
@@ -45,7 +46,7 @@ function Project(env::Environment, prms)
bc = prms["GRHD"]["bc"]
fixedprms = (; av_regularization=:covariant, id_smooth=true,
bernstein, slope_limiter_method, slope_limiter_tvb_M,
- c2p_dynamic_atm, atm_evolve,
+ c2p_dynamic_atm, atm_evolve, atm_equalize_on_interface,
c2p_set_atmosphere_on_failure, c2p_enforce_causal_atm, c2p_enforce_atm,
av_drag, av_sensor_abslog_D, problem, freeze_vars,
id, bc, hrsc, fv_numerical_flux)
@@ -479,7 +480,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 +491,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 ca7ac131ac544db34301d62fc4f9bb414554fbd7..ce0bfd81c01a6f52e861f86be485f74602d8a8b1 100644
--- a/src/GRHD/spherical1d.jl
+++ b/src/GRHD/spherical1d.jl
@@ -433,9 +433,41 @@ end
@with_signature function stop_atmosphere_evolution_spherical1d(eq::Equation)
- @accepts rhs_D, rhs_Sr, rhs_Ï„, c2p_reset_atm
- if c2p_reset_atm > 0
+ @accepts rhs_D, rhs_Sr, rhs_Ï„, c2p_reset_atm, c2p_freeze_atm
+ if c2p_reset_atm > 0 && c2p_freeze_atm > 0
rhs_D = rhs_Sr = rhs_Ï„ = 0.0
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
diff --git a/src/GRHD/valencia1d.jl b/src/GRHD/valencia1d.jl
index 075eb6aa73f81eda4ba20242b28748a21b7b9514..9a7153b46afaba647955b7acb693c3642175c9ef 100644
--- a/src/GRHD/valencia1d.jl
+++ b/src/GRHD/valencia1d.jl
@@ -440,8 +440,8 @@ end
@with_signature function stop_atmosphere_evolution_valencia1d(eq::Equation)
- @accepts rhs_rD, rhs_rSr, rhs_rτ, c2p_reset_atm
- if c2p_reset_atm > 0
+ @accepts rhs_rD, rhs_rSr, rhs_rτ, c2p_reset_atm, c2p_freeze_atm
+ if c2p_reset_atm > 0 && c2p_freeze_atm > 0
rhs_rD = rhs_rSr = rhs_rτ = 0.0
end
@returns rhs_rD, rhs_rSr, rhs_rτ