From e2b1a670a869cc12b00ad81ad86634bb65b63eb5 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Mon, 12 Aug 2024 16:32:11 +0000
Subject: [PATCH] GRHD: fix atmosphere freezing and evolution logic
(https://git.tpi.uni-jena.de/dg/dg1d.jl/-/merge_requests/199)
Fixes the following issues with the current implementation:
- `update_atm_domain_of_dependence!` did not distinguish between points
1. whose neighbors are also in the bulk,
2. where one neighbor is on an interface,
3. where the point itself is on an interface.
These cases are already properly handled in 2d.
- `stop_atmosphere_evolution_*` methods did not account for the `c2p_freeze_atm` flag, which essentially prevented the atmosphere to ever evolve.
- `update_atm_domain_of_dependence!` was activated by setting `c2p_dynamical_atm`, but it should be controlled by `c2p_enforce_causal_atm`
---
Additional fixes with less impact:
- remove unused parameter `id_atmosphere_density` from `GRHD` section
- `c2p_set_atmosphere_on_failure` should really suppress any errors appearing in the cons2prim
---
src/GRHD/GRHD.jl | 17 ++++----
src/GRHD/callbacks.jl | 69 ++++++++++++++++++++++++++------
src/GRHD/cartoon.jl | 4 +-
src/GRHD/cons2prim.jl | 4 +-
src/GRHD/doublecartoon.jl | 4 +-
src/GRHD/rescaled_spherical1d.jl | 4 +-
src/GRHD/rhs.jl | 43 +++++++++++---------
src/GRHD/setup.jl | 9 +++--
src/GRHD/spherical1d.jl | 36 ++++++++++++++++-
src/GRHD/valencia1d.jl | 4 +-
10 files changed, 141 insertions(+), 53 deletions(-)
diff --git a/src/GRHD/GRHD.jl b/src/GRHD/GRHD.jl
index f5c63495..a3174a8b 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 1b302c14..9192c61c 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 b81f5cac..37bad364 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 f36cc28a..588d2680 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 551b4f53..31caea98 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 f3e1c5d1..16697031 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 33984c27..abe5a67b 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 6e85e21a..3c76c8bb 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 ca7ac131..ce0bfd81 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 075eb6aa..9a7153b4 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τ
--
GitLab