From 6fcf974da3b78b296022d55208c5eed4d51f1018 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Tue, 10 Sep 2024 08:25:41 +0000
Subject: [PATCH] GRHD: Implement Bspline approximation
(https://git.tpi.uni-jena.de/dg/dg1d.jl/-/merge_requests/223)
Only implemented for `spherical1d`.
Works somehow, but results are not comparable to FV+TVB method.
---
src/GRHD/GRHD.jl | 27 ++++++++
src/GRHD/callbacks.jl | 59 ++++++++++++++++
src/GRHD/cons2prim.jl | 88 +++++++++++++++++++++++
src/GRHD/rhs.jl | 157 +++++++++++++++++++++++++++++++++++++++++-
src/GRHD/setup.jl | 44 ++++++++++--
5 files changed, 367 insertions(+), 8 deletions(-)
diff --git a/src/GRHD/GRHD.jl b/src/GRHD/GRHD.jl
index a3174a8b..8fdfe091 100644
--- a/src/GRHD/GRHD.jl
+++ b/src/GRHD/GRHD.jl
@@ -315,6 +315,33 @@ dg1d.@parameters GRHD begin
doublecartoon_disable_derivatives = false
@check doublecartoon_disable_derivatives isa Bool
+ """
+ If `true` use a time step that is tailored to the smallest grid spacing.
+ """
+ dt_min_grid_spacing = false
+ @check dt_min_grid_spacing isa Bool
+
+ """
+ If `true` enable the convex hull substep limiter.
+
+ Requires `Mesh.dg_kind = modal_bspline[1,2]`.
+ See also `limiter_limit_log_D, limiter_tci_log_D`.
+ """
+ limiter = false
+ @check limiter isa Bool
+
+ """
+ If `true` apply limiter to `log(D)` instead of `D`.
+ """
+ limiter_limit_log_D = false
+ @check limiter_limit_log_D isa Bool
+
+ """
+ If `true` apply TCI that is used to activate the limiter to `log(D)` instead of `D`.
+ """
+ limiter_tci_log_D = false
+ @check limiter_tci_log_D isa Bool
+
end
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index 66cf7115..10db122a 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -31,6 +31,21 @@ function make_callback(env, P)
end
+function make_substep_callback(env, P)
+ if P.prms.limiter
+ if env.mesh.element.kind ∉ (:modal_bspline1, :modal_bspline2)
+ error("Substep limiter only works with Mesh.dg_kind = modal_bspline[1,2].")
+ end
+ tci = TCI.ModalDecayAverage(env.mesh, 1.0, 3.0)
+ cbfn_limiter(u, t) = callback_limiter(env, P, env.mesh, tci)
+ cb_limiter = FunctionCallback(cbfn_limiter, CallbackTiming(always_on=true))
+ return CallbackSet(cb_limiter)
+ else
+ return nothing
+ end
+end
+
+
#######################################################################
# Equation #
#######################################################################
@@ -543,3 +558,47 @@ function callback_analyze0d(env, P::Project{:spherical1d}, mesh::Mesh1d)
end
end
end
+
+
+#######################################################################
+# limiter #
+#######################################################################
+
+
+callback_limiter(env, P, mesh) = TODO()
+function callback_limiter(env, P::Project{:spherical1d}, mesh::Mesh1d{<:DGElement}, tci)
+ @unpack D, Sr = get_dynamic_variables(mesh)
+ @unpack D_modal, Sr_modal, log_D, c2p_reset_atm = get_static_variables(mesh)
+ @unpack flag = get_cell_variables(mesh)
+
+ broadcast_volume!(impose_atmosphere_spherical1d, P.equation, mesh)
+
+ @unpack limiter_tci_log_D, limiter_limit_log_D = P.prms
+
+ if limiter_tci_log_D || limiter_limit_log_D
+ @. log_D = log(D)
+ end
+
+ K = mesh.tree.dims[1]
+ dg1d.enter(P.prms.workspace) do
+ tmp_flag = dg1d.borrow(P.prms.workspace, length(flag))
+ TCI.compute_indicator!(tmp_flag, limiter_tci_log_D ? log_D : D, P.tci)
+ for (i,f) in enumerate(tmp_flag)
+ flag[i] = f > 0.9
+ end
+ TCI.compute_indicator!(tmp_flag, Sr, P.tci)
+ for (i,f) in enumerate(tmp_flag)
+ flag[i] = max(flag[i], f > 0.9)
+ end
+ tmp_flag .= (flag .+ flag[end:-1:1])./2
+ flag .= tmp_flag
+ end
+
+ dg1d.limit_slope_convex_hull!(limiter_limit_log_D ? log_D : D, flag, P.prms.workspace, mesh)
+ dg1d.limit_slope_convex_hull!(Sr, flag, P.prms.workspace, mesh)
+
+ if limiter_tci_log_D || limiter_limit_log_D
+ @. D = exp(log_D)
+ end
+
+end
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index 07a708d6..ad547ea1 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -395,3 +395,91 @@ 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
+ Ï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
+
+
+@with_signature function impose_atmosphere_spherical1d(eq::Equation)
+ @accepts D, Sr, τ, γrr, c2p_reset_atm_v2
+
+ @unpack atm_density, atm_threshold, cold_eos = eq
+ Ïatm = atm_density
+ patm = cold_eos(Pressure, Density, Ïatm)
+ ϵatm = cold_eos(InternalEnergy, Density, Ïatm)
+ Ïmin = atm_density * atm_threshold
+ γuurr = 1/γrr
+
+ q = Ï„ / D
+ rr = Sr / D
+ r2 = γuurr*rr*rr
+ r = sqrt(r2)
+ k = r / (1.0+q)
+ c2p_reset_atm_v2 = 0.0
+ if D < Ïmin || !(0.0 ≤ k ≤ 1.0)
+ c2p_reset_atm_v2 = 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, Ï„, c2p_reset_atm_v2
+end
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 6a1b8636..ad1cafa6 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -34,8 +34,12 @@ function timestep(mesh, P::Project, hrsc::Maybe{HRSC.AbstractHRSC})
vmax = get_maxspeed(mesh, P.equation)
L, = dg1d.widths(mesh)
K, = mesh.tree.dims
- dL = L/K
- dt = dL / (vmax * dtfactor(mesh))
+ dl = if P.prms.dt_min_grid_spacing
+ min_grid_spacing(mesh)
+ else
+ L/K
+ end
+ dt = dl / (vmax * dtfactor(mesh))
return dt
end
dtfactor(mesh::Mesh1d{FVElement}) = 2
@@ -1197,6 +1201,15 @@ end
function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
+ if mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
+ modal_rhs!(mesh, P, hrsc)
+ else
+ nodal_rhs!(mesh, P, hrsc)
+ end
+end
+
+
+function nodal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
@unpack cache = mesh
@unpack equation = P
@@ -1285,6 +1298,146 @@ function rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
end
+function modal_rhs!(mesh::Mesh1d, P::Project{:spherical1d}, hrsc::Nothing)
+
+ @unpack cache = mesh
+ @unpack equation = P
+
+ @unpack D, Sr, Ï„ = get_dynamic_variables(cache)
+ @unpack rhs_D, rhs_Sr, rhs_Ï„ = get_rhs_variables(cache)
+ @unpack flr_D, flr_Sr, flr_Ï„,
+ max_v, vr, p,
+ src_D, src_Sr, src_Ï„,
+ D_modal, Sr_modal, Ï„_modal,
+ flr_D_modal, flr_Sr_modal, flr_Ï„_modal,
+ src_D_modal, src_Sr_modal, src_Ï„_modal = get_static_variables(cache)
+ @unpack nflr_D, nflr_Sr, nflr_Ï„,
+ bdry_D, bdry_Sr, bdry_Ï„,
+ bdry_max_v, bdry_vr, bdry_p = get_bdry_variables(cache)
+ @unpack flag = get_cell_variables(mesh)
+ @unpack c2p_reset_atm, c2p_reset_atm_v2, log_D = get_static_variables(mesh)
+
+ # evaluate nodal values
+ D_modal .= D
+ Sr_modal .= Sr
+ Ï„_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, P.prms.V_bspline, v_var_modal)
+ end
+ end
+
+ 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)
+
+ dg1d.interpolate_face_data!(mesh, D, bdry_D)
+ dg1d.interpolate_face_data!(mesh, Sr, bdry_Sr)
+ dg1d.interpolate_face_data!(mesh, Ï„, bdry_Ï„)
+ dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+ dg1d.interpolate_face_data!(mesh, vr, bdry_vr)
+ dg1d.interpolate_face_data!(mesh, p, bdry_p)
+
+ broadcast_volume!(flux_source_spherical1d, equation, mesh)
+ impose_symmetry_sources!(P, mesh)
+ broadcast_faces!(llf_spherical1d, equation, mesh)
+ id = P.prms.id
+ bc = P.prms.bc
+ if id == "bondi_accretion" || id == "bondi_accretion_infall"
+ if bc == "from_id"
+ # already accounts for inflow boundary
+ broadcast_bdry!(bdryllf_spherical1d_bondi, equation, mesh)
+ else
+ TODO("unknown boundary conditions for initial data $id: $bc")
+ end
+ elseif id == "tov"
+ if bc == "tov_symmetric_domain"
+ # nothing todo, atmosphering should take care of outer bdrys
+ elseif bc == "tov_reflective_origin"
+ broadcast_bdry!(bdryllf_spherical1d_tov, equation, mesh)
+ else
+ TODO("unknown boundary conditions for initial data $id: $bc")
+ end
+ else
+ TODO("unknown boundary conditions for initial data $id: $bc")
+ end
+
+
+ # evaluate modal coefficients
+ 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)
+ )
+ 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
+
+ if :D ∉ P.prms.freeze_vars
+ compute_rhs_weak_form!(rhs_D, flr_D_modal, src_D_modal, nflr_D, mesh)
+ end
+ if :Sr ∉ P.prms.freeze_vars
+ compute_rhs_weak_form!(rhs_Sr, flr_Sr_modal, src_Sr_modal, nflr_Sr, mesh)
+ end
+ if :τ ∉ P.prms.freeze_vars
+ compute_rhs_weak_form!(rhs_Ï„, flr_Ï„_modal, src_Ï„_modal, nflr_Ï„, mesh)
+ end
+
+ if !P.prms.atm_evolve
+ broadcast_volume!(determine_atmosphere, P.equation, mesh)
+ broadcast_volume!(stop_atmosphere_evolution_spherical1d, P.equation, mesh)
+ end
+
+ if bc == "tov_symmetric_domain"
+ # Need to symmetrize data around origin to avoid asymmetries to blow up there.
+ # 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)
+ # dg1d.enter(P.prms.workspace) do
+ # buf = dg1d.borrow(P.prms.workspace, mesh.element.Npts)
+ # 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 buf .= (v_rhs_D .+ v_rhs_D[end:-1:1])./2
+ # # v_rhs_D .= buf
+ # # @views buf .= (v_rhs_Sr .- v_rhs_Sr[end:-1:1])./2
+ # # v_rhs_Sr .= buf
+ # # @views buf .= (v_rhs_Ï„ .+ v_rhs_Ï„[end:-1:1])./2
+ # # v_rhs_Ï„ .= buf
+ # # @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
+ # end
+ K = mesh.tree.dims[1]
+ dg1d.enter(P.prms.workspace) do
+ buf = dg1d.borrow(P.prms.workspace, mesh.element.Npts*K)
+ @views begin
+ @. buf = (rhs_D + rhs_D[end:-1:1])/2
+ rhs_D .= buf
+ @. buf = (rhs_Sr - rhs_Sr[end:-1:1])/2
+ rhs_Sr .= buf
+ @. buf = (rhs_Ï„ + rhs_Ï„[end:-1:1])/2
+ rhs_Ï„ .= buf
+ end
+ end
+ end
+
+ D .= D_modal
+ Sr .= Sr_modal
+ Ï„ .= Ï„_modal
+
+ return
+end
+
+
function rhs!(mesh::Mesh1d, P::Project{:rescaled_spherical1d}, hrsc::Nothing)
@unpack cache = mesh
diff --git a/src/GRHD/setup.jl b/src/GRHD/setup.jl
index 3c76c8bb..7638364e 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -40,15 +40,35 @@ function Project(env::Environment, prms)
c2p_enforce_causal_atm = Bool(prms["GRHD"]["c2p_enforce_causal_atm"])
atm_evolve = Bool(prms["GRHD"]["atm_evolve"])
hrsc = Symbol(prms["GRHD"]["hrsc"])
+ limiter = prms["GRHD"]["limiter"]
+ limiter_limit_log_D = prms["GRHD"]["limiter_limit_log_D"]
+ limiter_tci_log_D = prms["GRHD"]["limiter_tci_log_D"]
freeze_vars = Symbol.(prms["GRHD"]["freeze_vars"])
fv_numerical_flux = prms["GRHD"]["fv_numerical_flux"]
id = prms["GRHD"]["id"]
bc = prms["GRHD"]["bc"]
+ dt_min_grid_spacing = prms["GRHD"]["dt_min_grid_spacing"]
+ workspace = dg1d.Workspace()
+ 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
+ V_bspline, invV_bspline = nothing, nothing
+ end
+ else
+ V_bspline, invV_bspline = nothing, nothing
+ end
fixedprms = (; av_regularization=:covariant, id_smooth=true,
bernstein, slope_limiter_method, slope_limiter_tvb_M,
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, dt_min_grid_spacing,
+ limiter, limiter_limit_log_D, limiter_tci_log_D,
id, bc, hrsc, fv_numerical_flux)
# construct Project
@@ -77,7 +97,9 @@ function Project(env::Environment, prms)
initialdata!(env, P)
projectcb = make_callback(env, P)
- append!(env.callbacks, CallbackSet(projectcb.callbacks))
+ append!(env.callbacks, projectcb)
+ substep_projectcb = make_substep_callback(env, P)
+ append!(env.callbacks_substep, substep_projectcb)
return P, nothing
end
@@ -169,11 +191,13 @@ function register_evolution!(mesh, P::Project{:spherical1d})
rhs_variablenames = (:rhs_D, :rhs_Sr, :rhs_Ï„), # rhs of conservatives
static_variablenames = (:flr_D, :flr_Sr, :flr_Ï„, # conservatives' fluxes
:src_D, :src_Sr, :src_Ï„, # sources
- :Ï, :vr, :p, :ϵ, # primitives
+ :Ï, :vr, :p, :ϵ, # primitives
+ :raw_Ï, :raw_vr, :raw_p, :raw_ϵ, # unlimited primitives
:max_v, :r, # maximum charactersitic speed, radius
:init_D, :init_Sr, :init_Ï„, :init_max_v, :init_vr, :init_p,
:init_Ï, :init_ϵ,
:E, :Em1, :flr_E, :flr_Em1, :EP,
+ :log_D,
:ÏhW2),
bdry_variablenames = (:bdry_D, :bdry_Sr, :bdry_Ï„,
:bdry_Ï, :bdry_ϵ, :bdry_vr, :bdry_p,
@@ -184,6 +208,12 @@ function register_evolution!(mesh, P::Project{:spherical1d})
)
@unpack r, x = get_static_variables(mesh)
r .= x
+ register_variables!(mesh,
+ static_variablenames = (:D_modal, :Sr_modal, :Ï„_modal,
+ :flr_D_modal, :flr_Sr_modal, :flr_Ï„_modal,
+ :src_D_modal, :src_Sr_modal, :src_Ï„_modal,
+ :c2p_failed)
+ )
end
function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
hrsc = P.prmsdb["GRHD"]["hrsc"]
@@ -480,8 +510,10 @@ end
function register_analysis!(mesh::Mesh1d, P)
register_variables!(mesh,
static_variablenames = (:c2p_reset_ϵ, :c2p_reset_atm, :c2p_limit_vr,
- :c2p_freeze_atm, :c2p_init_admissible, :v),
- bdry_variablenames = (:bdry_c2p_reset_atm,)
+ :c2p_freeze_atm, :c2p_init_admissible, :v,
+ :c2p_reset_atm_v2),
+ bdry_variablenames = (:bdry_c2p_reset_atm,),
+ cell_variablenames = (:flag,),
)
if "Mtot" in P.prmsdb["GRHD"]["variables0d_analyze"]
register_variables!(mesh, global_variablenames = (:Mtot,))
@@ -576,8 +608,8 @@ register_tci!(mesh, P, tci::TCI.AbstractTCI) = TODO(tci)
function register_tci!(mesh, P::Project{:spherical1d}, tci::TCI.AbstractTCI)
- TODO(tci)
- register_variables!(mesh, cell_variablenames = (:flag,:sD_flag,:sSr_flag,:stau_flag))
+ return
+ register_variables!(mesh, cell_variablenames = (:flag,))
end
function register_tci!(mesh, P::Project{:spherical1d}, tci::TCI.EntropyProduction)
--
GitLab