From f9d5fc7326ed73cb06c4005e572b46f1b8d255bf Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Thu, 29 Aug 2024 12:18:48 +0200
Subject: [PATCH] bspline1_grhd_tov_spherical1d is doing something
---
bspline1_grhd_tov_spherical1d.toml | 37 +++---
src/GRHD/callbacks.jl | 97 +++++++++++-----
src/GRHD/cons2prim.jl | 53 +++++++++
src/GRHD/rhs.jl | 180 ++++++++++++++++-------------
src/GRHD/setup.jl | 11 +-
src/callbacks.jl | 2 +-
src/evolve.jl | 11 +-
7 files changed, 254 insertions(+), 137 deletions(-)
diff --git a/bspline1_grhd_tov_spherical1d.toml b/bspline1_grhd_tov_spherical1d.toml
index 15fa6eda..5aaaacdd 100644
--- a/bspline1_grhd_tov_spherical1d.toml
+++ b/bspline1_grhd_tov_spherical1d.toml
@@ -7,42 +7,48 @@ eos = "polytrope"
cfl = 0.2
tend = 1000.0
# tend = 290.0
-method = "ssprk3"
+method = "rk4"
+# method = "ssprk4"
[Output]
+
# # variables0d = ["Mtot"]
-# aligned_ts = "$(collect(range(1.0,2000.0,step=1.0)))"
+
aligned_ts = "$(collect(range(1.0,2000.0,step=1.0)))"
-# every_iteration = 1
+# # every_iteration = 1
variables1d = ["D", "Sr", "Ï„",
"D_modal", "Sr_modal", "Ï„_modal",
"flr_D_modal", "flr_Sr_modal", "flr_Ï„_modal",
"flr_D", "flr_Sr", "flr_Ï„",
# "src_D_modal", "src_Sr_modal", "src_Ï„_modal",
"rhs_D", "rhs_Sr", "rhs_Ï„",
- "Ï", "p", "vr", "ÏhW2",
+ "Ï", "p", "vr",
+ "raw_Ï", "raw_p", "raw_vr",
# "D_err", "Ï_err",
"c2p_reset_atm", "c2p_freeze_atm", "c2p_failed"
]
-#
+
+
# interpolate_nodes = """$(collect(range(-20.0,20.0,step=0.01)))"""
# interpolate_variables = ["D", "Sr", "Ï„",
# "flr_D", "flr_Sr", "flr_Ï„",
# "src_D", "src_Sr", "src_Ï„",
# "rhs_D", "rhs_Sr", "rhs_Ï„"]
-# interpolate_every_iteration = 1
+# interpolate_aligned_ts = "$(collect(range(1.0,2000.0,step=1.0)))"
# substep_every_iteration = 1
+# # substep_every_dt = 0.01
# substep_variables = ["D", "Sr", "Ï„",
# "D_modal", "Sr_modal", "Ï„_modal",
-# "flr_D", "flr_Sr", "flr_Ï„",
-# "flr_D_modal", "flr_Sr_modal", "flr_Ï„_modal",
-# "src_D", "src_Sr", "src_Ï„",
-# "src_D_modal", "src_Sr_modal", "src_Ï„_modal",
-# "c2p_reset_atm", "c2p_freeze_atm",
-# "Ï", "p", "vr", "ÏhW2",
+# # "flr_D", "flr_Sr", "flr_Ï„",
+# # "flr_D_modal", "flr_Sr_modal", "flr_Ï„_modal",
+# # "src_D", "src_Sr", "src_Ï„",
+# # "src_D_modal", "src_Sr_modal", "src_Ï„_modal",
+# "c2p_reset_atm", "c2p_freeze_atm", "flag",
+# "Ï", "p", "vr", "ϵ", "ÏhW2",
+# "raw_Ï", "raw_p", "raw_vr", "raw_ϵ",
# "rhs_D", "rhs_Sr", "rhs_Ï„"
# ]
@@ -54,7 +60,7 @@ analyze_error_reference_solution = "from_id"
atm_threshold_factor = 100.001
# atm_threshold_factor = 1.001
id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
-atm_factor = 1.0e-10
+atm_factor = 1.0e-8
id = "tov"
formulation = "spherical1d"
bc = "tov_symmetric_domain"
@@ -63,13 +69,12 @@ atm_evolve = false
c2p_enforce_causal_atm = true
c2p_enforce_atm = false
c2p_set_atmosphere_on_failure = true
-c2p_threshold_velocity_limiter = 1e-8
-# perturbation_rho = true
+c2p_threshold_velocity_limiter = 1e-6
[Mesh]
periodic = false
range = [-20.0, 20.0]
k = 21
basis = "lgl"
-n = 5
+n = 11
dg_kind = "modal_bspline1"
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index d3229459..5148c88c 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -8,7 +8,7 @@ function make_callback(env, P::Project)
cbfn_hrsc(u, t) = callback_hrsc(u, t, env, P)
cb_hrsc = FunctionCallback(cbfn_hrsc,
CallbackTiming(every_iteration=1,every_dt=0))
- callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
+ callbackset = CallbackSet(cb_equation, #=cb_tci,=# cb_hrsc)
if !isempty(P.prmsdb["GRHD"]["variables0d_analyze_error"])
cbfn_analyzed0d_err(u, t) = callback_analyze0d_error(env, P, env.mesh)
cb = FunctionCallback(cbfn_analyzed0d_err,
@@ -32,10 +32,10 @@ end
function make_callbacks_substeps(env, P::Project)
- cbfn_limiter(u, t) = callback_limiter(env, P)
- cb_limiter = FunctionCallback(cbfn_limiter,
- CallbackTiming(every_iteration=1, every_dt=0))
- CallbackSet(cb_limiter)
+ # cbfn_limiter(u, t) = callback_limiter(env, P)
+ # cb_limiter = FunctionCallback(cbfn_limiter,
+ # CallbackTiming(every_iteration=1, every_dt=0))
+ # CallbackSet(cb_limiter)
end
@@ -74,23 +74,53 @@ function compute_atmosphere_mask(env, P, mesh::Mesh1d{DGElement})
@unpack D, Sr, Ï„ = get_dynamic_variables(mesh)
@unpack D_modal, Sr_modal, Ï„_modal = get_static_variables(mesh)
- D_modal .= D
- Sr_modal .= Sr
- Ï„_modal .= Ï„
- for (var,var_modal) in ( (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal) )
- for (k,(v_var, v_var_modal)) in enumerate(zip(eachcell(mesh,var),eachcell(mesh,var_modal)))
- mul!(v_var, P.prms.V_bspline, v_var_modal)
- end
- end
-
if formulation(P) === :valencia1d
broadcast_volume!(cons2prim_valencia1d, P.equation, env.mesh)
elseif formulation(P) === :spherical1d
-
broadcast_volume!(cons2prim_spherical1d, P.equation, env.mesh)
+
+ # slope limit modal data
+ @unpack p, Ï, ϵ, vr = get_static_variables(mesh)
+ @unpack raw_p, raw_Ï, raw_ϵ, raw_vr = get_static_variables(mesh)
+ @unpack c2p_reset_atm, c2p_freeze_atm = get_static_variables(mesh)
+ raw_p .= p
+ raw_Ï .= Ï
+ raw_vr .= vr
+ raw_ϵ .= ϵ
+ @unpack flag = get_cell_variables(mesh)
+ # TCI.compute_indicator!(flag, vr, P.prms.tci)
+ for (k,f) in enumerate(flag)
+ cidxs = cellindices(mesh, k)
+ rst = c2p_reset_atm[cidxs]
+ frz = c2p_freeze_atm[cidxs]
+ if all(r -> r>0, rst) || all(r -> r ≈ 0.0, rst)
+ flag[k] = 0.0
+ else
+ flag[k] = Float64( any(r->r>0, rst) && any(f->f≈0.0, frz) )
+ end
+ end
+ # for var in (Ï, ϵ, vr)
+ for var in (vr,)
+ limit_slope!(var, flag, P, mesh)
+ end
+ @unpack eos = P.equation
+ p .= eos.(Pressure, Density, Ï, InternalEnergy, ϵ)
+ # impose atmosphere
+ broadcast_volume!(prim2cons_spherical1d, P.equation, mesh)
+
+
+ # D_modal .= D
+ # Sr_modal .= Sr
+ # Ï„_modal .= Ï„
+ # for (var,var_modal) in ( (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal) )
+ # for (k,(v_var, v_var_modal)) in enumerate(zip(eachcell(mesh,var),eachcell(mesh,var_modal)))
+ # mul!(v_var, P.prms.V_bspline, v_var_modal)
+ # end
+ # end
+
if P.prms.atm_equalize_on_interface
broadcast_volume!(determine_atmosphere, P.equation, env.mesh)
@unpack c2p_reset_atm = get_static_variables(env.mesh)
@@ -122,14 +152,14 @@ function compute_atmosphere_mask(env, P, mesh::Mesh1d{DGElement})
- for (var,var_modal) in ((D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal))
- for (v_var,v_var_modal) in zip(eachcell(mesh,var),eachcell(mesh,var_modal))
- mul!(v_var_modal, P.prms.invV_bspline, v_var)
- end
- end
- D .= D_modal
- Sr .= Sr_modal
- Ï„ .= Ï„_modal
+ # for (var,var_modal) in ((D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal))
+ # for (v_var,v_var_modal) in zip(eachcell(mesh,var),eachcell(mesh,var_modal))
+ # mul!(v_var_modal, P.prms.invV_bspline, v_var)
+ # end
+ # end
+ # D .= D_modal
+ # Sr .= Sr_modal
+ # Ï„ .= Ï„_modal
end
@@ -504,8 +534,8 @@ function compute_tci(P, eq, mesh,
TCI.ModalDecayAverage, TCI.ModalDecayHighest})
var = tci_indicator_variable(P, mesh)
@unpack flag = get_cell_variables(mesh)
- TCI.compute_indicator!(flag, log10.(abs.(var)), tci)
- smooth_flag!(flag, mesh)
+ TCI.compute_indicator!(flag, var, tci)
+ # smooth_flag!(flag, mesh)
end
function tci_indicator_variable(::Project{:doublecartoon}, mesh)
@@ -594,18 +624,22 @@ callback_limiter(env, P::Project) = nothing
function callback_limiter(env, P::Project{:spherical1d})
+ return
@unpack mesh = env
@unpack D, Sr, Ï„ = get_dynamic_variables(mesh)
+ # Main.verbose && println("RUNNING LIMITER!")
+ # L = layout(mesh)
+ # display(dg1d.reshape(D, L))
for var in (D, Sr, Ï„)
+ # for var in (D,)
limit_slope!(var, P, mesh)
end
+ # display(dg1d.reshape(D, L))
+ # readline()
end
-# TODO Add mask!
-function limit_slope!(var::AbstractArray, P::Project, mesh::Mesh1d)
-
- # TODO Test with mask
+function limit_slope!(var::AbstractArray, mask::AbstractArray, P::Project, mesh::Mesh1d)
dg1d.enter(P.prms.workspace) do
@@ -622,7 +656,7 @@ function limit_slope!(var::AbstractArray, P::Project, mesh::Mesh1d)
A = reshape(mat_A,(Npts-1,Npts-1))
for k in 1:K
- if k == 1 || k == K
+ if k == 1 || k == K || mask[k] ≈ 0.0
# skip bdry cells, but backup data because we restore the whole vector at once later
idxs = cellindices(mesh, k)
view(tmp_var, idxs) .= view(var, idxs)
@@ -641,7 +675,6 @@ function limit_slope!(var::AbstractArray, P::Project, mesh::Mesh1d)
v1, v2 = vvar_rhs[1], vvar_rhs[2]
x1, x2 = vx_rhs[1], vx_rhs[2]
avgs[end] = (v1+v2)/2; δs[end] = (v2-v1)/(x2-x1)
- L = 0.0
for i in 1:Npts-1
v1, v2 = vvar[i], vvar[i+1]
x1, x2 = vx[i], vx[i+1]
@@ -653,7 +686,7 @@ function limit_slope!(var::AbstractArray, P::Project, mesh::Mesh1d)
# lim_δs[i-1] = TCI.minmod(δs[i-1],δs[i],δs[i+1])
x1, x2 = vx[i-1], vx[i]
dx = x2-x1
- M = 1.0
+ M = 0.01
δl, δc, δr = δs[i-1], δs[i], δs[i+1]
lim_δs[i-1] = TCI.minmod_M(dx, M, δc, δl, δr)
end
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index 07a708d6..0bd3a7b6 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -191,6 +191,7 @@ end
k = r / (1.0+q)
if !(0.0 ≤ k ≤ 1.0) # admissibility constraint from arXiv:1306.4953, Appendix
c2p_init_admissible = 1.0
+ @goto impose_atmosphere
# rescale Si such that r = 1+q
rescale = 0.99*(1+q)/r
Sx *= rescale
@@ -395,3 +396,55 @@ end
c2p_reset_atm = Float64(D<Ïmin)
@returns c2p_reset_atm
end
+
+
+#######################################################################
+# primitives to conservatives #
+#######################################################################
+
+
+@with_signature function prim2cons_spherical1d(eq::Equation)
+ @accepts Ï, p, ϵ, vr, γrr
+
+ @unpack atm_density, atm_threshold, cold_eos = eq
+ # atmosphere values, assuming cold EoS
+ Ïatm = atm_density
+ patm = cold_eos(Pressure, Density, Ïatm)
+ ϵatm = cold_eos(InternalEnergy, Density, Ïatm)
+ Ïmin = atm_density * atm_threshold
+
+ γuurr = 1/γrr
+ v2 = γuurr*vr*vr
+ h = 1.0+ϵ+p/Ï
+ W2 = 1.0/(1.0-v2)
+ W = sqrt(W2)
+ ÏhW2 = Ï*h*W2
+ D = W*Ï
+ Sr = ÏhW2*vr
+ Ï„ = ÏhW2 - p - D
+
+ q = Ï„ / D
+ rr = Sr / D
+ r2 = γuurr*rr*rr
+ r = sqrt(r2)
+ k = r / (1.0+q)
+ if D < Ïmin || !(0.0 ≤ k ≤ 1.0)
+ Ï = Ïatm
+ vr = 0.0
+ W = 1.0
+ p = patm
+ ϵ = ϵatm
+
+ v2 = γuurr*vr*vr
+ h = 1.0+ϵ+p/Ï
+ W2 = 1.0/(1.0-v2)
+ W = sqrt(W2)
+ ÏhW2 = Ï*h*W2
+ D = W*Ï
+ Sr = ÏhW2*vr
+ Ï„ = ÏhW2 - p - D
+ end
+
+
+ @returns D, Sr, Ï„, ÏhW2, Ï, p, ϵ, vr
+end
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index f432b661..add3fab2 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -499,6 +499,7 @@ function step_fv_slope_limiter!(mesh, P::Project{:spherical1d}, state)
rhs_D[1] = rhs_Sr[1] = rhs_Ï„[1] = rhs_D[end] = rhs_Sr[end] = rhs_Ï„[end] = 0
end
+ # TODO Is this needed for the FV evolutions?
if !P.prms.atm_evolve
broadcast_volume!(determine_atmosphere, P.equation, mesh)
broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
@@ -1337,76 +1338,62 @@ function modal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
src_D_modal, src_Sr_modal, src_Ï„_modal,
rhs_D_modal, rhs_Sr_modal, rhs_Ï„_modal = get_static_variables(cache)
- D_modal .= D
- Sr_modal .= Sr
- Ï„_modal .= Ï„
-
- # Ïatm = P.equation.atm_density
- # Ïmin = P.equation.atm_density * P.equation.atm_threshold
- # ϵatm = P.equation.cold_eos(InternalEnergy, Density, Ïatm)
- # patm = P.equation.cold_eos(Pressure, Density, Ïatm)
- # hatm = 1.0+ϵatm+patm/Ïatm
- # ÏhW2atm = Ïatm*hatm
- # Datm = Ïatm
- # Sratm = 0.0
- # Ï„atm = ÏhW2atm - patm - Datm
- # for k in 1:mesh.tree.dims[1]
- # cidxs = cellindices(mesh, k)
- # v_D_modal = view(D_modal, cidxs)
- # v_Sr_modal = view(Sr_modal, cidxs)
- # v_Ï„_modal = view(Ï„_modal, cidxs)
- # for i in 1:length(cidxs)
- # if v_D_modal[i] < Ïmin
- # v_D_modal[i] = Datm
- # v_Sr_modal[i] = Sratm
- # v_Ï„_modal[i] = Ï„atm
- # end
- # end
- # end
-
- # n_resamples = 0
- # @label resample_nodal
- # n_resamples += 1
- # if n_resamples > 3
- # error("too many attempts on resampling!")
- # end
- for (var,var_modal) in ( (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal) )
- for (k,(v_var, v_var_modal)) in enumerate(zip(eachcell(mesh,var),eachcell(mesh,var_modal)))
- mul!(v_var, P.prms.V_bspline, v_var_modal)
- end
- end
-
+ # c2p with modal data
if P.prms.c2p_enforce_causal_atm || P.prms.c2p_enforce_atm
broadcast_volume!(cons2prim_spherical1d_freeze_flags, equation, mesh)
else
broadcast_volume!(cons2prim_spherical1d, equation, mesh)
end
- broadcast_volume!(maxspeed_spherical1d, equation, mesh)
- # @unpack c2p_failed = get_static_variables(mesh)
- # # if any(>(0.0), c2p_failed)
- # # TODO("handle c2p failure")
- # # end
- # Ïatm = P.equation.atm_density
- # c2p_success = true
- # for (k,v_c2p_failed) in enumerate(eachcell(mesh,c2p_failed))
- # any(>(0.0), c2p_failed) || continue
- # c2p_success = false
- # cidxs = cellindices(mesh, k)
- # v_D = view(D, cidxs)
- # v_Sr = view(Sr, cidxs)
- # v_Ï„ = view(Ï„, cidxs)
- # v_D_modal = view(D_modal, cidxs)
- # v_Sr_modal = view(Sr_modal, cidxs)
- # v_Ï„_modal = view(Ï„_modal, cidxs)
- # for (i,f) in enumerate(v_c2p_failed)
- # f > 0.0 || continue
- # v_D_modal[i] = Ïatm
- # v_Sr_modal[i] = 0.0
- # v_Ï„_modal[i] = 0.0
+ # slope limit modal data
+ @unpack p, Ï, ϵ, vr = get_static_variables(mesh)
+ @unpack raw_p, raw_Ï, raw_ϵ, raw_vr = get_static_variables(mesh)
+ @unpack c2p_reset_atm, c2p_freeze_atm = get_static_variables(mesh)
+ raw_p .= p
+ raw_Ï .= Ï
+ raw_vr .= vr
+ raw_ϵ .= ϵ
+ @unpack flag = get_cell_variables(mesh)
+ for (k,f) in enumerate(flag)
+ cidxs = cellindices(mesh, k)
+ rst = c2p_reset_atm[cidxs]
+ frz = c2p_freeze_atm[cidxs]
+ if all(r -> r>0, rst) || all(r -> r ≈ 0.0, rst)
+ flag[k] = 0.0
+ else
+ flag[k] = Float64( any(r->r>0, rst) && any(f->f≈0.0, frz) )
+ end
+ end
+ for var in (vr,)
+ limit_slope!(var, flag, P, mesh)
+ end
+ @unpack eos = P.equation
+ p .= eos.(Pressure, Density, Ï, InternalEnergy, ϵ)
+ # impose atmosphere
+ broadcast_volume!(prim2cons_spherical1d, equation, mesh)
+
+ # slope limit modal data near origin
+ K = mesh.tree.dims[1]
+ fill!(flag, 0.0)
+ flag[K÷2+1] = 1.0
+ for var in (vr,)
+ limit_slope!(var, flag, P, mesh)
+ end
+ broadcast_volume!(prim2cons_spherical1d, equation, mesh)
+
+ # # only now sample nodal data and compute residuals
+ # D_modal .= D
+ # Sr_modal .= Sr
+ # Ï„_modal .= Ï„
+
+
+ # for (var,var_modal) in ( (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal) )
+ # for (k,(v_var, v_var_modal)) in enumerate(zip(eachcell(mesh,var),eachcell(mesh,var_modal)))
+ # mul!(v_var, P.prms.V_bspline, v_var_modal)
# end
# end
- # !c2p_success && @goto resample_nodal
+
+ broadcast_volume!(maxspeed_spherical1d, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
dg1d.interpolate_face_data!(mesh, Sr, bdry_Sr)
@@ -1442,7 +1429,7 @@ function modal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
for (var,var_modal) in ( (flr_D,flr_D_modal), (flr_Sr,flr_Sr_modal), (flr_Ï„,flr_Ï„_modal),
(src_D,src_D_modal), (src_Sr,src_Sr_modal), (src_Ï„,src_Ï„_modal),
- (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal)
+ # (D,D_modal), (Sr,Sr_modal), (Ï„,Ï„_modal)
)
for (v_var,v_var_modal) in zip(eachcell(mesh,var),eachcell(mesh,var_modal))
mul!(v_var_modal, P.prms.invV_bspline, v_var)
@@ -1469,28 +1456,27 @@ function modal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
# This is needed, because near the surface the solution becomes asymmetric,
# and those errors are amplified when propagating inwards.
# Not sure if this is caused by the LDG/AV method or by the DG method itself.
- # korigin = find_cell_index(0.0, mesh)
- # for (k,(v_rhs_D, v_rhs_Sr, v_rhs_Ï„)) in enumerate(zip(eachcell(mesh,rhs_D),
- # eachcell(mesh,rhs_Sr),
- # eachcell(mesh,rhs_Ï„)))
- # if korigin == k
- # # @views v_rhs_D .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
- # # @views v_rhs_D .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
- # # @views v_rhs_Sr .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
- # # @views v_rhs_Ï„ .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
- # break
- # end
- # end
- @views begin
- @. rhs_D = (rhs_D + rhs_D[end:-1:1])/2
- @. rhs_Sr = (rhs_Sr - rhs_Sr[end:-1:1])/2
- @. rhs_Ï„ = (rhs_Ï„ + rhs_Ï„[end:-1:1])/2
+ korigin = find_cell_index(0.0, mesh)
+ for (k,(v_rhs_D, v_rhs_Sr, v_rhs_Ï„)) in enumerate(zip(eachcell(mesh,rhs_D),
+ eachcell(mesh,rhs_Sr),
+ eachcell(mesh,rhs_Ï„)))
+ if korigin == k
+ @views v_rhs_D .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
+ @views v_rhs_Sr .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
+ @views v_rhs_Ï„ .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
+ break
+ end
end
+ # @views begin
+ # @. rhs_D = (rhs_D + rhs_D[end:-1:1])/2
+ # @. rhs_Sr = (rhs_Sr - rhs_Sr[end:-1:1])/2
+ # @. rhs_Ï„ = (rhs_Ï„ + rhs_Ï„[end:-1:1])/2
+ # end
end
- D .= D_modal
- Sr .= Sr_modal
- Ï„ .= Ï„_modal
+ # D .= D_modal
+ # Sr .= Sr_modal
+ # Ï„ .= Ï„_modal
return
end
@@ -1890,6 +1876,8 @@ end
function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialViscosity)
+ TODO()
+
@unpack cache = mesh
@unpack equation = P
@@ -1913,6 +1901,34 @@ function rhs!(mesh, P::Project{:spherical1d}, hrsc::HRSC.AbstractArtificialVisco
else
broadcast_volume!(cons2prim_spherical1d, equation, mesh)
end
+
+ @unpack p, Ï, ϵ, vr = get_static_variables(mesh)
+ @unpack raw_p, raw_Ï, raw_ϵ, raw_vr = get_static_variables(mesh)
+ raw_p .= p
+ raw_Ï .= Ï
+ raw_vr .= vr
+ raw_ϵ .= ϵ
+ @unpack flag = get_cell_variables(mesh)
+ TCI.compute_indicator!(flag, vr, P.prms.tci)
+ for (k,f) in enumerate(flag)
+ cidxs = cellindices(mesh, k)
+ rst = c2p_reset_atm[cidxs]
+ if all(r -> r>0, rst) || all(r -> r ≈ 0.0, rst)
+ flag[k] = 0.0
+ else
+ flag[k] = f > 0.0 ? 1.0 : 0.0
+ end
+ end
+ for var in (Ï, ϵ, vr)
+ HRSC.reconstruct!(var, flag, P.prms.bernstein, isperiodic=isperiodic(mesh))
+ end
+ # for var in (Ï, ϵ, vr)
+ # limit_slope!(var, P, mesh)
+ # end
+ @unpack eos = P.equation
+ p .= eos.(Pressure, Density, Ï, InternalEnergy, ϵ)
+ broadcast_volume!(prim2cons_spherical1d, equation, mesh)
+
broadcast_volume!(maxspeed_spherical1d, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index c0e99ee0..90d3f199 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -45,15 +45,17 @@ function Project(env::Environment, prms)
id = prms["GRHD"]["id"]
bc = prms["GRHD"]["bc"]
workspace = dg1d.Workspace()
+ tci = TCI.ModalDecayAverage(mesh, 1.0, 3.0)
if mesh.element isa DGElement
if mesh.element.kind === :modal_bspline2
V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs2)
+ invV_bspline = pinv(V_bspline)
elseif mesh.element.kind === :modal_bspline1
V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs1)
+ invV_bspline = pinv(V_bspline)
else
- TODO()
+ V_bspline, invV_bspline = nothing, nothing
end
- invV_bspline = pinv(V_bspline)
else
V_bspline, invV_bspline = nothing, nothing
end
@@ -62,7 +64,7 @@ function Project(env::Environment, prms)
c2p_dynamic_atm, atm_evolve, atm_equalize_on_interface,
c2p_set_atmosphere_on_failure, c2p_enforce_causal_atm, c2p_enforce_atm,
av_drag, av_sensor_abslog_D, problem, freeze_vars,
- V_bspline, invV_bspline, workspace,
+ V_bspline, invV_bspline, workspace, tci,
id, bc, hrsc, fv_numerical_flux)
# construct Project
@@ -506,7 +508,8 @@ function register_analysis!(mesh::Mesh1d, P)
register_variables!(mesh,
static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
:c2p_freeze_atm, :c2p_init_admissible, :v),
- bdry_variablenames = (:bdry_c2p_reset_atm,)
+ bdry_variablenames = (:bdry_c2p_reset_atm,),
+ cell_variablenames = (:flag,),
)
if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
register_variables!(mesh, global_variablenames = (:Mtot,))
diff --git a/src/callbacks.jl b/src/callbacks.jl
index f1c99d86..c9b8aad2 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -844,7 +844,7 @@ mutable struct ProgressCallback <: AbstractCallback
prog_buffer = IOBuffer(; sizehint=2048)
prog = Progress(100, # %
dt=0.0, # always update
- barglyphs=BarGlyphs('|','â–ˆ', ['â–' ,'â–‚' ,'â–ƒ' ,'â–„' ,'â–…' ,'â–†', 'â–‡'],' ','|',),
+ # barglyphs=BarGlyphs('|','â–ˆ', ['â–' ,'â–‚' ,'â–ƒ' ,'â–„' ,'â–…' ,'â–†', 'â–‡'],' ','|',),
showspeed=true, output=prog_buffer)
# accumulate report here
report_buffer = IOBuffer(; sizehint=2048)
diff --git a/src/evolve.jl b/src/evolve.jl
index 7aff7aea..d2e12df0 100644
--- a/src/evolve.jl
+++ b/src/evolve.jl
@@ -257,8 +257,11 @@ Stepper for Explicit Runge Kutta scheme.
function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::TimeStepAlgorithm, callback_substep)
@unpack nstages, a, b, c = coeffs
+ # Main.verbose && println("TAKING A STEP")
+
for i = 1:nstages
ti = t + c[i] * dt
+ # Main.verbose && println("SHOULD BE RUNNING RHS")
# assumes that a[1,:] are all zeros
F(tmp_k[i], i == 1 ? u0 : u!, ti)
i == nstages && break
@@ -266,7 +269,8 @@ function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::TimeStepAlgorithm, callback_
for j = 1:i
@. u! += a[i+1, j] * dt * tmp_k[j]
end
- success = callback_substep(u!, ti, i+1) # +1 because u! is value for next stage
+ # Main.verbose && println("SHOULD BE RUNNING LIMITER")
+ success = callback_substep(u!, ti, i+1, force=true) # +1 because u! is value for next stage
!success && return false
end
@@ -274,7 +278,10 @@ function ode_step!(u!, F, u0, t, dt, tmp_k, coeffs::TimeStepAlgorithm, callback_
for i = 1:nstages
@. u! += b[i] * dt * tmp_k[i]
end
- success = callback_substep(u!, t+dt, nstages)
+ # Main.verbose && println("SHOULD BE RUNNING LIMITER")
+ success = callback_substep(u!, t+dt, nstages, force=true)
+
+ # Main.verbose && println('='^10)
return success
end
--
GitLab