From 3fac171f1107590b04ac828921f99a5c6f4a5b75 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Fri, 6 Sep 2024 15:56:49 +0200
Subject: [PATCH] wip
---
src/GRHD/callbacks.jl | 147 ++++++++++++++++++++++++++++++
src/GRHD/cons2prim.jl | 87 +++++++++++++++---
src/GRHD/rhs.jl | 206 +++++++++++++++++++++++++++++++++++++++++-
src/GRHD/setup.jl | 28 +++++-
src/bspline.jl | 30 ++++++
5 files changed, 479 insertions(+), 19 deletions(-)
diff --git a/src/GRHD/callbacks.jl b/src/GRHD/callbacks.jl
index 66cf7115..56348160 100644
--- a/src/GRHD/callbacks.jl
+++ b/src/GRHD/callbacks.jl
@@ -61,6 +61,8 @@ function compute_atmosphere_mask(env, P, mesh::Mesh1d{DGElement})
enforce = P.prms.c2p_enforce_atm
enforce_causal = P.prms.c2p_enforce_causal_atm
(!enforce && !enforce_causal) && return
+
+
if formulation(P) === :valencia1d
broadcast_volume!(cons2prim_valencia1d, P.equation, env.mesh)
elseif formulation(P) === :spherical1d
@@ -71,6 +73,46 @@ function compute_atmosphere_mask(env, P, mesh::Mesh1d{DGElement})
broadcast_volume!(cons2prim_doublecartoon, P.equation, env.mesh)
elseif formulation(P) === :cartoon
broadcast_volume!(cons2prim_cartoon, 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)
+ 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) )
+ flag[k] = 1.0
+ end
+ end
+ K = mesh.tree.dims[1]
+ flag[K÷2+1] = 1.0
+ for var in (vr,)
+ # limit_slope!(var, flag, P, mesh)
+ end
+ @unpack eos = P.equation
+ # p .= eos.(Pressure, Density, Ï, InternalEnergy, ϵ)
+ # impose atmosphere on modal data
+ # broadcast_volume!(prim2cons_spherical1d, P.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)
+
if P.prms.equalize_atmosphere_spherical1d
broadcast_volume!(determine_atmosphere, P.equation, env.mesh)
@unpack c2p_reset_atm = get_static_variables(env.mesh)
@@ -543,3 +585,108 @@ function callback_analyze0d(env, P::Project{:spherical1d}, mesh::Mesh1d)
end
end
end
+
+
+#######################################################################
+# Limiter #
+#######################################################################
+
+
+callback_limiter(env, P::Project) = nothing
+
+
+function callback_limiter(env, P::Project{:spherical1d})
+ TODO()
+ @unpack mesh = env
+ @unpack D, Sr, Ï„ = get_dynamic_variables(mesh)
+ for var in (D, Sr, Ï„)
+ limit_slope!(var, P, mesh)
+ end
+end
+
+
+function limit_slope!(var::AbstractArray, mask::AbstractArray, P::Project, mesh::Mesh1d)
+ dg1d.enter(P.prms.workspace) do
+
+ @unpack x = get_static_variables(mesh)
+ Npts = mesh.element.Npts
+ K = mesh.tree.dims[1]
+ tmp_var = dg1d.borrow(P.prms.workspace, Npts*K)
+ avgs = dg1d.borrow(P.prms.workspace, Npts+1)
+ δs = dg1d.borrow(P.prms.workspace, Npts+1)
+ lim_avgs = dg1d.borrow(P.prms.workspace, Npts-1)
+ lim_δs = dg1d.borrow(P.prms.workspace, Npts-1)
+ mat_A = dg1d.borrow(P.prms.workspace, (Npts-1)*(Npts-1))
+ b = dg1d.borrow(P.prms.workspace, Npts-1)
+ A = reshape(mat_A,(Npts-1,Npts-1))
+
+ for k in 1: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)
+ continue
+ end
+
+ idxs_lhs, idxs, idxs_rhs = cellindices.(Ref(mesh), (k-1,k,k+1))
+ vvar_lhs, vvar, vvar_rhs = view.(Ref(var), (idxs_lhs,idxs,idxs_rhs))
+ vx_lhs, vx, vx_rhs = view.(Ref(x), (idxs_lhs,idxs,idxs_rhs))
+ vtmp_var = view(tmp_var, idxs)
+
+ # determine averages and slopes, including neighbors
+ v1, v2 = vvar_lhs[end-1], vvar_lhs[end]
+ x1, x2 = vx_lhs[end-1], vx_lhs[end]
+ avgs[1] = (v1+v2)/2; δs[1] = (v2-v1)/(x2-x1)
+ 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)
+ for i in 1:Npts-1
+ v1, v2 = vvar[i], vvar[i+1]
+ x1, x2 = vx[i], vx[i+1]
+ avgs[1+i] = (v1+v2)/2; δs[1+i] = (v2-v1)/(x2-x1)
+ end
+
+ # limit slopes
+ for i in 2:Npts
+ # 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 = 0.1
+ δl, δc, δr = δs[i-1], δs[i], δs[i+1]
+ lim_δs[i-1] = TCI.minmod_M(dx, M, δc, δl, δr)
+ end
+
+ # recompute avgs so that continuity and total avg is preserved
+ fill!(A, 0.0) # because we use lu! below
+ for i in 1:Npts-2
+ x1, x2, x3 = vx[i], vx[i+1], vx[i+2]
+ dx1, dx2 = x2-x1, x3-x2
+ A[i,i] = 1.0
+ A[i,i+1] = -1.0
+ b[i] = - 0.5*dx1*lim_δs[i] - 0.5*dx2*lim_δs[i+1]
+ end
+ A[end,:] .= 1.0
+ v_avgs = view(avgs, 2:Npts)
+ b[end] = sum(v_avgs)
+
+ # solve the equations
+ F = lu!(A)
+ ldiv!(lim_avgs, F, b)
+
+ # resample data
+ x1, x2 = vx[1], vx[2]
+ dx = x2-x1
+ vtmp_var[1] = lim_avgs[1] - 0.5*dx*lim_δs[1]
+ x1, x2 = vx[end-1], vx[end]
+ dx = x2-x1
+ vtmp_var[end] = lim_avgs[end] + 0.5*dx*lim_δs[end]
+ for i in 2:Npts-1
+ x1, x2 = vx[i], vx[i+1]
+ dx = x2-x1
+ vtmp_var[i] = lim_avgs[i] - 0.5*dx*lim_δs[i]
+ end
+ end
+
+ copyto!(var, tmp_var)
+ end
+end
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index 07a708d6..6a447b29 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -163,15 +163,17 @@ end
γuuxx = 1.0/γxx
γuuzz = 1.0/γzz
- if !c2p_set_atmosphere_on_failure && (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < Ïmin)
- @warn "Require atmosphere but denied by mask" D Ïmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
- TODO("Upsi")
- elseif (c2p_set_atmosphere_on_failure && D < Ïmin) ||
- (c2p_freeze_atm == 0.0 && D < Ïmin) ||
- (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
- if D < Ïmin
- c2p_init_admissible = 0.0
- end
+ # if !c2p_set_atmosphere_on_failure && (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && D < Ïmin)
+ # @warn "Require atmosphere but denied by mask" D Ïmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
+ # @goto impose_atmosphere
+ # TODO("Upsi")
+ # elseif (c2p_set_atmosphere_on_failure && D < Ïmin) ||
+ # (c2p_freeze_atm == 0.0 && D < Ïmin) ||
+ # (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
+ # if D < Ïmin
+ # c2p_init_admissible = 0.0
+ # end
+ if D < Ïmin
@label impose_atmosphere
Ï = Ïatm
vx = 0.0
@@ -191,6 +193,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
@@ -353,14 +356,18 @@ end
W = 1.0 / sqrt(1.0 - v^2)
Ï = D / W
end
- if !c2p_set_atmosphere_on_failure && (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && Ï < Ïmin)
- @warn "Require atmosphere but denied by mask" D Ïmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
- TODO("Upsi")
- elseif (c2p_set_atmosphere_on_failure && Ï < Ïmin) ||
- (c2p_freeze_atm == 0.0 && Ï < Ïmin) ||
- (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
+ if Ï < Ïmin
@goto impose_atmosphere
end
+ # if !c2p_set_atmosphere_on_failure && (c2p_freeze_atm > 0.0 && c2p_reset_atm == 0.0 && Ï < Ïmin)
+ # @warn "Require atmosphere but denied by mask" D Ïmin xcoord zcoord c2p_freeze_atm c2p_reset_atm
+ # @goto impose_atmosphere
+ # TODO("Upsi")
+ # elseif (c2p_set_atmosphere_on_failure && Ï < Ïmin) ||
+ # (c2p_freeze_atm == 0.0 && Ï < Ïmin) ||
+ # (c2p_freeze_atm > 0.0 && c2p_reset_atm > 0.0)
+ # @goto impose_atmosphere
+ # end
# compute thermodynamic vars
if eos isa ColdEquationOfState || y ≤ y0
ϵ = cold_eos(InternalEnergy, Density, Ï)
@@ -395,3 +402,53 @@ 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
diff --git a/src/GRHD/rhs.jl b/src/GRHD/rhs.jl
index 09187952..46e65f58 100644
--- a/src/GRHD/rhs.jl
+++ b/src/GRHD/rhs.jl
@@ -34,8 +34,10 @@ 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 = L/K
+ dl = min_grid_spacing(mesh)
+ # dt = dL / (vmax * dtfactor(mesh))
+ dt = dl / (vmax * dtfactor(mesh))
return dt
end
dtfactor(mesh::Mesh1d{FVElement}) = 2
@@ -1221,6 +1223,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
@@ -1309,6 +1320,197 @@ 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_Ï„ = 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 D_modal, Sr_modal, Ï„_modal,
+ # rhs_D_modal, rhs_Sr_modal, rhs_Ï„_modal = get_static_variables(cache)
+ @unpack flr_D_modal, flr_Sr_modal, flr_Ï„_modal,
+ src_D_modal, src_Sr_modal, src_Ï„_modal = get_static_variables(cache)
+
+ # 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
+
+ # L = layout(mesh)
+ # vD = dg1d.vreshape(D, L)
+ # vSr = dg1d.vreshape(Sr, L)
+ # display(vD)
+ # display(vSr)
+
+ # 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) )
+ flag[k] = 1.0
+ end
+ end
+ K = mesh.tree.dims[1]
+ # flag[K÷2+1] = 1.0
+ 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
+
+ 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
+
+
+ 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
+ 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
+
+ # vrhs_D = dg1d.vreshape(rhs_D, L)
+ # vrhs_Sr = dg1d.vreshape(rhs_Sr, L)
+ # display(vD)
+ # display(vSr)
+ # display(vrhs_D)
+ # display(vrhs_Sr)
+ # println('='^100)
+ # readline()
+
+ 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..f2ef7527 100644
--- a/src/GRHD/setup.jl
+++ b/src/GRHD/setup.jl
@@ -44,11 +44,26 @@ function Project(env::Environment, prms)
fv_numerical_flux = prms["GRHD"]["fv_numerical_flux"]
id = prms["GRHD"]["id"]
bc = prms["GRHD"]["bc"]
+ 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,
id, bc, hrsc, fv_numerical_flux)
# construct Project
@@ -169,7 +184,8 @@ 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_ϵ,
@@ -184,6 +200,13 @@ 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,
+ :rhs_D_modal, :rhs_Sr_modal, :rhs_Ï„_modal,
+ :c2p_failed)
+ )
end
function register_evolution!(mesh::Mesh1d{FVElement}, P::Project{:spherical1d})
hrsc = P.prmsdb["GRHD"]["hrsc"]
@@ -481,7 +504,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/bspline.jl b/src/bspline.jl
index 3919fd54..c53a6a26 100644
--- a/src/bspline.jl
+++ b/src/bspline.jl
@@ -283,6 +283,36 @@ function mass_matrix(bs::Bspline2)
end
+# function mass_matrix(bs1::Bspline1, bs2::Bspline2)
+# Np1, xs1 = bs1.Npts, bs1.xs
+# Np2, xs2 = bs2.Npts, bs2.xs
+# M = zeros(Float64, Np1, Np2+1)
+# zs, ws = dg1d.LG.rule(3) # 4th order accurate
+# o1, o2 = 1, 2 # offset bdry knots
+# xs, o, Np = if Np1 ≤ Np2
+# xs2, o2, Np2
+# else
+# xs1, o1, Np1
+# end
+# for k in 1:Np
+# il, ir = o+k, o+k+1
+# xl, xr = xs[il], xs[ir]
+# for i in 1:Np1, j in 1:Np2+1
+# mij = dg1d.LG.integrate(zs, ws) do z
+# x = ((xr-xl)*z + xr+xl)/2
+# dxdz = (xr-xl)/2
+# B1i = polynomial(x, i, bs1)
+# B2j = polynomial(x, j, bs2)
+# B1i*B2j*dxdz
+# end
+# M[i,j] += mij
+# end
+# end
+# return M
+# end
+# mass_matrix(bs2::Bspline2, bs1::Bspline1) = transpose(mass_matrix(bs1, bs2))
+
+
"""
stiffness_matrix(bs::Bspline1)
--
GitLab