From 5d18c5364b56e34e374bd7a8b6b47a6125658076 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Sun, 19 Nov 2023 20:46:33 +0100
Subject: [PATCH] wip -- perhaps slopy merged
---
src/SRHD/callbacks.jl | 8 +++++--
src/SRHD/cons2prim.jl | 48 ++++++++++++++++++++++-------------------
src/SRHD/equation.jl | 10 +++++++--
src/SRHD/initialdata.jl | 7 ++++++
src/SRHD/rhs.jl | 8 +++----
src/mesh.jl | 5 +++++
6 files changed, 56 insertions(+), 30 deletions(-)
diff --git a/src/SRHD/callbacks.jl b/src/SRHD/callbacks.jl
index 30d72cca..34fcd335 100644
--- a/src/SRHD/callbacks.jl
+++ b/src/SRHD/callbacks.jl
@@ -93,9 +93,13 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
max_v[k] = maximum(view(mat_max_v, :, k))
end
+ # HRSC.compute_viscosity!(
+ # mu, tmp_mu,
+ # (D,S,tau), cellmax_v,
+ # hrsc)
HRSC.compute_viscosity!(
- mu, tmp_mu,
- (D,S,tau), cellmax_v,
+ mu,
+ D, cellmax_v,
hrsc)
end
diff --git a/src/SRHD/cons2prim.jl b/src/SRHD/cons2prim.jl
index 7902f6b3..4fba2d38 100644
--- a/src/SRHD/cons2prim.jl
+++ b/src/SRHD/cons2prim.jl
@@ -57,15 +57,15 @@
# 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
+ # 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)
@@ -100,6 +100,7 @@ end
q = tau / D
S2 = (S/B)^2
Si = S
+ # @assert sign(Si) > 0
ri = Si / D
# r = sqrt(S2) / D # =^= r in paper, r already denotes radius here
z0i = ri / h0
@@ -194,6 +195,7 @@ end
# compute thermodynamic vars
phat = eos(Pressure, Density, rhohat, InternalEnergy, epshat)
+ # phat = eos(Pressure, Density, rhohat)
# hhat = eos(Enthalpy, Density, rhohat, InternalEnergy, epshat)
ahat = phat / (rhohat * (1.0 + epshat))
@@ -230,19 +232,20 @@ end
# eq (15) [Kastaun] : rho*W = D
rhohW2 = D / mu
- # regularize velocity according to Dumbser et al arxiv:2307.06629
- # eqs (82)-(84)
- y, y0 = rhohW2, 1e-4
- if y > y0
- # v = S / mu
- # compute v accoriding to eq (68) or (40) [Kastaun]
- v = sign(Si) * min(mu * r, v0)
- else
- # @warn "limiting v"
- # limit velocity ratio v/S
- fy = 2.0 * (y/y0)^3 - 3.0 * (y/y0)^2 + 1.0
- v = Si * y / (y^2 + fy * 5e-9)
- end
+ v = sign(Si) * min(mu * r, v0)
+ # # regularize velocity according to Dumbser et al arxiv:2307.06629
+ # # eqs (82)-(84)
+ # y, y0 = rhohW2, 1e-4
+ # if y > y0
+ # # v = S / mu
+ # # compute v accoriding to eq (68) or (40) [Kastaun]
+ # v = sign(Si) * min(mu * r, v0)
+ # else
+ # # @warn "limiting v"
+ # # limit velocity ratio v/S
+ # fy = 2.0 * (y/y0)^3 - 3.0 * (y/y0)^2 + 1.0
+ # v = Si * y / (y^2 + fy * 5e-9)
+ # end
# eq (68) or (40)
# v = min(mu * r, v0)
@@ -273,6 +276,7 @@ end
# compute thermodynamic vars
p = eos(Pressure, Density, rho, InternalEnergy, eps)
+ # p = eos(Pressure, Density, rho)
if p <= 0
error((rho,eps))
end
diff --git a/src/SRHD/equation.jl b/src/SRHD/equation.jl
index 655a2989..90638c22 100644
--- a/src/SRHD/equation.jl
+++ b/src/SRHD/equation.jl
@@ -86,14 +86,20 @@ end
nflx = nx*flx_D
bdry_nflx = nx*bdry_flx_D
nflx_D = LLF(nflx, bdry_nflx, D, bdry_D, vmax)
+ # nflx_D = 0.5 * (nflx + bdry_nflx)
+ # display((flx_D, bdry_flx_D, nx, nflx, bdry_nflx, nflx_D))
nflx = nx*flx_S
bdry_nflx = nx*bdry_flx_S
nflx_S = LLF(nflx, bdry_nflx, S, bdry_S, vmax)
+ # nflx_S = 0.5 * (nflx + bdry_nflx)
+ # display((flx_S, bdry_flx_S, nx, nflx, bdry_nflx, nflx_S))
nflx = nx*flx_tau
bdry_nflx = nx*bdry_flx_tau
nflx_tau = LLF(nflx, bdry_nflx, tau, bdry_tau, vmax)
+ # nflx_tau = 0.5 * (nflx + bdry_nflx)
+ # display((flx_D, bdry_flx_D, nx, nflx, bdry_nflx, nflx_D))
@returns [bdry] nflx_D, nflx_S, nflx_tau
end
@@ -180,8 +186,8 @@ end
@with_signature [legacy=false] function av_flux(equations::Equation)
@accepts flx_D, flx_S, flx_tau, ldg_D, ldg_S, ldg_tau, smoothed_mu
- flx_D += smoothed_mu * ldg_D
- flx_S += smoothed_mu * ldg_S
+ flx_D += smoothed_mu * ldg_D
+ flx_S += smoothed_mu * ldg_S
flx_tau += smoothed_mu * ldg_tau
@returns flx_D, flx_S, flx_tau
end
diff --git a/src/SRHD/initialdata.jl b/src/SRHD/initialdata.jl
index 7c98b4ce..64b4d547 100644
--- a/src/SRHD/initialdata.jl
+++ b/src/SRHD/initialdata.jl
@@ -168,6 +168,13 @@ function initialdata_equation_of_state(
xl, xr = -0.5, 0.5
+ Ï0 = [ (xi < xl || xi > xr) ? 1.0 : 0.125 for xi in x ]
+ v0 = [ (xi < xl || xi > xr) ? 0.0 : 0.0 for xi in x ]
+ p0 = [ (xi < xl || xi > xr) ? 1.0 : 0.1 for xi in x ]
+
+ (Ï0, v0, p0)
+end
+
function initialdata_equation_of_state(
::Val{:simple_wave}, eos::EquationOfState.Polytrope, x)
diff --git a/src/SRHD/rhs.jl b/src/SRHD/rhs.jl
index c646a364..2fa10378 100644
--- a/src/SRHD/rhs.jl
+++ b/src/SRHD/rhs.jl
@@ -245,8 +245,8 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
# HRSC.reconstruct!(S, flag, bernstein, isperiodic=true, limit_with_boundaries=false)
# HRSC.reconstruct!(tau, flag, bernstein, isperiodic=true, m=1e-10, limit_with_boundaries=false)
- dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
- dg1d.new_broadcast_volume!(new_speed, equation, mesh)
+ # dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+ # dg1d.new_broadcast_volume!(new_speed, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
dg1d.interpolate_face_data!(mesh, S, bdry_S)
@@ -258,8 +258,8 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
dg1d.interpolate_face_data!(mesh, p, bdry_p)
# # # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
- # dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
- # dg1d.new_broadcast_volume!(new_speed, equation, mesh)
+ dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+ dg1d.new_broadcast_volume!(new_speed, equation, mesh)
dg1d.new_broadcast_volume!(flux, equation, mesh)
dg1d.new_broadcast_volume!(av_flux, equation, mesh)
diff --git a/src/mesh.jl b/src/mesh.jl
index 8bce5113..f88f1a72 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -12,6 +12,11 @@ struct Mesh{N_Dim,N_Sides} <: AbstractMesh
extends::NTuple{N_Dim,Tuple{Float64,Float64}} # total extend of mesh
element::SpectralElement
cache::Cache
+ # TODO Shall we replace these with start:stride:stride. Will require change in @with_signature
+ # to loop over ranges. start should include also cell offsets, because then these
+ # vectors don't have to be as long as the number of cells we have (which we require otherwise
+ # to find the right offset from mesh.offsets)
+ # TODO This should be benchmarked!
bulkfaceindices::Vector{Int64}
bulkbdryindices::Vector{Int64}
faceindices::Vector{Int64}
--
GitLab