From 4531207ea0f36a71768f6e96d7c6d960e107aa96 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Mon, 19 Feb 2024 12:38:35 +0000
Subject: [PATCH] Improvements to TOV evolutions, pt 2 (dg/dg1d.jl!113)
* wip: stash some tov examples
* update grhd ref tests to use c2p_freeze_atm_reset = false to restore previous behavior
* init c2p_freeze_atm
* fixup logic of freezing atm
* introduce c2p_freeze_atm flag which allows to define when an atmosphere can propagate
* fix time step computation to use maximum of smoothed_mu
* rescale smoothed viscosity by the extrema of D
* add more example parfiles for tov case
* enable c2p_freeze_atm_reset option
* add cons2prim_rescaled_spherical1d_freeze_flags
* add c2p_freeze_atm_reset option
* document perturbation_rho parameter
* update examples
* experiment with partially frozen c2p flags
* don't bernstein in AV computation
* add c2p with partially frozen reset flags
* use buffers for logarithmic AV indicators
* remove isinvalid
* remove dead code
* avoid splat constructor from CallbackSet
* add punctuation to some parameter doc strings
---
examples/grhd_tov/avmda_spherical1d.toml | 52 +++++--
examples/grhd_tov/avmda_spherical1d_k17.toml | 85 ++++++++++
.../grhd_tov/avmda_spherical1d_k17_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k18.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k18_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k19.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k19_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k20.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k20_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k21.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k21_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k22.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k22_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k23.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k23_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k24.toml | 83 ++++++++++
.../avmda_spherical1d_k24_nofreeze.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_k24_pt.toml | 86 ++++++++++
examples/grhd_tov/avmda_spherical1d_k48.toml | 83 ++++++++++
.../avmda_spherical1d_k48_nofreeze.toml | 84 ++++++++++
examples/grhd_tov/avmda_spherical1d_k62.toml | 83 ++++++++++
.../avmda_spherical1d_k62_nofreeze.toml | 84 ++++++++++
.../grhd_tov/avmda_spherical1d_mirrored.toml | 60 ++++---
examples/grhd_tov/dg_rescaled.toml | 40 +++++
examples/grhd_tov/dg_spherical1d_k.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k17.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k17_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k18.toml | 38 +++++
examples/grhd_tov/dg_spherical1d_k18_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k19.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k19_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k20.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k20_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k21.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k21_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k22.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k22_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k23.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k23_pt.toml | 32 ++++
examples/grhd_tov/dg_spherical1d_k24.toml | 31 ++++
examples/grhd_tov/dg_spherical1d_k24_pt.toml | 32 ++++
src/GRHD/GRHD.jl | 10 +-
src/GRHD/callbacks.jl | 77 ++++++---
src/GRHD/cons2prim.jl | 147 ++++++------------
src/GRHD/initialdata.jl | 18 +++
src/GRHD/rhs.jl | 48 ++++--
src/GRHD/setup.jl | 9 +-
src/callbacks.jl | 3 +-
src/evolve.jl | 18 +--
src/parameters.jl | 12 +-
.../fv_grhd_bondi_accretion.toml | 1 +
.../fv_superbee_grhd_bondi_accretion.toml | 2 +-
.../grhd_bondi_accretion.toml | 1 +
.../grhd_tov_spherical1d.toml | 1 +
54 files changed, 2631 insertions(+), 189 deletions(-)
create mode 100644 examples/grhd_tov/avmda_spherical1d_k17.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k17_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k18.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k18_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k19.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k19_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k20.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k20_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k21.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k21_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k22.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k22_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k23.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k23_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k24.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k24_nofreeze.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k24_pt.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k48.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k48_nofreeze.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k62.toml
create mode 100644 examples/grhd_tov/avmda_spherical1d_k62_nofreeze.toml
create mode 100644 examples/grhd_tov/dg_rescaled.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k17.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k17_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k18.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k18_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k19.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k19_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k20.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k20_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k21.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k21_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k22.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k22_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k23.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k23_pt.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k24.toml
create mode 100644 examples/grhd_tov/dg_spherical1d_k24_pt.toml
diff --git a/examples/grhd_tov/avmda_spherical1d.toml b/examples/grhd_tov/avmda_spherical1d.toml
index a74ed718..805f9fbe 100644
--- a/examples/grhd_tov/avmda_spherical1d.toml
+++ b/examples/grhd_tov/avmda_spherical1d.toml
@@ -1,5 +1,3 @@
-# runs nicely on: d4f313169d021fbc87e9911fca6a19a0ac8d412e
-
[EquationOfState]
eos = "polytrope"
polytrope_k = 100.0
@@ -12,6 +10,7 @@ id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
atm_factor = 1e-8
formulation = "spherical1d"
av_drag = 0.99
+# c2p_freeze_atm_reset = false
# training_wheels = true
[Mesh]
@@ -21,39 +20,58 @@ k = 16
basis = "lgl"
periodic = false
-
[Output]
# every_iteration = 1
-# every_dt = 2.5
-# every_dt = 5.0
-aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
-# aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
variables = [ "D", "Sr", "Ï„",
- # "p", "Ï", "ϵ", "vr",
"p", "Ï", "vr", "ÏhW2",
- # "p", "Ï", "vr", "v",
- # "src_D", "src_Sr", "src_Ï„",
- # "rhs_D", "rhs_Sr", "rhs_Ï„",
- # "flr_D", "flr_Sr", "flr_Ï„",
- # "smoothed_mu",
"mu", "smoothed_mu",
"c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
"max_v",
- # "E", #"flr_E"
- # "A", "∂Ar", "B", "∂Br"
]
-enable1d = true
[Evolution]
cfl = 0.8
# tend = 500.0
tend = 1000.0
+# tend = 12.0
[HRSC]
method = "av"
av_method = "mda"
mda_cmax = 100.0
-mda_n_low = 2.0
+mda_n_low = 3.0
mda_n_high = 4.0
# mda_n_low = 2.0
# mda_n_high = 3.0
diff --git a/examples/grhd_tov/avmda_spherical1d_k17.toml b/examples/grhd_tov/avmda_spherical1d_k17.toml
new file mode 100644
index 00000000..3ab403f8
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k17.toml
@@ -0,0 +1,85 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 17
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k17_pt.toml b/examples/grhd_tov/avmda_spherical1d_k17_pt.toml
new file mode 100644
index 00000000..8bcc3747
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k17_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 17
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k18.toml b/examples/grhd_tov/avmda_spherical1d_k18.toml
new file mode 100644
index 00000000..5aa66c8a
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k18.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 18
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k18_pt.toml b/examples/grhd_tov/avmda_spherical1d_k18_pt.toml
new file mode 100644
index 00000000..8e7562b2
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k18_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 18
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k19.toml b/examples/grhd_tov/avmda_spherical1d_k19.toml
new file mode 100644
index 00000000..42706682
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k19.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 19
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k19_pt.toml b/examples/grhd_tov/avmda_spherical1d_k19_pt.toml
new file mode 100644
index 00000000..2699bca2
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k19_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 19
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k20.toml b/examples/grhd_tov/avmda_spherical1d_k20.toml
new file mode 100644
index 00000000..f5a5d217
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k20.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 20
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k20_pt.toml b/examples/grhd_tov/avmda_spherical1d_k20_pt.toml
new file mode 100644
index 00000000..4b9db624
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k20_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 20
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k21.toml b/examples/grhd_tov/avmda_spherical1d_k21.toml
new file mode 100644
index 00000000..5496d490
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k21.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 21
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k21_pt.toml b/examples/grhd_tov/avmda_spherical1d_k21_pt.toml
new file mode 100644
index 00000000..58ae54f0
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k21_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 21
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k22.toml b/examples/grhd_tov/avmda_spherical1d_k22.toml
new file mode 100644
index 00000000..4f68d6ba
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k22.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 22
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k22_pt.toml b/examples/grhd_tov/avmda_spherical1d_k22_pt.toml
new file mode 100644
index 00000000..33749340
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k22_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 22
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k23.toml b/examples/grhd_tov/avmda_spherical1d_k23.toml
new file mode 100644
index 00000000..919c14fd
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k23.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 23
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k23_pt.toml b/examples/grhd_tov/avmda_spherical1d_k23_pt.toml
new file mode 100644
index 00000000..8826fa1d
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k23_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 23
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k24.toml b/examples/grhd_tov/avmda_spherical1d_k24.toml
new file mode 100644
index 00000000..7334400e
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k24.toml
@@ -0,0 +1,83 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 24
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k24_nofreeze.toml b/examples/grhd_tov/avmda_spherical1d_k24_nofreeze.toml
new file mode 100644
index 00000000..d9bfc250
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k24_nofreeze.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 24
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k24_pt.toml b/examples/grhd_tov/avmda_spherical1d_k24_pt.toml
new file mode 100644
index 00000000..8e1ea535
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k24_pt.toml
@@ -0,0 +1,86 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+perturbation_rho = true
+# c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 24
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k48.toml b/examples/grhd_tov/avmda_spherical1d_k48.toml
new file mode 100644
index 00000000..96846dab
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k48.toml
@@ -0,0 +1,83 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 48
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k48_nofreeze.toml b/examples/grhd_tov/avmda_spherical1d_k48_nofreeze.toml
new file mode 100644
index 00000000..72ef2e9e
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k48_nofreeze.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 48
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k62.toml b/examples/grhd_tov/avmda_spherical1d_k62.toml
new file mode 100644
index 00000000..ffbdf451
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k62.toml
@@ -0,0 +1,83 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 62
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_k62_nofreeze.toml b/examples/grhd_tov/avmda_spherical1d_k62_nofreeze.toml
new file mode 100644
index 00000000..c205b638
--- /dev/null
+++ b/examples/grhd_tov/avmda_spherical1d_k62_nofreeze.toml
@@ -0,0 +1,84 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+av_drag = 0.99
+c2p_freeze_atm_reset = false
+# training_wheels = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 62
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "mu", "smoothed_mu",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.8
+# tend = 500.0
+tend = 1000.0
+# tend = 12.0
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 100.0
+mda_n_low = 3.0
+mda_n_high = 4.0
+# mda_n_low = 2.0
+# mda_n_high = 3.0
+# av_method = "entropy"
+# entropy_ce = 100.0
+# entropy_cmax = 100.0
+# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/avmda_spherical1d_mirrored.toml b/examples/grhd_tov/avmda_spherical1d_mirrored.toml
index 47f685ef..5f5e9f82 100644
--- a/examples/grhd_tov/avmda_spherical1d_mirrored.toml
+++ b/examples/grhd_tov/avmda_spherical1d_mirrored.toml
@@ -1,5 +1,3 @@
-# runs nicely on: d4f313169d021fbc87e9911fca6a19a0ac8d412e
-
[EquationOfState]
eos = "polytrope"
polytrope_k = 100.0
@@ -12,48 +10,67 @@ id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
atm_factor = 1e-8
formulation = "spherical1d"
av_drag = 0.99
-training_wheels = true
+# training_wheels = true
[Mesh]
range = [ -20.0, 20.0 ]
n = 5
-k = 31
+k = 49
basis = "lgl"
periodic = false
-
[Output]
# every_iteration = 1
-# every_dt = 2.5
-# every_dt = 5.0
-aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
-# aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# # every_dt = 2.5
+# # every_dt = 5.0
+# # aligned_ts = "$(collect(range(0.1,1000.0,step=0.1)))"
+# # aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# # aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+# variables = [ "D", "Sr", "Ï„",
+# # "p", "Ï", "ϵ", "vr",
+# "p", "Ï", "vr", "ÏhW2",
+# # "p", "Ï", "vr", "v",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "rhs_D", "rhs_Sr", "rhs_Ï„",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "smoothed_mu",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# # "E", #"flr_E"
+# # "A", "∂Ar", "B", "∂Br"
+# ]
+# enable1d = true
+# substep_every_dt = 5.0
+
+# substep_every_iteration = 1
+# substep_variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "mu", "smoothed_mu",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+
+# every_iteration = 1
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
variables = [ "D", "Sr", "Ï„",
- # "p", "Ï", "ϵ", "vr",
"p", "Ï", "vr", "ÏhW2",
- # "p", "Ï", "vr", "v",
- # "src_D", "src_Sr", "src_Ï„",
- # "rhs_D", "rhs_Sr", "rhs_Ï„",
- # "flr_D", "flr_Sr", "flr_Ï„",
- # "smoothed_mu",
- # "mu", "smoothed_mu",
+ "mu", "smoothed_mu",
"c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
"max_v",
- # "E", #"flr_E"
- # "A", "∂Ar", "B", "∂Br"
]
-enable1d = true
[Evolution]
cfl = 0.8
# tend = 500.0
tend = 1000.0
+# tend = 12.0
[HRSC]
method = "av"
av_method = "mda"
mda_cmax = 100.0
-mda_n_low = 2.0
+mda_n_low = 3.0
mda_n_high = 4.0
# mda_n_low = 2.0
# mda_n_high = 3.0
@@ -61,3 +78,6 @@ mda_n_high = 4.0
# entropy_ce = 100.0
# entropy_cmax = 100.0
# av_smoother_order = 1
+
+# [Log]
+# progress_stdout = false
diff --git a/examples/grhd_tov/dg_rescaled.toml b/examples/grhd_tov/dg_rescaled.toml
new file mode 100644
index 00000000..205e41b4
--- /dev/null
+++ b/examples/grhd_tov/dg_rescaled.toml
@@ -0,0 +1,40 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "rescaled_spherical1d"
+c2p_freeze_atm_reset = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 16
+basis = "lgl"
+periodic = false
+
+[Output]
+every_iteration = 1
+# every_dt = 2.5
+# every_dt = 5.0
+# aligned_ts = "$(collect(range(5.0,5000.0,step=5.0)))"
+# aligned_ts = "$(collect(range(0.2,120.0,step=0.2)))"
+variables = [ "D", "Sr", "Ï„",
+ # "p", "Ï", "ϵ", "vr",
+ "p", "Ï", "vr",
+ # "p", "Ï", "vr", "v",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "c2p_freeze_atm",
+ "max_v",
+ "E",
+ ]
+enable1d = true
+
+[Evolution]
+cfl = 0.5
+tend = 2000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k.toml b/examples/grhd_tov/dg_spherical1d_k.toml
new file mode 100644
index 00000000..36b6c358
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 17
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k17.toml b/examples/grhd_tov/dg_spherical1d_k17.toml
new file mode 100644
index 00000000..2797e7be
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k17.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 17
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k17_pt.toml b/examples/grhd_tov/dg_spherical1d_k17_pt.toml
new file mode 100644
index 00000000..36b6c358
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k17_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 17
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k18.toml b/examples/grhd_tov/dg_spherical1d_k18.toml
new file mode 100644
index 00000000..ed7b35b5
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k18.toml
@@ -0,0 +1,38 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 18
+basis = "lgl"
+periodic = false
+
+[Output]
+
+# aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+# variables = [ "D", "Sr", "Ï„",
+# "p", "Ï", "vr", "ÏhW2",
+# "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+# "max_v",
+# ]
+substep_every_iteration = 1
+substep_variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+# tend = 1000.0
+tend = 100.0
diff --git a/examples/grhd_tov/dg_spherical1d_k18_pt.toml b/examples/grhd_tov/dg_spherical1d_k18_pt.toml
new file mode 100644
index 00000000..cbe08744
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k18_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 18
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k19.toml b/examples/grhd_tov/dg_spherical1d_k19.toml
new file mode 100644
index 00000000..ac9eaecf
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k19.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 19
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k19_pt.toml b/examples/grhd_tov/dg_spherical1d_k19_pt.toml
new file mode 100644
index 00000000..9093fca4
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k19_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 19
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k20.toml b/examples/grhd_tov/dg_spherical1d_k20.toml
new file mode 100644
index 00000000..0d835a58
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k20.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 20
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k20_pt.toml b/examples/grhd_tov/dg_spherical1d_k20_pt.toml
new file mode 100644
index 00000000..d0430cf8
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k20_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 20
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k21.toml b/examples/grhd_tov/dg_spherical1d_k21.toml
new file mode 100644
index 00000000..9754dc7d
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k21.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 21
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k21_pt.toml b/examples/grhd_tov/dg_spherical1d_k21_pt.toml
new file mode 100644
index 00000000..f6090281
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k21_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 21
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k22.toml b/examples/grhd_tov/dg_spherical1d_k22.toml
new file mode 100644
index 00000000..f9b5de4e
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k22.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 22
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k22_pt.toml b/examples/grhd_tov/dg_spherical1d_k22_pt.toml
new file mode 100644
index 00000000..0031befd
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k22_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 22
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k23.toml b/examples/grhd_tov/dg_spherical1d_k23.toml
new file mode 100644
index 00000000..8858323c
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k23.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 23
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k23_pt.toml b/examples/grhd_tov/dg_spherical1d_k23_pt.toml
new file mode 100644
index 00000000..dbf16ad1
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k23_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 23
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k24.toml b/examples/grhd_tov/dg_spherical1d_k24.toml
new file mode 100644
index 00000000..769eb5fa
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k24.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 24
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/examples/grhd_tov/dg_spherical1d_k24_pt.toml b/examples/grhd_tov/dg_spherical1d_k24_pt.toml
new file mode 100644
index 00000000..d2268039
--- /dev/null
+++ b/examples/grhd_tov/dg_spherical1d_k24_pt.toml
@@ -0,0 +1,32 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+perturbation_rho = true
+
+[Mesh]
+range = [ 0.0, 20.0 ]
+n = 5
+k = 24
+basis = "lgl"
+periodic = false
+
+[Output]
+
+aligned_ts = "$(collect(range(1.0,1000.0,step=1.0)))"
+variables = [ "D", "Sr", "Ï„",
+ "p", "Ï", "vr", "ÏhW2",
+ "c2p_reset_ϵ", "c2p_reset_atm", "c2p_limit_vr", "c2p_init_admissible",
+ "max_v",
+ ]
+
+[Evolution]
+cfl = 0.25
+tend = 1000.0
diff --git a/src/GRHD/GRHD.jl b/src/GRHD/GRHD.jl
index 458cebdf..f35ff19b 100644
--- a/src/GRHD/GRHD.jl
+++ b/src/GRHD/GRHD.jl
@@ -99,7 +99,8 @@ dg1d.@parameters GRHD begin
@check training_wheels isa Bool
"""
- perturbation_rho = true
+ If `true` adds a random perturbation with amplitude `1e-3*max(Ï)` to the rest mass-energy density Ï.
+ TODO Add parameter to control amplitude.
"""
perturbation_rho = false
@check perturbation_rho isa Bool
@@ -133,6 +134,13 @@ dg1d.@parameters GRHD begin
slope_limiter_tvb_M = 150.0
@check slope_limiter_tvb_M > 0.0
+ """
+ If `true` then atmosphere is only enforced after full time steps.
+ This has a stabilizing effect onto AV regularizations.
+ """
+ c2p_freeze_atm_reset = true
+ @check c2p_freeze_atm_reset isa Bool
+
end
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index f06ec21a..7cb1f5d5 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -18,8 +18,34 @@ end
#######################################################################
+function update_atm_domain_of_dependence!(mesh::Mesh1d)
+ @unpack c2p_freeze_atm, c2p_reset_atm = get_static_variables(mesh)
+ freeze, reset = c2p_freeze_atm, c2p_reset_atm
+ for i in 2:length(reset)-1
+ rl, rc, rr = reset[i-1], reset[i], reset[i+1]
+ if rl == rc == rr
+ # atmosphere can't spread within next timestep
+ freeze[i-1] = freeze[i] = freeze[i+1] = 1
+ continue
+ else
+ freeze[i-1] = freeze[i] = freeze[i+1] = 0
+ end
+ end
+end
+
+
function callback_equation(state_u, state_t, env, P, isperiodic, eq::Equation)
wrap_dynamic_variables!(env.cache, state_u)
+ if P.prms.c2p_freeze_atm_reset
+ if formulation(P) === :spherical1d
+ broadcast_volume!(cons2prim, eq, env.mesh)
+ elseif formulation(P) === :rescaled_spherical1d
+ broadcast_volume!(cons2prim_rescaled_spherical1d, eq, env.mesh)
+ else
+ TODO()
+ end
+ update_atm_domain_of_dependence!(env.mesh)
+ end
callback_equation(state_t, isperiodic, P.equation, env.mesh)
end
function callback_equation(state_t, isperiodic, eq::Equation, mesh::Mesh1d)
@@ -101,14 +127,20 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
end
# display(extrema(cellmax_v))
+ @unpack abslog_D, abslog_Ï„ = get_static_variables(cache)
+ @. abslog_D = log10.(abs.(D))
+ @. abslog_Ï„ = log10.(abs.(Ï„))
+
fill!(tmp_mu, 0.0)
HRSC.compute_viscosity!(
mu, tmp_mu,
+ (abslog_D, abslog_Ï„), cellmax_v,
+ # (log10.(abs.(D)),log10.(abs.(Ï„))), cellmax_v,
# (log10.(abs.(D)),Sr,log10.(abs.(Ï„))), cellmax_v,
# (log10.(abs.(D)),), cellmax_v,
# (log10.(abs.(D)),log10.(abs.(Sr)),log10.(abs.(Ï„))), cellmax_v,
# (D,Sr,Ï„), cellmax_v,
- (D,Ï„), cellmax_v,
+ # (D,Ï„), cellmax_v,
# (vr,), cellmax_v,
hrsc)
# Ïmin = equation.atm_density * equation.atm_factor
@@ -132,20 +164,20 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
# end
mu[1] = mu[2] = mu[end] = 0 # ignore boundaries
- prms = parameters(:HRSC)
- prms["method"] = "filter"
- prms["filter_method"] = "bernstein"
- dg1d.check_parameters(:HRSC, prms)
- bernstein_hrsc = HRSC.make_HRSC(mesh, prms)
- flag = mu .> 0.0
- @unpack D, Sr, Ï„ = get_dynamic_variables(mesh)
- logD = log.(abs.(D))
- logτ = log.(abs.(τ))
- HRSC.reconstruct!(logD, flag, bernstein_hrsc, isperiodic=false)
- HRSC.reconstruct!(Sr, flag, bernstein_hrsc, isperiodic=false)
- HRSC.reconstruct!(logτ, flag, bernstein_hrsc, isperiodic=false)
- D .= exp.(logD)
- τ .= exp.(logτ)
+ # prms = parameters(:HRSC)
+ # prms["method"] = "filter"
+ # prms["filter_method"] = "bernstein"
+ # dg1d.check_parameters(:HRSC, prms)
+ # bernstein_hrsc = HRSC.make_HRSC(mesh, prms)
+ # flag = mu .> 0.0
+ # @unpack D, Sr, Ï„ = get_dynamic_variables(mesh)
+ # logD = log.(abs.(D))
+ # logτ = log.(abs.(τ))
+ # HRSC.reconstruct!(logD, flag, bernstein_hrsc, isperiodic=false)
+ # HRSC.reconstruct!(Sr, flag, bernstein_hrsc, isperiodic=false)
+ # HRSC.reconstruct!(logτ, flag, bernstein_hrsc, isperiodic=false)
+ # D .= exp.(logD)
+ # τ .= exp.(logτ)
# broadcast_volume!(cons2prim, P.equation, mesh)
# broadcast_volume!(maxspeed, P.equation, mesh)
# broadcast_volume!(flux_source_spherical1d, P.equation, mesh)
@@ -170,9 +202,18 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Smoothed
@unpack cache, mesh = env
@unpack smoothed_mu = get_static_variables(cache)
hrsc.smoother(smoothed_mu, mu, mesh, false)
- max_D = maximum(D)
- min_D = minimum(D)
- smoothed_mu .*= (D .- min_D) / max_D
+ # min_D, max_D = extrema(D)
+ # max_D = maximum(D)
+ # min_D = minimum(D)
+ # smoothed_mu .*= (D .- min_D) / max_D
+ # @. smoothed_mu *= D * min_D / max_D
+ K, = env.mesh.tree.dims
+ mat_D = dg1d.vreshape(D, layout(env.mesh))
+ mat_smoothed_mu = dg1d.vreshape(smoothed_mu, layout(env.mesh))
+ for k in 1:K
+ min_D, max_D = extrema(@views mat_D[:,k])
+ @. @views mat_smoothed_mu[:,k] *= mat_D[:,k] * min_D / max_D
+ end
end
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index a1e90876..fa6c5e25 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -1,121 +1,59 @@
#######################################################################
-# Conservative fixing #
+# cons2prim a la Kastaun et al 2021 #
#######################################################################
-@with_signature function conservative_fixing(equation::AbstractEquation)
- @accepts D, S, Ï„, r
-
- TODO()
-
- @unpack atm_density, atm_threshold, eos = equation
- # atm_density = equation.atmosphere.Ï
- # atm_threshold = equation.atmosphere.threshold
- Ïmin = atm_density * atm_threshold
-
- # B = 1.0 # in special relativity
-
- #####
- # conservative variable fixing following Deppe et al. arXiv:2109.12033
- ######
-
- fixed = false
-
- if D < Ïmin || Ï„ < 0.0
- D = Ïmin
- Ï„ = 1e-12
- fixed = true
- end
-
- # if D < 0.001 * Ïmin
- # # @warn "fixing D $D @ r = $r"
- # D = Ïmin
- # fixed = true
- # end
- #
- # if Ï„ < 0.0
- # # @warn "fixing Ï„ $Ï„ @ r = $r"
- # Ï„ = 1e-12
- # fixed = true
- # end
- # fix Ï„ such that v < 0.999
- # if Ï„ < 0.0
- # @warn "fixing Ï„ $Ï„ @ r = $r"
- # # What = 1 + Ï„ / D = 1 / sqrt(1 - v^2)
- # v2 = 0.999
- # Ï„ = max(D * (1.0 / sqrt(1.0 - v2) - 1), 1e-12)
- # display(Ï„)
- # fixed = true
- # end
-
- # (48)-(50) in arXiv:2109.12033 we put B=0, ̂μ = 0
- # solution to (52) is then
- What = 1.0 + Ï„ / D
-
- S2 = S^2
- D2 = D^2
- ϵs = 1e-12
- factor = sqrt( (1.0-ϵs) * (What^2-1.0) * D2 / (S2 + D2 * 1e-16) )
- # @show factor, What
- if factor < 1.0 && What > 1.0
- # @warn "fixing ..."
- Sold = S
- S = 0.99 * sign(S) * sqrt( (1.0-ϵs) * (What^2-1.0) * D2 - D2 * 1e-16 )
- # @warn "S = $Sold => $S @ r = $r"
- S2 = S^2
- factor = max(1.0, sqrt( (1.0-ϵs) * (What^2-1.0) * D2 / (S2 + D2 * 1e-16) ))
- # @show factor
- @assert factor >= 1.0
- fixed = true
- end
-
- # h0 = 1e-3
- # S2 = S^2 / B^2
- # r2 = S2 / D^2
- # z02 = r2 / h0^2
- # v02 = z02 / (1.0 + z02)
- # max_v2 = 0.999
- # if v02 > max_v2
- # S = sign(S) * sqrt(max_v2/(1-max_v2) * h0^2 * D^2)
- # end
-
- flag_consfixed = Float64(fixed)
-
- @returns D, S, Ï„, flag_consfixed
+@with_signature [simd=false] function cons2prim_freeze_flags(equation::AbstractEquation)
+ @accepts D, Sr, Ï„, r
+ @accepts c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm
+ @accepts γrr
+ D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, _ =
+ cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr, c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm)
+ @returns D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible
+end
+@with_signature [simd=false] function cons2prim_rescaled_spherical1d_freeze_flags(equation::AbstractEquation)
+ @accepts D, Sr, Ï„, r
+ @accepts c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm
+ @accepts B
+ γrr = 1/B^2
+ D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, _ =
+ cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr, c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm)
+ # in case something was adjusted by cons2prim
+ rD = r*B^3*D
+ rSr = r*B^3*Sr
+ rτ = r*B^3*τ
+ @returns D, Sr, Ï„, rD, rSr, rÏ„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible
end
-
-
-#######################################################################
-# cons2prim a la Kastaun et al 2021 #
-#######################################################################
-
@with_signature [simd=false] function cons2prim(equation::AbstractEquation)
@accepts D, Sr, Ï„, r
+ @accepts c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm
@accepts γrr
- D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible =
- cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
- @returns D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible
+ D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, c2p_freeze_atm =
+ cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr, c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm)
+ @returns D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, c2p_freeze_atm
end
@with_signature [simd=false] function cons2prim_rescaled_spherical1d(equation::AbstractEquation)
@accepts D, Sr, Ï„, r
+ @accepts c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm
@accepts B
γrr = 1/B^2
- D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible =
- cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr)
+ D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, c2p_freeze_atm =
+ cons2prim_kastaun_impl(equation, D, Sr, τ, r, γrr, c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm)
# in case something was adjusted by cons2prim
rD = r*B^3*D
rSr = r*B^3*Sr
rτ = r*B^3*τ
- @returns D, Sr, Ï„, rD, rSr, rÏ„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible
+ @returns D, Sr, Ï„, rD, rSr, rÏ„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, c2p_freeze_atm
end
-@inline function cons2prim_kastaun_impl(equation, D, Sr, τ, rcoord, γrr)
+@inline function cons2prim_kastaun_impl(equation, D, Sr, τ, rcoord, γrr,
+ c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm)
verbose = false
- c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible = 0.0, 0.0, 0.0, 1.0
+ c2p_reset_ϵ, c2p_init_admissible = 0.0, 1.0
@unpack atm_density, atm_threshold, eos = equation
Ïmin = atm_density * atm_threshold
@@ -234,7 +172,10 @@ end
error("Invalid state encountered!")
end
- if D < Ïmin
+ if (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < Ïmin)
+ @show D, Ïmin, rcoord, 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)
@assert eos isa EquationOfState.Polytrope
@unpack K, Gamma = eos
@@ -421,7 +362,10 @@ end
# enforce EOS bounds
- if D < Ïmin
+ if (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < Ïmin)
+ @show D, Ïmin, rcoord, 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)
# atmosphere II
error("Atmosphering after root finding")
@@ -467,5 +411,14 @@ end
W2 = 1/(1-v^2)
ÏhW2 = Ï*h*W2
- return D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible
+ return D, Sr, Ï„, Ï, vr, v, ϵ, p, ÏhW2, c2p_reset_ϵ, c2p_reset_atm, c2p_limit_vr, c2p_init_admissible, c2p_freeze_atm
+end
+
+
+@with_signature function determine_atmosphere(equation::Equation)
+ @accepts D
+ @unpack atm_density, atm_threshold, eos = equation
+ Ïmin = atm_density * atm_threshold
+ c2p_reset_atm = Float64(D<Ïmin)
+ @returns c2p_reset_atm
end
diff --git a/src/GRHD/initialdata.jl b/src/GRHD/initialdata.jl
index d8bc940a..ccf70919 100644
--- a/src/GRHD/initialdata.jl
+++ b/src/GRHD/initialdata.jl
@@ -122,6 +122,12 @@ function initialdata_equation!(::Val{:bondi_accretion}, env, P::Project{:spheric
@. Sr = W^2 * Ï * h * vr
@. Ï„ = W^2 * Ï * h - p - D
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(determine_atmosphere, equation, mesh)
+ else
+ @unpack c2p_freeze_atm = get_static_variables(mesh)
+ c2p_freeze_atm .= 0.0 # always on
+ end
broadcast_volume!(cons2prim, equation, mesh)
broadcast_volume!(maxspeed, equation, mesh)
@@ -303,6 +309,12 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:spherical1d}, prms)
@. Sr = W^2 * Ï * h * vr
@. Ï„ = W^2 * Ï * h - p - D
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(determine_atmosphere, equation, mesh)
+ else
+ @unpack c2p_freeze_atm = get_static_variables(mesh)
+ c2p_freeze_atm .= 0.0 # always on
+ end
broadcast_volume!(cons2prim, equation, mesh)
broadcast_volume!(maxspeed, equation, mesh)
@@ -457,6 +469,12 @@ function initialdata_equation!(::Val{:tov}, env, P::Project{:rescaled_spherical1
broadcast_volume!(unscale_conservatives, equation, mesh)
impose_symmetry!(P, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(determine_atmosphere, equation, mesh)
+ else
+ @unpack c2p_freeze_atm = get_static_variables(mesh)
+ c2p_freeze_atm .= 0.0 # always on
+ end
broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index c416a367..5d4e7992 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -37,8 +37,10 @@ function timestep(mesh, P::Project, hrsc::HRSC.AbstractArtificialViscosity)
compute_maxspeed(P, P.equation, mesh)
vmax = get_maxspeed(mesh, P.equation)
- @unpack mu = get_cell_variables(mesh)
- mumax = maximum(abs, mu)
+ # @unpack mu = get_cell_variables(mesh)
+ # mumax = maximum(abs, mu)
+ @unpack smoothed_mu = get_static_variables(mesh)
+ mumax = maximum(abs, smoothed_mu)
@unpack element = mesh
@unpack N = element
@@ -88,7 +90,11 @@ function rhs_fv_central!(P::Project{:spherical1d}, mesh::Mesh1d{FVElement})
bdry_D, bdry_Sr, bdry_Ï„ = get_bdry_variables(cache)
@unpack dt = get_global_variables(cache)
- broadcast_volume!(cons2prim, equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim, equation, mesh)
+ end
broadcast_volume!(maxspeed, equation, mesh)
broadcast_volume!(flux_source_spherical1d, equation, mesh)
broadcast_bdry!(fv_bdry_flux, equation, mesh)
@@ -116,7 +122,11 @@ function rhs_fv_central!(P::Project{:rescaled_spherical1d}, mesh::Mesh1d{FVEleme
broadcast_volume!(unscale_conservatives, equation, mesh)
impose_symmetry!(P, mesh)
- broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_rescaled_spherical1d_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+ end
broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
broadcast_volume!(flux_source_rescaled_spherical1d, equation, mesh)
broadcast_bdry!(fv_bdry_flux_rescaled_spherical1d, equation, mesh)
@@ -190,7 +200,11 @@ function step_fv_muscl!(P::Project, mesh::Mesh1d{FVElement}, state)
D .= sD
Sr .= sSr
τ .= sτ
- broadcast_volume!(cons2prim, P.equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim, equation, mesh)
+ end
broadcast_volume!(flux_source_spherical1d, P.equation, mesh)
broadcast_volume!(maxspeed, P.equation, mesh)
D .= bD
@@ -349,7 +363,11 @@ 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)
- broadcast_volume!(cons2prim, equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim, equation, mesh)
+ end
broadcast_volume!(maxspeed, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
@@ -394,7 +412,11 @@ function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
broadcast_volume!(unscale_conservatives, equation, mesh)
impose_symmetry!(P, mesh)
- broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_rescaled_spherical1d_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim_rescaled_spherical1d, equation, mesh)
+ end
broadcast_volume!(maxspeed_rescaled_spherical1d, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
@@ -435,7 +457,11 @@ 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)
- broadcast_volume!(cons2prim, equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim, equation, mesh)
+ end
broadcast_volume!(maxspeed, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
@@ -482,7 +508,11 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
bdry_rhs_D, bdry_rhs_Sr, bdry_rhs_Ï„,
bdry_smoothed_mu = get_bdry_variables(cache)
- broadcast_volume!(cons2prim, equation, mesh)
+ if P.prms.c2p_freeze_atm_reset
+ broadcast_volume!(cons2prim_freeze_flags, equation, mesh)
+ else
+ broadcast_volume!(cons2prim, equation, mesh)
+ end
broadcast_volume!(maxspeed, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index 4fa62b94..c10ed3a0 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -14,9 +14,11 @@ function Project(env::Environment, prms)
hrsc = Symbol(prms["GRHD"]["hrsc"])
slope_limiter_method = Symbol(prms["GRHD"]["slope_limiter_method"])
slope_limiter_tvb_M = prms["GRHD"]["slope_limiter_tvb_M"]
+ c2p_freeze_atm_reset = prms["GRHD"]["c2p_freeze_atm_reset"]
bernstein = BernsteinReconstruction(env.mesh)
fixedprms = (; av_regularization=:covariant, id_smooth=true,
bernstein, hrsc, slope_limiter_method, slope_limiter_tvb_M,
+ c2p_freeze_atm_reset,
av_drag, problem)
# construct Project
@@ -37,7 +39,7 @@ function Project(env::Environment, prms)
# setup callbacks
projectcb = make_callback(env, P, isperiodic(env.mesh))
- append!(env.callbacks, CallbackSet(projectcb.callbacks...))
+ append!(env.callbacks, CallbackSet(projectcb.callbacks))
return P, nothing
end
@@ -188,7 +190,7 @@ end
function register_analysis!(mesh)
register_variables!(mesh,
static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
- :c2p_init_admissible, :v)
+ :c2p_freeze_atm, :c2p_init_admissible, :v)
)
end
@@ -204,7 +206,8 @@ register_hrsc_rescaled_spherical1d!(mesh, hrsc::AbstractHRSC) = TODO(hrsc)
function register_hrsc_spherical1d!(mesh, hrsc::HRSC.AbstractArtificialViscosity)
register_variables!(mesh,
- static_variablenames = (:ldg_D, :ldg_Sr, :ldg_Ï„),
+ static_variablenames = (:ldg_D, :ldg_Sr, :ldg_Ï„,
+ :abslog_D, :abslog_Ï„),
bdry_variablenames = (:bdry_ldg_D, :bdry_ldg_Sr, :bdry_ldg_Ï„,
:bdry_smoothed_mu,
:bdry_rhs_D, :bdry_rhs_Sr, :bdry_rhs_Ï„),
diff --git a/src/callbacks.jl b/src/callbacks.jl
index a6a2a776..b7d5d445 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -184,8 +184,9 @@ Callbacks are called in the order they are queued.
"""
struct CallbackSet
callbacks::Vector{AbstractCallback}
- CallbackSet(callbacks...) = new([ cb for cb in callbacks ])
end
+# TODO Deprecate this constructor because of unnecessary compilation!
+CallbackSet(callbacks::AbstractCallback...) = CallbackSet([ cb for cb in callbacks ])
function (cbs::CallbackSet)(u, t, i=0; force=false)
diff --git a/src/evolve.jl b/src/evolve.jl
index 1356ebb4..965c8478 100644
--- a/src/evolve.jl
+++ b/src/evolve.jl
@@ -108,7 +108,7 @@ function step!(evo::Evolution)
@unpack up1, u, t, tend, timestep, stepper!, cb_fullstep = evo
dt = timestep(u, t, 0)
- if isinvalid(dt)
+ if !isfinite(dt)
@info """Invalid timestep proposed at t = $(t): dt = $dt"""
return :timestep_invalid
end
@@ -123,7 +123,7 @@ function step!(evo::Evolution)
return :stepper_failed
end
- if isinvalid(up1)
+ if any(!isfinite, up1)
println()
@warn """NaN or Inf encountered in evolution at t = $(t)!"""
return :invalid_state
@@ -305,17 +305,3 @@ for substeps.
function make_RK_stages(size, n)
return [ zeros(Float64, size) for _ in 1:n ]
end
-
-
-"""
- isinvalid(u)
-
-Return false if any component of `u` is either `NaN` or `Inf`.
-"""
-function isinvalid(u)
- for ui in u
- ( isnan(ui) || isinf(ui) ) && return true
- end
- return false
-end
-isinvalid(u::Real) = isnan(u) || isinf(u)
diff --git a/src/parameters.jl b/src/parameters.jl
index 86afb6d7..61f1391b 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -545,11 +545,11 @@ end
@parameters Output begin
- "Output data after `every_iteration` full time steps"
+ "Output data after `every_iteration` full time steps passed."
every_iteration = 0
@check every_iteration >= 0
- "Output data after `every_dt` simulation time"
+ "Output data after `every_dt` simulation time passed."
every_dt = 0.0
@check every_dt >= 0
@@ -581,16 +581,16 @@ end
@check all(t -> t > 0, interpolate_aligned_ts)
@check issorted(interpolate_aligned_ts)
- "Save timer log after `timer_every_dt_walltime` walltime (in seconds)"
+ "Save timer log after `timer_every_dt_walltime` walltime (in seconds)."
timer_every_dt_walltime = 1
@check timer_every_dt_walltime >= 0
- "Print progress to console after `progress_every_dt` walltime (in seconds)"
+ "Print progress to console after `progress_every_dt` walltime (in seconds)."
progress_every_dt_walltime = 0.2
@check progress_every_dt_walltime >= 0
"""
- List of variables that should be outputted
+ List of variables that should be outputted.
Uses options `every_dt, every_iteration, aligned_ts`.
"""
variables = String[]
@@ -669,7 +669,7 @@ end
@parameters Log begin
- "Enable progress bar logging to stdout"
+ "Enable progress bar logging to stdout."
progress_stdout = true
end
diff --git a/test/IntegrationTests/refs/fv_grhd_bondi_accretion/fv_grhd_bondi_accretion.toml b/test/IntegrationTests/refs/fv_grhd_bondi_accretion/fv_grhd_bondi_accretion.toml
index b6a2c871..cc239b9e 100644
--- a/test/IntegrationTests/refs/fv_grhd_bondi_accretion/fv_grhd_bondi_accretion.toml
+++ b/test/IntegrationTests/refs/fv_grhd_bondi_accretion/fv_grhd_bondi_accretion.toml
@@ -7,6 +7,7 @@ polytrope_gamma = "$(4/3)"
bc = "from_id"
id = "bondi_accretion"
id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"bondi_accretion.h5\"))"
+c2p_freeze_atm_reset = false
[Mesh]
range = [ 1.8, 10.0 ]
diff --git a/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/fv_superbee_grhd_bondi_accretion.toml b/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/fv_superbee_grhd_bondi_accretion.toml
index e1f7e086..6c34c867 100644
--- a/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/fv_superbee_grhd_bondi_accretion.toml
+++ b/test/IntegrationTests/refs/fv_superbee_grhd_bondi_accretion/fv_superbee_grhd_bondi_accretion.toml
@@ -9,7 +9,7 @@ id = "bondi_accretion"
id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"bondi_accretion.h5\"))"
hrsc = "slope_limiter"
slope_limiter_method = "superbee"
-# slope_limiter_method = "none"
+c2p_freeze_atm_reset = false
[Mesh]
range = [ 1.8, 10.0 ]
diff --git a/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml b/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml
index 1a2307c9..37c3ada4 100644
--- a/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml
+++ b/test/IntegrationTests/refs/grhd_bondi_accretion/grhd_bondi_accretion.toml
@@ -7,6 +7,7 @@ polytrope_gamma = "$(4/3)"
bc = "from_id"
id = "bondi_accretion"
id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"bondi_accretion.h5\"))"
+c2p_freeze_atm_reset = false
[Mesh]
range = [ 1.8, 10.0 ]
diff --git a/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml b/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml
index 7f8831e0..591a4a8f 100644
--- a/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml
+++ b/test/IntegrationTests/refs/grhd_tov_spherical1d/grhd_tov_spherical1d.toml
@@ -9,6 +9,7 @@ id = "tov"
id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
atm_factor = 1e-8
formulation = "spherical1d"
+c2p_freeze_atm_reset = false
[Mesh]
range = [ 0.0, 20.0 ]
--
GitLab