diff --git a/examples/simple_wave/avmda.toml b/examples/simple_wave/avmda_covariant.toml
similarity index 57%
rename from examples/simple_wave/avmda.toml
rename to examples/simple_wave/avmda_covariant.toml
index e8692e91fad81f3981a8b9103fae7d367a98d825..6c1b86aebf29326f0773d38b669bc0aa33b06240 100644
--- a/examples/simple_wave/avmda.toml
+++ b/examples/simple_wave/avmda_covariant.toml
@@ -6,23 +6,22 @@ polytrope_k = 100.0
[SRHD]
id = "simple_wave"
bc = "from_id"
+av_regularization = "mono_approx_covariant"
[Mesh]
range = [ -1.5, 1.5 ]
n = 5
-k = 200
+k = 45
basis = "lgl"
-periodic = true
+periodic = false
[Output]
-variables = [ "D", "S", "tau", "p", "rho", "eps", "v", "EP", "smoothed_mu", "E", "flx_E" ]
-aligned_ts = "$(collect(range(0.01,6.4,step=0.01)))"
-enable1d = true
+variables1d = [ "D", "S", "tau" ]
+aligned_ts = "$(collect(range(0.01,1.2,step=0.01)))"
[Evolution]
-cfl = 0.125
-tend = 6.4
-# tend = 0.2
+cfl = 0.2
+tend = 1.2
[Log]
progress_stdout = true
diff --git a/examples/simple_wave/avmda_mono.toml b/examples/simple_wave/avmda_mono.toml
new file mode 100644
index 0000000000000000000000000000000000000000..c99ebe5162fcb1f9609d769ea5abb594c7153abe
--- /dev/null
+++ b/examples/simple_wave/avmda_mono.toml
@@ -0,0 +1,34 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+av_regularization = "mono"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 45
+basis = "lgl"
+periodic = false
+
+[Output]
+variables1d = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu" ]
+aligned_ts = "$(collect(range(0.01,1.2,step=0.01)))"
+
+[Evolution]
+cfl = 0.2
+tend = 1.2
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "mda"
+
+[TCI]
+method = "entropy_production"
diff --git a/src/GRHD/cons2prim.jl b/src/GRHD/cons2prim.jl
index 74fb83d6e8fa244c01d3996eb4095f3e15f7d2d4..81747833e9601bc7b801c951bf9ec23b3114f565 100644
--- a/src/GRHD/cons2prim.jl
+++ b/src/GRHD/cons2prim.jl
@@ -133,9 +133,10 @@ end
parameters,
D, Sx, Sz, Ï„,
xcoord, zcoord, γxx, γzz,
- c2p_reset_atm, #=TODO remove this arg=#c2p_limit_vr, c2p_freeze_atm)
+ c2p_reset_atm, c2p_limit_vr, c2p_freeze_atm)
verbose = false
+ # TODO Remove these
c2p_reset_ϵ, c2p_init_admissible = 0.0, 1.0
# setup some constants
@unpack atm_density, atm_threshold, c2p_threshold_velocity_limiter, eos, cold_eos = parameters
@@ -356,9 +357,9 @@ end
p = cold_eos(Pressure, Density, Ï)
ϵ = cold_eos(InternalEnergy, Density, Ï)
end
- if verbose && abs(ϵ) > 10.0
- @warn "abs(ϵ) > 1.0" Sx Sz Ï„ Ï W D ϵ y v vx vz xcoord zcoord c2p_limit_vr c2p_freeze_atm c2p_reset_atm c2p_init_admissible v2 (q-mu*r2) mu q r2
- end
+ # if verbose && abs(ϵ) > 10.0
+ # @warn "abs(ϵ) > 1.0" Sx Sz Ï„ Ï W D ϵ y v vx vz xcoord zcoord c2p_limit_vr c2p_freeze_atm c2p_reset_atm c2p_init_admissible v2 (q-mu*r2) mu q r2
+ # end
# check from PHYSICAL REVIEW D 87, 084006 (2013), appendix A.2
# seems to be always be satisfied here
diff --git a/src/SRHD/SRHD.jl b/src/SRHD/SRHD.jl
index 8b506a6625d5ca32e5ba413f54c198dd16b974e6..dc23c9339ab255c8851a2c96f17c1c523a8c34f6 100644
--- a/src/SRHD/SRHD.jl
+++ b/src/SRHD/SRHD.jl
@@ -43,10 +43,11 @@ dg1d.@parameters SRHD begin
atm_threshold_factor = 100.0
@check atm_threshold_factor > 0.0
+ # TODO I think the @check macro is not working properly...
"""
Initial data configuration. Available options:
- - `id_smooth`
- - `id_smooth_diagonal`
+ - `id_smooth_flow`
+ - `id_smooth_flow_diagonal`
- `id_periodic_relativistic_shocktube`
- `id_periodic_relativistic_shocktube_2`
- `id_periodic_shocktube_pressure_jump`
@@ -54,23 +55,31 @@ dg1d.@parameters SRHD begin
- `id_periodic_sod_shocktube`
- `id_sod_shocktube`
- `id_shocktube`
+ - `id_simple_wave`
"""
- id = "id_smooth"
- # TODO I think the @check macro is not working properly...
- @check id in [ "id_smooth",
- "id_smooth_diagonal",
+ id = "id_smooth_flow"
+ @check id in [ "id_smooth_flow",
+ "id_smooth_flow_diagonal",
"id_periodic_relativistic_shocktube",
"id_periodic_relativistic_shocktube_2",
"id_periodic_shocktube_pressure_jump",
"id_periodic_shocktube_less_extreme",
"id_periodic_sod_shocktube",
"id_sod_shocktube",
- "id_shocktube_2d" ]
+ "id_shocktube_2d",
+ "id_simple_wave"
+ ]
"Flow speed for id_smooth test"
id_smooth_flow_speed = 0.2
@check id_smooth_flow_speed > 0.0
+ """
+ If `true` smooth initial data using a Bernstein filter.
+ """
+ id_smooth = false
+ @check id_smooth isa Bool
+
"""
Boundary conditions. Available options:
- `periodic`
@@ -108,6 +117,37 @@ dg1d.@parameters SRHD begin
slope_limiter_tvb_M = 150.0
@check slope_limiter_tvb_M > 0.0
+ """
+ Type of artificial viscosity regularization. Available options:
+ - `none`
+ - `mono` - monolithic regularization, spatial derivatives only
+ - `mono_approx_covariant` - monolithic regularization with approximate covariant terms of the form ∂x(μ(∂t - ∂x)u)
+ """
+ av_regularization = "none"
+ @check av_regularization in ["none", "mono", "mono_covariant"]
+
+ """
+ The parameters `[K, Γ]` of a polytropic EoS that is used as a fall back to model cold matter near the fluid-atmosphere interface.
+ """
+ c2p_cold_eos_parameters = [ 100.0, 2.0 ]
+ @check length(c2p_cold_eos_parameters) == 2
+ @check c2p_cold_eos_parameters[1] > 0.0
+ @check c2p_cold_eos_parameters[2] > 0.0
+
+ """
+ Determines how quickly artificial viscosity can alter between time steps.
+ If `av_drag = 0.0` then no dragging occurs.
+ If `av_drag = 1.0` then viscosity is frozen to its initial value.
+ """
+ av_drag = 0.0
+ @check 0.0 <= av_drag <= 1.0
+
+ """
+ If `true` any artificial viscosity related callbacks are run on every substep of the evolution.
+ """
+ av_recompute_substeps = false
+ @check av_recompute_substeps isa Bool
+
end
diff --git a/src/SRHD/callbacks.jl b/src/SRHD/callbacks.jl
index b849aaaf782a6bbd8650f21b585bd48464459cc1..8987a60a9f41b29024e633780f1f15570710107e 100644
--- a/src/SRHD/callbacks.jl
+++ b/src/SRHD/callbacks.jl
@@ -8,8 +8,11 @@ function make_callback(env, P::Project, isperiodic)
cbfn_hrsc(u, t) = callback_hrsc(u, t, env, P, isperiodic, P.hrsc, env.mesh)
cb_hrsc = FunctionCallback(cbfn_hrsc,
CallbackTiming(every_iteration=1,every_dt=0))
- callbackset = CallbackSet(cb_equation, cb_tci, cb_hrsc)
- return callbackset
+ if P.prms.av_recompute_substeps && P.hrsc isa AbstractArtificialViscosity
+ return CallbackSet(cb_equation, cb_tci), CallbackSet(cb_hrsc)
+ else
+ return CallbackSet(cb_equation, cb_tci, cb_hrsc), nothing
+ end
end
@@ -129,10 +132,14 @@ end
function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity, mesh)
- callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc.av, mesh)
+ wrap_dynamic_variables!(env.cache, state_u)
@unpack cache = env
- @unpack mu = get_cell_variables(cache)
+ @unpack mu, mu_m1 = get_cell_variables(cache)
@unpack smoothed_mu = get_static_variables(cache)
+ @unpack av_drag = P.prms
+ mu_m1 .= mu
+ callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc.av, mesh)
+ @. mu = ((1.0-av_drag)*mu+av_drag*mu_m1)
hrsc.smoother(smoothed_mu, mu, mesh, false)
end
function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.SmoothedArtificialViscosity, mesh::Mesh2d)
diff --git a/src/SRHD/cons2prim.jl b/src/SRHD/cons2prim.jl
index 48188350d7640d4283766cde3931150946bfa21d..ed623aafd10e75b0d13cdeabf07371ae8c64247c 100644
--- a/src/SRHD/cons2prim.jl
+++ b/src/SRHD/cons2prim.jl
@@ -1,301 +1,12 @@
-#######################################################################
-# Conservative fixing #
-#######################################################################
-
-
-@with_signature function conservative_fixing(equation::AbstractEquation)
- @accepts D, S, tau, x
-
- @unpack eos = equation
- atm_density = equation.atmosphere.rho
- atm_threshold = equation.atmosphere.threshold
- rhomin = atm_density * atm_threshold
-
- # B = 1.0 # in special relativity
-
- #####
- # conservative variable fixing following Deppe et al. arXiv:2109.12033
- ######
-
- fixed = false
-
- if D < rhomin || tau < 0.0
- D = rhomin
- tau = 1e-12
- fixed = true
- end
-
- # if D < 0.001 * rhomin
- # # @warn "fixing D $D @ x = $x"
- # D = rhomin
- # fixed = true
- # end
- #
- # if tau < 0.0
- # # @warn "fixing tau $tau @ x = $x"
- # tau = 1e-12
- # fixed = true
- # end
- # fix tau such that v < 0.999
- # if tau < 0.0
- # @warn "fixing tau $tau @ x = $x"
- # # What = 1 + tau / D = 1 / sqrt(1 - v^2)
- # v2 = 0.999
- # tau = max(D * (1.0 / sqrt(1.0 - v2) - 1), 1e-12)
- # display(tau)
- # fixed = true
- # end
-
- # (48)-(50) in arXiv:2109.12033 we put B=0, ̂μ = 0
- # solution to (52) is then
- What = 1.0 + tau / 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 @ x = $x"
- 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, tau, flag_consfixed
-end
-
-
-#######################################################################
-# cons2prim a la Kastaun et al 2021 #
-#######################################################################
-
-
@with_signature [simd=false] function cons2prim_kastaun(equation::AbstractEquation)
@accepts D, S, tau, x
-
- @unpack eos = equation
- atm_density = equation.atmosphere.rho
- atm_threshold = equation.atmosphere.threshold
- rhomin = atm_density * atm_threshold
- rhomax = Inf
- rhoatm = atm_density
- # epsmin adjusted for IdealGas
- epsmin, epsmax = 0.0, Inf
- h0 = 1e-3
-
- ####
- # cons2prim from Kastaun et al. arXiv:2005.01821
- ####
-
- q = tau / D
- S2 = S^2
- Si = S
- ri = Si / D
- # r = sqrt(S2) / D # =^= r in paper, r already denotes radius here
- z0i = ri / h0
- v0i = z0i / sqrt(1.0 + z0i^2)
- h02 = h0^2
- r2 = ri^2
- r = abs(ri)
- v0 = abs(v0i)
-
- if v0 > 1
- error((v0,z0i,Si,D))
- end
-
- # What is not a question but rather W^hat, cf. (40)
- fn_What = let r = r, v0 = v0
- mu -> begin
- vhat = min(mu * r, v0)
- return 1.0 / sqrt(1.0 - vhat^2)
- end
- end
-
- # setup bracket
- mu_a = 0.0
- # with no EM fields we can find the root of (49) analytically, because r = const.
- mu_plus = 1.0 / sqrt(h02 + r2)
- mu_b = r < h0 ? 1.0 / h0 : mu_plus
-
- What_a = fn_What(mu_a)
- What_b = fn_What(mu_b)
- if D / What_b > rhomax
- error("Invalid state encountered!")
- end
-
- if D < rhomin
- @warn "atmosphering I"
-
- error((D,S,tau))
-
- @assert eos isa EquationOfState.IdealGas
- @unpack Gamma = eos
-
- # employ atmosphere values by taking the zero temperature limit
- # for an ideal gas eos, we have that eos.Gamma becomes the Gamma for
- # a polytrope, see grhd notes for derivation
-
- rho = rhoatm
- v = 0.0
- W = 1.0
- # p = K * rho^Gamma
- # eps = Gamma * rho^
- # D = W * rho
- # S = rho * (1+eps+p/rho) * W^2 * v
- # tau = rho * (1+eps+p/rho) * W^2 - p - D
-
- # if eos isa ColdEquationOfState
- # eps = eos.cold_eos(InternalEnergy, Density, rho)
- # p = eos.cold_eos(Pressure, Density, rho)
- # else
- # eps = eos(InternalEnergy, Density, rho)
- # p = eos(Pressure, Density, rho)
- # end
-
- else
-
- # root bracket adjustment, cf. Appendix
- if D / What_b < rhomax < D
- @warn "adjusting upper root bracket"
-
- # infamous closure bug ...
- # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
- fn_rhomax = let D = D, rhomax = rhomax
- mu -> D / fn_What(mu) - rhomax
- end
-
- # mu_b = find_zero(fn_rhomax, (mu_a, mu_plus), Roots.A42())
- mu_b = find_zero(fn_rhomax, (mu_a, mu_plus), Roots.Brent())
-
- if isnan(mu_b)
- error("Failed to correct upper interval bound")
- end
-
- println("New interval: ($mu_a, $mu_b)")
- end
-
- if D / What_b < rhomin < D # yes, we need What_b here
- @warn "adjusting lower root bracket"
-
- # infamous closure bug ...
- # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
- fn_rhomin = let D = D, rhomin = rhomin
- mu -> D / fn_What(mu) - rhomin
- end
-
- # mu_a = find_zero(fn_rhomin, (mu_a, mu_b), Roots.A42())
- mu_a = find_zero(fn_rhomin, (mu_a, mu_b), Roots.Brent())
-
- if isnan(mu_b)
- error("Failed to correct lower interval bound")
- end
-
- println("New interval: ($mu_a, $mu_b)")
- end
-
- # infamous closure bug ...
- # https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured
- master_fn = let D = D, v0 = v0, r = r, q = q, r2 = r2
- mu -> begin # eq (44)
-
- _vhat = min(mu * r, v0)
- _What = 1.0 / sqrt(1.0 - _vhat^2)
- rhohat = D / _What
-
- # eq (42)
- epshat = _What * (q - mu * r2) + _vhat^2 * _What^2 / (1.0 + _What)
-
- # TODO enforce EOS bounds - needed?
-
- # compute thermodynamic vars
- phat = eos(Pressure, Density, rhohat, InternalEnergy, epshat)
- # hhat = eos(Enthalpy, Density, rhohat, InternalEnergy, epshat)
- ahat = phat / (rhohat * (1.0 + epshat))
-
- # eq (46)-(48)
- nu_A = (1.0 + ahat) * (1.0 + epshat) / _What
- nu_B = (1.0 + ahat) * (1.0 + q - mu * r2)
- nu = max(nu_A, nu_B)
-
- # eq (44)
- f = mu - 1.0 / (nu + r2 * mu)
- return f
- end
- end
-
- # verify if there is a solution inside bracket
- mf_a = master_fn(mu_a)
- mf_b = master_fn(mu_b)
- mf_ab = mf_a*mf_b
- if mf_a * mf_b >= 0 || !isfinite(mf_ab)
- println((; D, S, tau, mf_a, mf_b, mu_a, mu_b, x))
- error("Master function is ill conditioned!")
- end
-
- # mu = find_zero(master_fn, (mu_a, mu_b), Roots.A42() #= =^= algorithm 748 =#)
- mu = find_zero(master_fn, (mu_a, mu_b), Roots.Brent())
- if isnan(mu)
- error("Failed to find root")
- end
-
- #####
- # evaluate primitives
- #####
-
- # eq (36) [Kastaun] : h*W = 1/mu
- # eq (15) [Kastaun] : rho*W = D
- rhohW2 = D / mu
-
- v = sign(Si) * min(mu * r, v0)
-
- # eq (68) or (40)
- W = 1.0 / sqrt(1.0 - v^2)
- # eq (41)
- rho = D / W
- # eq (42)
- eps = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W)
- # eq (38)
- # eps2 = W * (q - mu * r2) + W - 1
-
- # enforce EOS bounds
-
- @assert rho > rhomin
- if eps < epsmin
- # adjusting tau such that energy density is within bounds
- # cf. 3rd paragraph in sec. III.A in Kastaun paper
- TODO("really need that for SRHD?")
- eps = epsmin
- taub4 = tau
- # invert this for tau: eps = W * (q - mu * r2) + v^2 * W^2 / (1.0 + W), q = tau / D
- tau = (eps / W - v^2 * W / (1.0 + W) + mu * r2) * D
- end
-
- # compute thermodynamic vars
- p = eos(Pressure, Density, rho, InternalEnergy, eps)
- if p < 0
- error((rho,eps))
- end
-
- end
-
- # @returns D, S, tau, rho, v, eps, p
- @returns tau, rho, v, eps, p
+ prms = (; eos=equation.eos, atm_density=1e-10, atm_threshold=100.0,
+ c2p_threshold_velocity_limiter=-1.0 #=disable limiter=#,
+ cold_eos=equation.cold_eos)
+ D, S, _, tau, rho, v, _, _, eps, p, _, _, _, _, _, _ =
+ GRHD.cons2prim_kastaun_impl(prms,
+ D, #=Sx=#S, #=Sz=#0.0, tau,
+ x, #=z=#0.0, #=γxx=#1.0, #=γzz=#1.0,
+ #=c2p_reset/limit/freeze=# 0, 0, 0)
+ @returns D, S, tau, rho, v, eps, p
end
diff --git a/src/SRHD/equation.jl b/src/SRHD/equation.jl
index ed87ba7661d93a1c4083059eb3aa366b96113844..22d7e2558fc190f1cec5c8c68badea6da3b260f8 100644
--- a/src/SRHD/equation.jl
+++ b/src/SRHD/equation.jl
@@ -1,10 +1,10 @@
-function make_Equation(eos, prms, m::Mesh)
+function make_Equation(eos, cold_eos, prms, m::Mesh)
@unpack atm_density, atm_threshold_factor = prms
atm = Atmosphere(atm_density, atm_threshold_factor)
if m isa Mesh1d
- return Equation(eos, atm)
+ return Equation(eos, cold_eos, atm)
elseif m isa Mesh2d
- return Equation2d(eos, atm)
+ return Equation2d(eos, cold_eos, atm)
else
error("unknown mesh type encountered")
end
@@ -178,7 +178,7 @@ end
end
-@with_signature function av_flux(equations::Equation)
+@with_signature function av_flux_mono(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
@@ -187,7 +187,7 @@ end
end
-@with_signature function av_flux_covariant(equations::Equation)
+@with_signature function av_flux_approx_mono_covariant(equations::Equation)
@accepts rhs_D, rhs_S, rhs_tau, flx_D, flx_S, flx_tau, ldg_D, ldg_S, ldg_tau, smoothed_mu
flx_D += -smoothed_mu * (rhs_D - ldg_D)
flx_S += -smoothed_mu * (rhs_S - ldg_S)
@@ -245,7 +245,7 @@ end
# end
-@with_signature function av_nflux(eq::Equation)
+@with_signature function av_nflux_mono(eq::Equation)
@accepts ldg_D, ldg_S, ldg_tau, smoothed_mu
@accepts [bdry] bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau, bdry_smoothed_mu
@accepts [bdry] nflx_D, nflx_S, nflx_tau, nx
@@ -264,7 +264,7 @@ end
@returns [bdry] nflx_D, nflx_S, nflx_tau
end
-# @with_signature function bdry_av_nflux(eq::Equation)
+# @with_signature function bdry_av_nflux_mono(eq::Equation)
# @accepts Prefix(lhs), ldg_D, ldg_S, ldg_tau, smoothed_mu
# @accepts Prefix(rhs), ldg_D, ldg_S, ldg_tau, smoothed_mu
# @accepts nflx_D, nflx_S, nflx_tau, nx
@@ -284,7 +284,7 @@ end
# end
-@with_signature function av_nflux_covariant(eq::Equation)
+@with_signature function av_nflux_approx_mono_covariant(eq::Equation)
@accepts rhs_D, rhs_S, rhs_tau, ldg_D, ldg_S, ldg_tau, smoothed_mu
@accepts [bdry] bdry_rhs_D, bdry_rhs_S, bdry_rhs_tau, bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau, bdry_smoothed_mu
@accepts [bdry] nflx_D, nflx_S, nflx_tau, nx
@@ -303,7 +303,7 @@ end
@returns [bdry] nflx_D, nflx_S, nflx_tau
end
-# @with_signature function bdry_av_nflux_covariant(eq::Equation)
+# @with_signature function bdry_av_nflux_approx_mono_covariant(eq::Equation)
# @accepts Prefix(lhs), ldg_D, ldg_S, ldg_tau, smoothed_mu
# @accepts Prefix(rhs), ldg_D, ldg_S, ldg_tau, smoothed_mu
# @accepts nflx_D, nflx_S, nflx_tau, nx
@@ -374,40 +374,6 @@ module Formulae2d
end
-function conservative_fixing(D, Sx, Sy, tau, equation)
-
- @unpack atmosphere, eos = equation
- atm_density, atm_threshold = atmosphere.rho, atmosphere.threshold
- rhomin = atm_density * atm_threshold
- rhoatm = atm_density
-
- if D < rhomin
- D = rhoatm
- end
-
- if tau < 0
- tau = 1e-12
- end
-
- # (48)-(50) in Deppe paper where we put B=0
- # that = tau / D
- # muhat = 0
- # Solution to (52) when B=0
- What = 1 + tau / D
-
- # in spherical symmetry: S^2 = S_r^2 / B^2
- # Hence, Smax^2/S^2 = Srmax^2 / Sr^2
- S2 = Sx^2 + Sy^2
- factor = sqrt( (1-1e-12) * (What^2-1) * D^2 / (S2 + D^2 * 1e-16) )
- if factor > 1
- Sx /= factor * (1 + 1e-8)
- Sy /= factor * (1 + 1e-8)
- end
-
- return D, Sx, Sy, tau
-end
-
-
# @with_signature function flux(equation::Equation2d)
# @accepts D, Sx, Sy, rho, eps, vx, vy, p, tau
# h = 1.0 + eps + p / rho
@@ -620,7 +586,7 @@ end
# end
#
#
-# @with_signature function av_flux_2d(eq::Equation2d)
+# @with_signature function av_flux_mono(eq::Equation2d)
# @accepts flx_D_x, flx_D_y, ldg_D_x, ldg_D_y,
# flx_Sx_x, flx_Sx_y, ldg_Sx_x, ldg_Sx_y,
# flx_Sy_x, flx_Sy_y, ldg_Sy_x, ldg_Sy_y,
@@ -642,7 +608,7 @@ end
# @returns flx_D_x, flx_D_y, flx_Sx_x, flx_Sx_y,
# flx_Sy_x, flx_Sy_y, flx_tau_x, flx_tau_y
# end
-# @with_signature function av_nflux_2d(eq::Equation2d)
+# @with_signature function av_nflux_mono(eq::Equation2d)
# @accepts Prefix(lhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu
# @accepts Prefix(rhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu
# @accepts nflx_D, nflx_Sx, nflx_Sy, nflx_tau, nx, ny
@@ -665,7 +631,7 @@ end
#
# @returns nflx_D, nflx_Sx, nflx_Sy, nflx_tau
# end
-# @with_signature function bdry_av_nflux_2d(eq::Equation2d)
+# @with_signature function bdry_av_nflux_mono(eq::Equation2d)
# @accepts Prefix(lhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu
# @accepts Prefix(rhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu
# @accepts nflx_D, nflx_Sx, nflx_Sy, nflx_tau, nx, ny
@@ -686,7 +652,7 @@ end
# end
-# @meshloop function av_nflux_2d(out, in1, in2, in3, mesh, eq::Equation2d)
+# @meshloop function av_nflux_mono(out, in1, in2, in3, mesh, eq::Equation2d)
# @accepts Prefix(lhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu = in1
# @accepts Prefix(rhs), ldg_D_x, ldg_D_y, ldg_Sx_x, ldg_Sx_y, ldg_Sy_x, ldg_Sy_y, ldg_tau_x, ldg_tau_y, smoothed_mu = ini2
# @accepts nflx_D, nflx_Sx, nflx_Sy, nflx_tau, nx, ny = in3
diff --git a/src/SRHD/initialdata.jl b/src/SRHD/initialdata.jl
index f2369635992a4517d51c571d0615f87935d139eb..6d53f4094e98ec12ec696e8b1f58ca54fbb5aa86 100644
--- a/src/SRHD/initialdata.jl
+++ b/src/SRHD/initialdata.jl
@@ -44,19 +44,22 @@ function initialdata_equation_of_state!(env, P::Project, eos, prms, mesh::Mesh1d
@. eps = ϵ0
@. rho = Ï0
- # broadcast_volume!(conservative_fixing, P.equation, mesh)
broadcast_volume!(cons2prim_kastaun, P.equation, mesh)
+ # TODO()
+
return
end
+@nospecialize
function initialdata_equation_of_state(::Val{type}, eos, mesh) where type
error("initialdata_equation_of_state: Unknown initial data type '$type'")
end
+@specialize
-function initialdata_equation_of_state(::Val{:smooth}, eos::EquationOfState.IdealGas, mesh)
+function initialdata_equation_of_state(::Val{:smooth_flow}, eos::EquationOfState.IdealGas, mesh)
@unpack x = get_static_variables(mesh.cache)
(xmin, xmax), = mesh.extends
@@ -191,7 +194,6 @@ function initialdata_equation_of_state(
@unpack x = get_static_variables(mesh.cache)
(xmin, xmax), = mesh.extends
- @toggled_assert isapprox(xmin, -1.5) && isapprox(xmax, 1.5)
@toggled_assert eos.K ≈ 100
@toggled_assert eos.Gamma ≈ 5/3
X, a, Γ, K = 0.3, 0.5, eos.Gamma, eos.K
@@ -270,7 +272,7 @@ function initialdata_equation_of_state!(env, P, eos::EquationOfState.IdealGas, p
Ly = xmax-xmin
@unpack gamma = eos
- Ï0, vx0, vy0, p0 = if id == "smooth_diagonal"
+ Ï0, vx0, vy0, p0 = if id == "smooth_flow_diagonal"
@unpack id_smooth_flow_speed = prms
@toggled_assert gamma == 5/3
@@ -322,7 +324,6 @@ function initialdata_equation_of_state!(env, P, eos::EquationOfState.IdealGas, p
@. Sy = Formulae2d.Sy(Ï0, vx0, vy0, ϵ0, p0)
@. tau = Formulae2d.Ï„(Ï0, vx0, vy0, ϵ0, p0)
- # broadcast_volume!(conservative_fixing, P.equation, mesh)
broadcast_volume!(cons2prim_kastaun, P.equation, mesh)
@unpack vx, vy = get_static_variables(cache)
@@ -411,6 +412,7 @@ function initialdata_hrsc!(env, P::Project2d, hrsc::HRSC.AbstractArtificialVisco
end
function initialdata_hrsc!(env, P::Project, hrsc::Union{HRSC.EntropyViscosity,HRSC.MDAViscosity}, prms, ::Mesh1d)
+ return
@unpack mu = get_cell_variables(env.cache)
@unpack E, Em1, flx_E, flx_Em1 = get_static_variables(env.cache)
@unpack t, tm1 = get_global_variables(env.cache)
diff --git a/src/SRHD/rhs.jl b/src/SRHD/rhs.jl
index c24e61aed693b9bdcc5ea9dbdd5165ac2fbddf18..058634c15ebd760c8fd0e1e9901eefffab51064a 100644
--- a/src/SRHD/rhs.jl
+++ b/src/SRHD/rhs.jl
@@ -33,9 +33,6 @@ function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh:
@unpack equation = P
@unpack max_v = get_static_variables(cache)
- # broadcast_volume!(conservative_fixing, equation, mesh)
- # broadcast_volume!(cons2prim_kastaun, equation, mesh)
- # broadcast_volume!(maxspeed, equation, mesh)
vmax = maximum(abs, max_v)
vmax_limit = 0.9999999
if vmax > vmax_limit
@@ -367,7 +364,6 @@ function rhs!(env, P::Project, hrsc::Nothing, mesh::Mesh1d)
bdry_max_v, bdry_v,
bdry_rho, bdry_eps, bdry_p = get_bdry_variables(cache)
- # broadcast_volume!(conservative_fixing, equation, mesh)
broadcast_volume!(cons2prim_kastaun, equation, mesh)
broadcast_volume!(maxspeed, equation, mesh)
@@ -409,7 +405,6 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction, ::Mesh1d)
HRSC.reconstruct!(S, flag, hrsc, isperiodic=isperiodic(mesh), limit_with_boundaries=false)
HRSC.reconstruct!(tau, flag, hrsc, isperiodic=isperiodic(mesh), m=1e-10, limit_with_boundaries=false)
- # broadcast_volume!(conservative_fixing, equation, mesh)
broadcast_volume!(cons2prim_kastaun, equation, mesh)
dg1d.interpolate_face_data!(mesh, D, bdry_D)
dg1d.interpolate_face_data!(mesh, S, bdry_S)
@@ -465,10 +460,10 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d)
if reg === :mono
rhs_mono!(cache, mesh, equation, P)
- elseif reg === :covariant
- rhs_covariant!(cache, mesh, equation, P)
+ elseif reg === :mono_approx_covariant
+ rhs_mono_approx_covariant!(cache, mesh, equation, P)
else
- TODO(reg)
+ TODO("artificial viscosity with av_regularization = $reg")
end
end
@@ -503,7 +498,6 @@ function rhs_mono!(cache, mesh, equation, P)
# fd_diff_1d!(ldg_tau, tau, x, mesh)
- # broadcast_volume!(conservative_fixing, equation, mesh)
broadcast_volume!(cons2prim_kastaun, equation, mesh)
broadcast_volume!(maxspeed, equation, mesh)
@@ -516,12 +510,8 @@ function rhs_mono!(cache, mesh, equation, P)
dg1d.interpolate_face_data!(mesh, eps, bdry_eps)
dg1d.interpolate_face_data!(mesh, p, bdry_p)
- # # # broadcast_volume!(conservative_fixing, equation, mesh)
- # broadcast_volume!(cons2prim_kastaun, equation, mesh)
- # broadcast_volume!(maxspeed, equation, mesh)
-
broadcast_volume!(flux, equation, mesh)
- broadcast_volume!(av_flux, equation, mesh)
+ broadcast_volume!(av_flux_mono, equation, mesh)
dg1d.interpolate_face_data!(mesh, ldg_D, bdry_ldg_D)
dg1d.interpolate_face_data!(mesh, ldg_S, bdry_ldg_S)
@@ -530,7 +520,7 @@ function rhs_mono!(cache, mesh, equation, P)
# numerical fluxes
broadcast_faces!(llf, equation, mesh)
- broadcast_faces!(av_nflux, equation, mesh)
+ broadcast_faces!(av_nflux_mono, equation, mesh)
broadcast_bdry!(bdryllf, equation, mesh)
compute_rhs_weak_form!(rhs_D, flx_D, nflx_D, mesh)
@@ -541,7 +531,7 @@ function rhs_mono!(cache, mesh, equation, P)
end
-function rhs_covariant!(cache, mesh, equation, P)
+function rhs_mono_approx_covariant!(cache, mesh, equation, P)
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack flx_D, flx_S, flx_tau,
@@ -560,7 +550,6 @@ function rhs_covariant!(cache, mesh, equation, P)
# compute rhs for ∂t u + ∂x f(x) = 0
- # broadcast_volume!(conservative_fixing, equation, mesh)
broadcast_volume!(cons2prim_kastaun, equation, mesh)
broadcast_volume!(maxspeed, equation, mesh)
@@ -600,7 +589,7 @@ function rhs_covariant!(cache, mesh, equation, P)
broadcast_faces!(llf, equation, mesh)
broadcast_bdry!(bdryllf, equation, mesh)
- broadcast_volume!(av_flux_covariant, equation, mesh)
+ broadcast_volume!(av_flux_approx_mono_covariant, equation, mesh)
dg1d.interpolate_face_data!(mesh, ldg_D, bdry_ldg_D)
dg1d.interpolate_face_data!(mesh, ldg_S, bdry_ldg_S)
@@ -610,7 +599,7 @@ function rhs_covariant!(cache, mesh, equation, P)
dg1d.interpolate_face_data!(mesh, rhs_tau, bdry_rhs_tau)
dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
- broadcast_faces!(av_nflux_covariant, equation, mesh)
+ broadcast_faces!(av_nflux_approx_mono_covariant, equation, mesh)
compute_rhs_weak_form!(rhs_D, flx_D, nflx_D, mesh)
compute_rhs_weak_form!(rhs_S, flx_S, nflx_S, mesh)
@@ -706,9 +695,9 @@ function dg1d.rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, me
broadcast_faces!(llf, equation, mesh)
broadcast_bdry!(bdry_llf, equation, P.bdrycond, mesh)
- broadcast_volume!(av_flux_2d, equation, mesh)
- broadcast_faces!(av_nflux_2d, equation, mesh)
- broadcast_bdry!(bdry_av_nflux_2d, equation, P.bdrycond, mesh)
+ broadcast_volume!(av_flux_mono, equation, mesh)
+ broadcast_faces!(av_nflux_mono, equation, mesh)
+ broadcast_bdry!(bdry_av_nflux_mono, equation, P.bdrycond, mesh)
compute_rhs_weak_form!(rhs_D, flx_D_x, flx_D_y, nflx_D, mesh)
compute_rhs_weak_form!(rhs_Sx, flx_Sx_x, flx_Sx_y, nflx_Sx, mesh)
diff --git a/src/SRHD/setup.jl b/src/SRHD/setup.jl
index c40d9ccd27ec4376025fa3ed32c0e2f1104963d4..172a0c15cf2e99d153f2d4228468f27e1535f28f 100644
--- a/src/SRHD/setup.jl
+++ b/src/SRHD/setup.jl
@@ -7,12 +7,18 @@ function Project(env::Environment, prms, mesh::Mesh1d)
slope_limiter_method = Symbol(prms["SRHD"]["slope_limiter_method"])
slope_limiter_tvb_M = prms["SRHD"]["slope_limiter_tvb_M"]
bernstein = BernsteinReconstruction(mesh)
- fixedprms = (; av_regularization=:covariant, id_smooth=true,
+ av_regularization = Symbol(prms["SRHD"]["av_regularization"])
+ av_drag = prms["SRHD"]["av_drag"]
+ av_recompute_substeps = prms["SRHD"]["av_recompute_substeps"]
+ K, Γ = prms["SRHD"]["c2p_cold_eos_parameters"]
+ fixedprms = (; av_regularization, av_drag, av_recompute_substeps,
+ id_smooth=prms["SRHD"]["id_smooth"],
bernstein, hrsc, slope_limiter_method, slope_limiter_tvb_M)
# construct Project
eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
- equation = make_Equation(eos, prms["SRHD"], mesh)
+ cold_eos = Polytrope(K, Γ)
+ equation = make_Equation(eos, cold_eos, prms["SRHD"], mesh)
tci = TCI.make_TCI(env.mesh, prms["TCI"])
hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
@@ -37,9 +43,10 @@ function Project(env::Environment, prms, mesh::Mesh1d)
initialdata!(env, P, prms["SRHD"])
# setup callbacks
- projectcb = make_callback(env, P, isperiodic(mesh))
+ cbs, cbs_substep = make_callback(env, P, isperiodic(mesh))
- append!(env.callbacks, CallbackSet(projectcb.callbacks...))
+ append!(env.callbacks, cbs)
+ append!(env.callbacks_substep, cbs_substep)
return P, (bdryconds, ldg_bdryconds, av_bdryconds)
end
@@ -47,8 +54,11 @@ end
function Project(env::Environment, prms, mesh::Mesh2d)
+ K, Γ = prms["SRHD"]["c2p_cold_eos_parameters"]
+
# construct Project
eos = EquationOfState.make_EquationOfState(prms["EquationOfState"])
+ cold_eos = Polytrope(K, Γ)
hrsc = HRSC.make_HRSC(mesh, prms["HRSC"])
equation = make_Equation(eos, prms["SRHD"], mesh)
bdrycond = dg1d.DirichletBC2()
@@ -173,7 +183,7 @@ function register!(mesh, hrsc::HRSC.AbstractArtificialViscosity, prms)
register_variables!(mesh.cache,
static_variablenames = (:ldg_D, :ldg_S, :ldg_tau,
:flx_ldg_D, :flx_ldg_S, :flx_ldg_tau),
- cell_variablenames = (:mu,:tmp_mu), # 'one viscosity to rule them all'
+ cell_variablenames = (:mu,:mu_m1,:tmp_mu), # 'one viscosity to rule them all'
)
end
@@ -186,7 +196,7 @@ function register!(mesh, hrsc::HRSC.SmoothedArtificialViscosity, prms)
:bdry_smoothed_mu),
)
reg = prms.av_regularization
- if reg == :covariant
+ if reg == :mono_approx_covariant
register_variables!(mesh.cache,
bdry_variablenames = (:bdry_rhs_D, :bdry_rhs_S, :bdry_rhs_tau)
)
@@ -201,7 +211,7 @@ function register!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity, prms)
:ldg_Sy_x, :ldg_Sy_y, :ldg_tau_x, :ldg_tau_y),
bdry_variablenames = (:nflx_D_x, :nflx_D_y, :nflx_Sx_x, :nflx_Sx_y,
:nflx_Sy_x, :nflx_Sy_y, :nflx_tau_x, :nflx_tau_y),
- cell_variablenames = (:mu,),
+ cell_variablenames = (:mu,:mu_m1,),
global_variablenames = (:tm1,)
)
end
diff --git a/src/SRHD/types.jl b/src/SRHD/types.jl
index 0c7accd2c2faaf1ad4464bb18e1459e4b41838af..a8eeb494c67980b0b82d705242035972bb70e439 100644
--- a/src/SRHD/types.jl
+++ b/src/SRHD/types.jl
@@ -6,12 +6,14 @@ end
struct Equation{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEquation
eos::T_EoS
+ cold_eos::Polytrope
atmosphere::Atmosphere
end
struct Equation2d{T_EoS<:EquationOfState.AbstractEquationOfState} <: AbstractEquation
eos::T_EoS
+ cold_eos::Polytrope
atmosphere::Atmosphere
end
diff --git a/src/callbacks.jl b/src/callbacks.jl
index 7e84b9964e8df2f7e359094aeca96f24dec22aea..f1c99d86b1238d3c36c1d344e4b2edaa5aae3db1 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -207,6 +207,7 @@ Base.push!(cbs::CallbackSet, cb::AbstractCallback) = push!(cbs.callbacks, cb)
Base.push!(cbs::CallbackSet, ::Nothing) = nothing
Base.append!(cbs::CallbackSet, cbs2::AbstractVector{<:AbstractCallback}) = append!(cbs.callbacks, cbs2)
Base.append!(cbs::CallbackSet, cbs2::CallbackSet) = append!(cbs.callbacks, cbs2.callbacks)
+Base.append!(cbs::CallbackSet, ::Nothing) = nothing
function Base.show(io::IO, cbs::CallbackSet)
diff --git a/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5 b/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5
index 1b2a7a7d05d421d85d6d5db1f333b025c796a594..9a4df45cfa893b4ca513f922fb23259118d19a45 100644
Binary files a/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5 and b/test/IntegrationTests/refs/srhd_simple_wave_avmda/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/srhd_simple_wave_avmda/srhd_simple_wave_avmda.toml b/test/IntegrationTests/refs/srhd_simple_wave_avmda/srhd_simple_wave_avmda.toml
index bc0cb0aa69d7ff3205d4e54dabe68724f3ee2a7c..82666c3bd5fff93f41aa1180c1f20f74b509bdd0 100644
--- a/test/IntegrationTests/refs/srhd_simple_wave_avmda/srhd_simple_wave_avmda.toml
+++ b/test/IntegrationTests/refs/srhd_simple_wave_avmda/srhd_simple_wave_avmda.toml
@@ -6,6 +6,7 @@ polytrope_k = 100.0
[SRHD]
id = "simple_wave"
bc = "from_id"
+av_regularization = "mono_approx_covariant"
[Mesh]
range = [ -1.5, 1.5 ]
diff --git a/test/IntegrationTests/refs/srhd_smooth/output1d.h5 b/test/IntegrationTests/refs/srhd_smooth_flow/output1d.h5
similarity index 100%
rename from test/IntegrationTests/refs/srhd_smooth/output1d.h5
rename to test/IntegrationTests/refs/srhd_smooth_flow/output1d.h5
diff --git a/test/IntegrationTests/refs/srhd_smooth/srhd_smooth.toml b/test/IntegrationTests/refs/srhd_smooth_flow/srhd_smooth_flow.toml
similarity index 93%
rename from test/IntegrationTests/refs/srhd_smooth/srhd_smooth.toml
rename to test/IntegrationTests/refs/srhd_smooth_flow/srhd_smooth_flow.toml
index d422f742ec7e2aa1ac52d9097db0542c09bd6e8f..ad4c528f53f08bb178d6849db1d1913ffeb9b756 100644
--- a/test/IntegrationTests/refs/srhd_smooth/srhd_smooth.toml
+++ b/test/IntegrationTests/refs/srhd_smooth_flow/srhd_smooth_flow.toml
@@ -3,7 +3,7 @@ eos = "idealgas"
idealgas_gamma = 1.6666666666666666666667
[SRHD]
-id = "smooth"
+id = "smooth_flow"
[Mesh]
range = [ -1.0, 1.0 ]
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index 99ad41a94c2de81c0c23236dd3c9ceac54d31d82..5f96a686864962059f7f2ce2b01c12ba4879ef5f 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -92,7 +92,7 @@ variables1d = [ "rho", "q", "E" ]
# [euler_sod_shock_tube_entropyav]
# variables = [ "rho", "q", "E" ]
-[srhd_smooth]
+[srhd_smooth_flow]
# stops before discontinuity appears
variables1d = [ "D", "S", "tau" ]