diff --git a/examples/simple_wave/aventropy.toml b/examples/simple_wave/aventropy.toml
new file mode 100644
index 0000000000000000000000000000000000000000..7beb0d854df7c42aecfb3927fcd0dc745d4cb324
--- /dev/null
+++ b/examples/simple_wave/aventropy.toml
@@ -0,0 +1,38 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 125
+basis = "lgl"
+periodic = true
+
+[Output]
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v", "EP", "smoothed_mu", "E", "flx_E" ]
+# variables = [ "D", "S", "tau" ]
+aligned_ts = "$(collect(range(0.01,6.4,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 0.125
+tend = 6.4
+# tend = 0.2
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_cmax = 1.0
+entropy_ce = 1.0
+
+[TCI]
+method = "entropy_production"
diff --git a/examples/simple_wave/avmda.toml b/examples/simple_wave/avmda.toml
new file mode 100644
index 0000000000000000000000000000000000000000..e8692e91fad81f3981a8b9103fae7d367a98d825
--- /dev/null
+++ b/examples/simple_wave/avmda.toml
@@ -0,0 +1,35 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 200
+basis = "lgl"
+periodic = true
+
+[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
+
+[Evolution]
+cfl = 0.125
+tend = 6.4
+# tend = 0.2
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "mda"
+
+[TCI]
+method = "entropy_production"
diff --git a/examples/simple_wave/dg.toml b/examples/simple_wave/dg.toml
new file mode 100644
index 0000000000000000000000000000000000000000..028cd5323b3f758eba2216da9177234b237c09f9
--- /dev/null
+++ b/examples/simple_wave/dg.toml
@@ -0,0 +1,38 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+n = 5
+k = 200
+basis = "lgl"
+periodic = false
+
+[Output]
+# variables = [ "D", "S", "tau", "p", "rho", "eps", "E", "Em1", "flx_E", "flx_Em1", "v", "EP", "smoothed_mu" ]
+variables = [ "D", "S", "tau" ]
+aligned_ts = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 0.2
+tend = 0.6
+# tend = 0.2
+
+[Log]
+progress_stdout = true
+
+# [HRSC]
+# method = "av"
+# av_method = "entropy"
+# entropy_cmax = 0.1
+# entropy_ce = 1.0
+
+[TCI]
+method = "entropy_production"
diff --git a/examples/simple_wave/fv.toml b/examples/simple_wave/fv.toml
new file mode 100644
index 0000000000000000000000000000000000000000..97629dd52ad5c346f652a410bea02e829526f229
--- /dev/null
+++ b/examples/simple_wave/fv.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+eos = "polytrope"
+polytrope_gamma = 1.6666666666666666667
+polytrope_k = 100.0
+
+[SRHD]
+id = "simple_wave"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.5, 1.5 ]
+k = 2056
+periodic = false
+scheme = "FV"
+
+[Output]
+variables = [ "D", "S", "tau" ]
+aligned_ts = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 1.0
+tend = 0.6
diff --git a/examples/srhd_blastwave1/avmda.toml b/examples/srhd_blastwave1/avmda.toml
new file mode 100644
index 0000000000000000000000000000000000000000..4615b4468fb4d10b769cdefc7c545184eb7d8c95
--- /dev/null
+++ b/examples/srhd_blastwave1/avmda.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave1"
+bc = "from_id"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+# n = 8
+# k = 255
+n = 5
+k = 205
+basis = "lgl"
+periodic = false
+
+[Output]
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu" ]
+aligned_ts = "$(collect(range(0.01,0.4,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 0.4
+tend = 0.4
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 20.0
+
+[TCI]
+method = "entropy_production"
diff --git a/examples/srhd_blastwave1/fv.toml b/examples/srhd_blastwave1/fv.toml
new file mode 100644
index 0000000000000000000000000000000000000000..6d1977868b9cde0429f26275f0a69895f63f48cd
--- /dev/null
+++ b/examples/srhd_blastwave1/fv.toml
@@ -0,0 +1,22 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave1"
+bc = "from_id"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+k = 1024
+scheme = "FV"
+periodic = false
+
+[Output]
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v" ]
+aligned_ts = "$(collect(range(0.01,0.4,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 1.0
+tend = 0.4
diff --git a/examples/srhd_blastwave2/avmda.toml b/examples/srhd_blastwave2/avmda.toml
new file mode 100644
index 0000000000000000000000000000000000000000..34e785f0327494bedb7d99aa7d3de5d15a7edffa
--- /dev/null
+++ b/examples/srhd_blastwave2/avmda.toml
@@ -0,0 +1,31 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave2"
+bc = "from_id"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+n = 5
+k = 205
+basis = "lgl"
+periodic = false
+
+[Output]
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu" ]
+aligned_ts = "$(collect(range(0.01,0.4,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 0.2
+tend = 0.4
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 20.0
+
+[TCI]
+method = "entropy_production"
diff --git a/examples/srhd_blastwave2/fv.toml b/examples/srhd_blastwave2/fv.toml
new file mode 100644
index 0000000000000000000000000000000000000000..f48beab387050a7b95725f6e852f5118dc86fadc
--- /dev/null
+++ b/examples/srhd_blastwave2/fv.toml
@@ -0,0 +1,22 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.6666666666666666666667
+
+[SRHD]
+id = "blastwave2"
+bc = "from_id"
+
+[Mesh]
+range = [ -0.5, 0.5 ]
+k = 1024
+scheme = "FV"
+periodic = false
+
+[Output]
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v" ]
+aligned_ts = "$(collect(range(0.01,0.4,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 1.0
+tend = 0.4
diff --git a/examples/srhd_sod_shocktube/avmda.toml b/examples/srhd_sod_shocktube/avmda.toml
new file mode 100644
index 0000000000000000000000000000000000000000..298da8b2fe209eaf30663b78aa7cbff0d3c6494f
--- /dev/null
+++ b/examples/srhd_sod_shocktube/avmda.toml
@@ -0,0 +1,35 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.4
+
+[SRHD]
+id = "sod_shocktube"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.0, 1.0 ]
+n = 4
+k = 205
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# variables = [ "D", "S", "tau", "p", "rho", "eps", "v", "smoothed_mu",
+# variables = [ "D", "S", "tau", "smoothed_mu",
+# "ldg_D", "ldg_S", "ldg_tau", "rhs_D", "rhs_S", "rhs_tau" ]
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v" ]
+aligned_ts = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 0.4
+tend = 0.6
+
+[HRSC]
+method = "av"
+av_method = "mda"
+mda_cmax = 20.0
+
+[TCI]
+method = "entropy_production"
diff --git a/examples/srhd_sod_shocktube/fv.toml b/examples/srhd_sod_shocktube/fv.toml
new file mode 100644
index 0000000000000000000000000000000000000000..eecf4daee1cae595ccc1e7a3a632a83619743d5a
--- /dev/null
+++ b/examples/srhd_sod_shocktube/fv.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 1.4
+
+[SRHD]
+id = "sod_shocktube"
+bc = "from_id"
+
+[Mesh]
+range = [ -1.0, 1.0 ]
+k = 1024
+scheme = "FV"
+periodic = false
+
+[Output]
+# every_iteration = 1
+variables = [ "D", "S", "tau", "p", "rho", "eps", "v" ]
+aligned_ts = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d = true
+
+[Evolution]
+cfl = 1.0
+tend = 0.6
diff --git a/src/HRSC/mda.jl b/src/HRSC/mda.jl
index 2827b4215349d3d04ee019e26aeac55c38ecc1df..0337216e73a44ee3f310707df5aa2e95f1c4fded 100644
--- a/src/HRSC/mda.jl
+++ b/src/HRSC/mda.jl
@@ -41,7 +41,7 @@ function compute_viscosity!(
# h = L / (K * N)
h = L / K
- u = reshape(u, layout(mesh))
+ u = dg1d.vreshape(u, layout(mesh))
uhat = invV * u # modal coefficients
# perfect decay modal coefficients, added to enhance sensitivity to
@@ -65,19 +65,33 @@ function compute_viscosity!(
rng_N = 1:N
data_x = @. - log(rng_N)
+ mod_uh = zeros(Float64, N)
+ skylined_mod_uh = zeros(Float64, N)
+ data_y = zeros(Float64, N)
+
for j in 1:K
# (3.11)
# only consider coefficients 2:Np (which corresponds to 1:N), because
# the 1st (0th) one corresponds to a constant mode that is irrelevant
# for the smoothness esimation.
- mod_uh = @. sqrt( uhat[2:end,j]^2 + cellnorm2_u[j] * vhat^2 )
+ @. mod_uh = sqrt( uhat[2:end,j]^2 + cellnorm2_u[j] * vhat^2 )
# (3.12)
# max(1,N-1) inserted to avoid a 0:end range which can happen when only using two points
- skylined_mod_uh = [ maximum(mod_uh[min(n,max(1,N-1)):end]) for n in 1:N ]
+ # skylined_mod_uh = [ maximum(mod_uh[min(n,max(1,N-1)):end]) for n in 1:N ]
+ for n = 1:N
+ skylined_mod_uh[n] = @views maximum(mod_uh[min(n,max(1,N-1)):end])
+ end
# data points in log scale
- data_y = @. log(skylined_mod_uh)
- data_y[isinf.(data_y)] .= -20.0
+ # data_y = @. log(skylined_mod_uh)
+ @. data_y = log(skylined_mod_uh)
+ for n = 1:N
+ y = data_y[n]
+ if isinf(y)
+ data_y[n] = -20.0
+ end
+ end
+ # data_y[isinf.(data_y)] .= -20.0
_, _, tau = dg1d.least_sqr_linear(data_x, data_y)
isnan(tau) && error("undesired nan encountered")
diff --git a/src/SRHD/callbacks.jl b/src/SRHD/callbacks.jl
index 1c27293b619f5f9a9fc809fc727e6a0fc5b25496..0599bc2e8b0361e1d82a3a974d3a7a1ce405c70e 100644
--- a/src/SRHD/callbacks.jl
+++ b/src/SRHD/callbacks.jl
@@ -49,6 +49,7 @@ function callback_equation(state_t, isperiodic, eq, mesh)
Em1 .= E
flx_Em1 .= flx_E
dg1d.new_broadcast_volume!(entropy_variables, eq, mesh)
+ dg1d.new_broadcast_volume!(maxspeed, eq, mesh)
dt <= 0.0 && return
LL = layout(mesh)
@@ -58,7 +59,6 @@ function callback_equation(state_t, isperiodic, eq, mesh)
Em1 = dg1d.vreshape(Em1, LL)
flx_Em1 = dg1d.vreshape(flx_Em1, LL)
invdetJ = dg1d.vreshape(invdetJ, LL)
- v = dg1d.vreshape(v, LL)
# compute entropy residual
@turbo @. flx_Em1 += flx_E
@@ -78,6 +78,7 @@ function callback_equation(state_t, isperiodic, eq::Equation, mesh::Mesh1d{FVEle
Em1 .= E
flx_Em1 .= flx_E
dg1d.new_broadcast_volume!(entropy_variables, eq, mesh)
+ dg1d.new_broadcast_volume!(maxspeed, eq, mesh)
end
@@ -103,12 +104,13 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack t, tm1 = get_global_variables(cache)
- max_v = cellmax_v
- t[1] <= 0.0 && tm1[1] <= 0 && return
- Npts, K = layout(mesh)
- mat_max_v = reshape(view(v, :), (Npts, K))
- for k = 1:K
- max_v[k] = maximum(view(mat_max_v, :, k))
+ # prevent AV recomputation when callbacks are forced on initial data
+ # need this for shock tube tests so that we have non-zero viscosity in troublesome places
+ # although we start with smoothed initial data
+ t[1] <= tm1[1] && return
+
+ for (k,vi) in enumerate(eachcell(mesh, max_v))
+ cellmax_v[k] = maximum(vi)
end
HRSC.compute_viscosity!(
diff --git a/src/SRHD/cons2prim.jl b/src/SRHD/cons2prim.jl
index 7902f6b376836f14ed2b84a5d32ce57db2c61ee1..d8f99cc9a57b4270a57f1fcb6ab4e4a12134ac8c 100644
--- a/src/SRHD/cons2prim.jl
+++ b/src/SRHD/cons2prim.jl
@@ -11,7 +11,7 @@
atm_threshold = equation.atmosphere.threshold
rhomin = atm_density * atm_threshold
- B = 1.0 # in special relativity
+ # B = 1.0 # in special relativity
#####
# conservative variable fixing following Deppe et al. arXiv:2109.12033
@@ -19,17 +19,23 @@
fixed = false
- if D < rhomin
- # @warn "fixing D $D @ x = $x"
+ if D < rhomin || tau < 0.0
D = rhomin
- fixed = true
- end
-
- if tau < 0.0
- # @warn "fixing tau $tau @ x = $x"
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"
@@ -40,33 +46,37 @@
# fixed = true
# end
- # (48)-(50) in arXiv:2109.12033 where we put B=0
- # that = tau / D
- # muhat = 0
- # Solution to (52) when B=0
- # What = 1.0 + tau / D
+ # (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
- # factor = sqrt( (1.0-1e-12) * (What^2-1.0) * D^2 / (S2 + D^2 * 1e-16) )
- # if factor <= 1.0
- # @warn "fixing S $S @ x = $x"
- # S = sign(S) * sqrt( (1.0-1e-12) * (What^2-1.0) * D^2 - D^2 * 5e-12 )
- # S2 = S^2 / B^2
- # factor = sqrt( (1.0-1e-12) * (What^2-1.0) * D^2 / (S2) )
- # @assert factor >= 1.0
- # fixed = true
+ # 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)
@returns D, S, tau, flag_consfixed
@@ -78,7 +88,7 @@ end
#######################################################################
-@with_signature [legacy=false] function cons2prim_kastaun(equation::AbstractEquation)
+@with_signature [simd=false,legacy=false] function cons2prim_kastaun(equation::AbstractEquation)
@accepts D, S, tau, x
@unpack eos = equation
@@ -88,7 +98,7 @@ end
rhomax = Inf
rhoatm = atm_density
# epsmin adjusted for IdealGas
- epsmin, epsmax = 1e-10, Inf
+ epsmin, epsmax = 0.0, Inf
h0 = 1e-3
B = 1.0 # in special relativity
@@ -109,7 +119,7 @@ end
r = abs(ri)
v0 = abs(v0i)
- if v0 >= 1
+ if v0 > 1
error((v0,z0i,Si,D))
end
@@ -134,10 +144,31 @@ 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
- eps = eos.cold_eos(InternalEnergy, Density, rho)
v = 0.0
- p = eos.cold_eos(Pressure, Density, rho)
+ 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
@@ -233,16 +264,16 @@ end
# regularize velocity according to Dumbser et al arxiv:2307.06629
# eqs (82)-(84)
y, y0 = rhohW2, 1e-4
- if y > y0
+ # 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
+ # 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)
@@ -261,23 +292,38 @@ end
# @assert v^2 < 1e-4
# end
- if rho < rhomin || eps < epsmin
- # @warn "atmosphering II @ x = $x ($rho, $eps)"
- # atmosphering
- rho, eps, v, W = max(rhomin, D), epsmin, 0.0, 1.0
- # rho, eps, v, W = rhomin, epsmin, 0.0, 1.0
- elseif rho > rhomax || eps > epsmax
- @warn "rhomax || epsmax reached"
- rho, eps, v = rhomax, epsmax, Si*mu/D
+ @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
+ 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
+ # if rho < rhomin || eps < epsmin
+ # @warn "atmosphering II @ x = $x : (rho,eps) = ($rho, $eps)"
+ # flag_cons2prim_atm2 = 1
+ # # atmosphering
+ # rho, eps, v, W = max(rhomin, D), epsmin, 0.0, 1.0
+ # D = rho
+ # tau =
+ # # rho, eps, v, W = rhomin, epsmin, 0.0, 1.0
+ # elseif rho > rhomax || eps > epsmax
+ # @warn "rhomax || epsmax reached"
+ # flag_cons2prim_atm2 = 1
+ # rho, eps, v = rhomax, epsmax, Si*mu/D
+ # D = rho
+ # end
# compute thermodynamic vars
p = eos(Pressure, Density, rho, InternalEnergy, eps)
- if p <= 0
+ if p < 0
error((rho,eps))
end
end
- @returns rho, v, eps, p
+ # @returns D, S, tau, rho, v, eps, p
+ @returns tau, rho, v, eps, p
end
diff --git a/src/SRHD/equation.jl b/src/SRHD/equation.jl
index fa68c416d11c867374cd5572aabc1d7c48e66629..8521d6cf59d8ec2ea6346029adc0f240508db543 100644
--- a/src/SRHD/equation.jl
+++ b/src/SRHD/equation.jl
@@ -189,6 +189,15 @@ end
end
+@with_signature [legacy=false] function av_flux_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)
+ flx_tau += -smoothed_mu * (rhs_tau - ldg_tau)
+ @returns flx_D, flx_S, flx_tau
+end
+
+
# @with_signature [legacy=false] function av_nflux(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
@@ -277,6 +286,45 @@ end
end
+@with_signature [legacy=false] function av_nflux_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
+
+ nflx = -nx*(rhs_D - ldg_D)
+ bdry_nflx = -nx*(bdry_rhs_D - bdry_ldg_D)
+ nflx_D += 0.5 * (smoothed_mu*nflx + bdry_smoothed_mu*bdry_nflx)
+
+ nflx = -nx*(rhs_S - ldg_S)
+ bdry_nflx = -nx*(bdry_rhs_S - bdry_ldg_S)
+ nflx_S += 0.5 * (smoothed_mu*nflx + bdry_smoothed_mu*bdry_nflx)
+
+ nflx = -nx*(rhs_tau - ldg_tau)
+ bdry_nflx = -nx*(bdry_rhs_tau - bdry_ldg_tau)
+ nflx_tau += 0.5 * (smoothed_mu*nflx + bdry_smoothed_mu*bdry_nflx)
+
+ @returns [bdry] nflx_D, nflx_S, nflx_tau
+end
+@with_signature function bdry_av_nflux_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
+
+ TODO()
+
+ lhs_nflx = nx*lhs_ldg_D
+ nflx_D += lhs_smoothed_mu*lhs_nflx
+
+ lhs_nflx = nx*lhs_ldg_S
+ nflx_S += lhs_smoothed_mu*lhs_nflx
+
+ lhs_nflx = nx*lhs_ldg_tau
+ nflx_tau += lhs_smoothed_mu*lhs_nflx
+
+ @returns nflx_D, nflx_S, nflx_tau
+end
+
+
@with_signature [legacy=false] function fv_bdry_flux(equation::Equation)
@accepts D, S, tau, v, p
@accepts init_D, init_S, init_tau, init_v, init_p
diff --git a/src/SRHD/initialdata.jl b/src/SRHD/initialdata.jl
index 377269f426fa7ef2ab71441ee333623a85f50210..c1bdd32ad3930fea2e7c54d38193f52c67d2253d 100644
--- a/src/SRHD/initialdata.jl
+++ b/src/SRHD/initialdata.jl
@@ -6,6 +6,9 @@ function initialdata!(env, P::Project, prms, mesh::Mesh1d)
initialdata_hrsc!(env, P, P.hrsc, prms, mesh)
initialdata_tci!(env, P, P.tci, prms, mesh)
initialdata_bdrycond!(env, P, dg1d.DirichletBC2(), mesh)
+ if get(P.prms, :id_smooth, false)
+ initialdata_smooth!(env, P, mesh)
+ end
return
end
@@ -163,19 +166,16 @@ end
function initialdata_equation_of_state(
::Val{:sod_shocktube}, eos::EquationOfState.IdealGas, mesh)
- TODO()
-
+ # from EFL paper: arxiv 2202.08839
@unpack x = get_static_variables(mesh.cache)
(xmin, xmax), = mesh.extends
@toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
- @toggled_assert eos.gamma == 5/3
+ @toggled_assert eos.gamma ≈ 1.4
x0 = 0.0
- # Ï0 = [ (xi < x0) ? 1.0 : 0.125 for xi in x ]
- # v0 = [ (xi < x0) ? 0.0 : 0.0 for xi in x ]
- # p0 = [ (xi < x0) ? 1.0 : 0.1 for xi in x ]
-
- xl, xr = -0.5, 0.5
+ Ï0 = [ (xi < x0) ? 1.0 : 0.125 for xi in x ]
+ v0 = [ (xi < x0) ? 0.0 : 0.0 for xi in x ]
+ p0 = [ (xi < x0) ? 1.0 : 0.1 for xi in x ]
(Ï0, v0, p0)
end
@@ -225,7 +225,7 @@ function initialdata_equation_of_state(
@unpack x = get_static_variables(mesh.cache)
(xmin, xmax), = mesh.extends
- @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert isapprox(xmin, -0.5) && isapprox(xmax, 0.5)
@toggled_assert eos.gamma == 5/3
x0 = 0.0
Ï0 = [ (xi < x0) ? 10.0 : 1.0 for xi in x ]
@@ -242,7 +242,7 @@ function initialdata_equation_of_state(
@unpack x = get_static_variables(mesh.cache)
(xmin, xmax), = mesh.extends
- @toggled_assert isapprox(xmin, -1.0) && isapprox(xmax, 1.0)
+ @toggled_assert isapprox(xmin, -0.5) && isapprox(xmax, 0.5)
@toggled_assert eos.gamma == 5/3
x0 = 0.0
Ï0 = [ (xi < x0) ? 1.0 : 1.0 for xi in x ]
@@ -430,7 +430,7 @@ function initialdata_hrsc!(env, P::Project, hrsc::Union{HRSC.EntropyViscosity,HR
@unpack mu, cellmax_v = get_cell_variables(env.cache)
@unpack max_v = get_static_variables(env.cache)
- @unpack D = get_dynamic_variables(env.cache)
+ @unpack D, S, tau = get_dynamic_variables(env.cache)
dg1d.new_broadcast_volume!(maxspeed, P.equation, env.mesh)
Npts, K = layout(env.mesh)
mat_max_v = reshape(view(max_v, :), (Npts, K))
@@ -438,8 +438,8 @@ function initialdata_hrsc!(env, P::Project, hrsc::Union{HRSC.EntropyViscosity,HR
cellmax_v[k] = maximum(view(mat_max_v, :, k))
end
HRSC.compute_viscosity!(
- mu,
- D, cellmax_v,
+ mu, deepcopy(mu),
+ (D,S,tau), cellmax_v,
mda_hrsc.av)
end
@@ -516,3 +516,26 @@ function initialdata_bdrycond!(env, P::Project2d, bdrycond::dg1d.DirichletBC2, m
init_flx_tau_y .= flx_tau_y
init_max_v .= max_v
end
+
+
+#######################################################################
+# Initial data smoothing #
+#######################################################################
+
+
+initialdata_smooth!(env, P::Project, mesh::Mesh{FVElement}) = nothing
+function initialdata_smooth!(env, P::Project, mesh)
+ !isregistered(mesh.cache, :mu) && return
+ @unpack D, S, tau = get_dynamic_variables(mesh.cache)
+ @unpack mu = get_cell_variables(mesh.cache)
+ bernstein = BernsteinReconstruction(mesh)
+ flag = zeros(mesh.K)
+ for k = 1:mesh.K
+ flag[k] = Float64(mu[k] > 0.0)
+ end
+ HRSC.reconstruct!(D, flag, bernstein, isperiodic=true)
+ HRSC.reconstruct!(S, flag, bernstein, isperiodic=true)
+ HRSC.reconstruct!(tau, flag, bernstein, isperiodic=true)
+ dg1d.new_broadcast_volume!(cons2prim_kastaun, P.equation, mesh)
+ dg1d.new_broadcast_volume!(maxspeed, P.equation, mesh)
+end
diff --git a/src/SRHD/rhs.jl b/src/SRHD/rhs.jl
index 2ff0e524fd27bd03f0bac22d961fe652fdd84bc0..f8bccd664309dbb59182506aef293b68857a06d5 100644
--- a/src/SRHD/rhs.jl
+++ b/src/SRHD/rhs.jl
@@ -6,9 +6,9 @@ timestep(env, P, hrsc) = timestep(env, P, hrsc, env.mesh)
function compute_maxspeed(mesh, equation, cache)
- @unpack v = get_static_variables(cache)
+ @unpack max_v = get_static_variables(cache)
dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
- vmax = dg1d.absolute_maximum(v)
+ vmax = dg1d.absolute_maximum(max_v)
vmax_limit = 1e4
if vmax > vmax_limit
@warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
@@ -31,15 +31,15 @@ dtfactor(mesh::Mesh1d{SpectralElement}) = mesh.element.N
function timestep(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, mesh::Mesh1d)
@unpack cache, mesh = env
@unpack equation = P
+ @unpack max_v = get_static_variables(cache)
# dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
- dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
- dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
- max_v = get_variable(cache, :max_v)
+ # dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+ # dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
vmax = dg1d.absolute_maximum(view(max_v,:))
- vmax_limit = 0.99
+ vmax_limit = 0.9999999
if vmax > vmax_limit
- @warn "Limiting timestep due to maximum speed exceeding $vmax_limit"
+ @warn "Limiting timestep due to maximum speed $vmax exceeding $vmax_limit"
vmax = vmax_limit
end
smoothed_mu = get_variable(cache, :smoothed_mu)
@@ -206,11 +206,53 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractReconstruction,
end
+function fd_diff_1d!(_df, _f, _x, mesh)
+ @unpack Npts = mesh.element
+ for (df, f, x) in zip(eachcell(mesh,_df),eachcell(mesh,_f),eachcell(mesh,_x))
+ for n in 2:Npts-1
+ fl,fm,fr = f[n-1],f[n],f[n+1]
+ xl,xm,xr = x[n-1],x[n],x[n+1]
+ dll = (xm-xr)/((xl-xm)*(xl-xr))
+ dlm = (2*xm-xl-xr)/((xm-xl)*(xm-xr))
+ dlr = (xm-xl)/((xr-xl)*(xr-xm))
+ df[n] = fl*dll + fm*dlm + fr*dlr
+ end
+ fl,fm,fr = f[1],f[2],f[3]
+ xl,xm,xr = x[1],x[2],x[3]
+ dll = (2*xl-xr-xm)/((xl-xm)*(xl-xr))
+ dlm = (xl-xr)/((xm-xl)*(xm-xr))
+ dlr = (xl-xm)/((xr-xl)*(xr-xm))
+ df[1] = fl*dll + fm*dlm + fr*dlr
+ fl,fm,fr = f[Npts-2],f[Npts-1],f[Npts]
+ xl,xm,xr = x[Npts-2],x[Npts-1],x[Npts]
+ dll = (xr-xm)/((xl-xm)*(xl-xr))
+ dlm = (xr-xl)/((xm-xl)*(xm-xr))
+ dlr = (2*xr-xl-xm)/((xr-xl)*(xr-xm))
+ df[Npts] = fl*dll + fm*dlm + fr*dlr
+ end
+end
+
+
function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
bdryconds, ldg_bdryconds, av_bdryconds, ::Mesh1d)
- @unpack cache, mesh = env
- @unpack equation = P
+ @unpack cache, mesh = env
+ @unpack equation, prms = P
+ reg = prms.av_regularization
+
+ if reg === :mono
+ rhs_mono!(cache, mesh, equation, P)
+ elseif reg === :covariant
+ rhs_covariant!(cache, mesh, equation, P)
+ else
+ TODO(reg)
+ end
+
+end
+
+
+function rhs_mono!(cache, mesh, equation, P)
+
@unpack D, S, tau = get_dynamic_variables(cache)
@unpack flx_D, flx_S, flx_tau,
ldg_D, ldg_S, ldg_tau,
@@ -224,6 +266,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
bdry_ldg_D, bdry_ldg_S,
bdry_ldg_tau,
bdry_smoothed_mu = get_bdry_variables(cache)
+ @unpack x = get_static_variables(cache)
## solve auxiliary equation: q + ∂x u = 0
dg1d.new_broadcast_faces!(ldg_nflux, equation, mesh)
@@ -232,7 +275,12 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
compute_rhs_weak_form!(ldg_D, D, nflx_D, mesh)
compute_rhs_weak_form!(ldg_S, S, nflx_S, mesh)
compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
+ # fd_diff_1d!(ldg_D, D, x, mesh)
+ # fd_diff_1d!(ldg_S, S, x, mesh)
+ # fd_diff_1d!(ldg_tau, tau, x, mesh)
+
+ # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
@@ -270,6 +318,85 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity,
end
+function rhs_covariant!(cache, mesh, equation, P)
+
+ @unpack D, S, tau = get_dynamic_variables(cache)
+ @unpack flx_D, flx_S, flx_tau,
+ ldg_D, ldg_S, ldg_tau,
+ smoothed_mu, max_v,
+ rho, v, eps, p = get_static_variables(cache)
+ @unpack rhs_D, rhs_S, rhs_tau = get_rhs_variables(cache)
+ @unpack nflx_D, nflx_S, nflx_tau,
+ bdry_D, bdry_S, bdry_tau,
+ bdry_rho, bdry_v, bdry_eps,
+ bdry_p, bdry_max_v,
+ bdry_ldg_D, bdry_ldg_S, bdry_ldg_tau,
+ bdry_rhs_D, bdry_rhs_S, bdry_rhs_tau,
+ bdry_smoothed_mu = get_bdry_variables(cache)
+ @unpack x = get_static_variables(cache)
+
+ # compute rhs for ∂t u + ∂x f(x) = 0
+
+ # dg1d.new_broadcast_volume!(conservative_fixing, equation, mesh)
+ dg1d.new_broadcast_volume!(cons2prim_kastaun, equation, mesh)
+ dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
+
+ dg1d.interpolate_face_data!(mesh, D, bdry_D)
+ dg1d.interpolate_face_data!(mesh, S, bdry_S)
+ dg1d.interpolate_face_data!(mesh, tau, bdry_tau)
+ dg1d.interpolate_face_data!(mesh, max_v, bdry_max_v)
+ dg1d.interpolate_face_data!(mesh, rho, bdry_rho)
+ dg1d.interpolate_face_data!(mesh, v, bdry_v)
+ dg1d.interpolate_face_data!(mesh, eps, bdry_eps)
+ dg1d.interpolate_face_data!(mesh, p, bdry_p)
+
+ dg1d.new_broadcast_volume!(flux, equation, mesh)
+
+ # numerical fluxes
+ dg1d.new_broadcast_faces!(llf, equation, mesh)
+ dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
+
+ compute_rhs_weak_form!(rhs_D, flx_D, nflx_D, mesh)
+ compute_rhs_weak_form!(rhs_S, flx_S, nflx_S, mesh)
+ compute_rhs_weak_form!(rhs_tau, flx_tau, nflx_tau, mesh)
+
+ # now compute rhs for ∂t u + ∂x f(x) = ∂x( μ (-∂t+∂x) u )
+
+ ## solve auxiliary equation: q + ∂x u = 0
+ dg1d.new_broadcast_faces!(ldg_nflux, equation, mesh)
+ dg1d.new_broadcast_bdry!(ldg_bdryllf, equation, mesh)
+ compute_rhs_weak_form!(ldg_D, D, nflx_D, mesh)
+ compute_rhs_weak_form!(ldg_S, S, nflx_S, mesh)
+ compute_rhs_weak_form!(ldg_tau, tau, nflx_tau, mesh)
+ # fd_diff_1d!(ldg_D, D, x, mesh)
+ # fd_diff_1d!(ldg_S, S, x, mesh)
+ # fd_diff_1d!(ldg_tau, tau, x, mesh)
+
+ # need to recompute, because ldg computation overwrites nflx_{D,S,tau}
+ # numerical fluxes
+ dg1d.new_broadcast_faces!(llf, equation, mesh)
+ dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
+
+ dg1d.new_broadcast_volume!(av_flux_covariant, equation, mesh)
+
+ dg1d.interpolate_face_data!(mesh, ldg_D, bdry_ldg_D)
+ dg1d.interpolate_face_data!(mesh, ldg_S, bdry_ldg_S)
+ dg1d.interpolate_face_data!(mesh, ldg_tau, bdry_ldg_tau)
+ dg1d.interpolate_face_data!(mesh, rhs_D, bdry_rhs_D)
+ dg1d.interpolate_face_data!(mesh, rhs_S, bdry_rhs_S)
+ dg1d.interpolate_face_data!(mesh, rhs_tau, bdry_rhs_tau)
+ dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
+
+ dg1d.new_broadcast_faces!(av_nflux_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)
+ compute_rhs_weak_form!(rhs_tau, flx_tau, nflx_tau, mesh)
+
+ return
+end
+
+
#######################################################################
# 2d #
#######################################################################
diff --git a/src/SRHD/setup.jl b/src/SRHD/setup.jl
index f6658523aa089bcd3e610b1d91d5f6f79683df3e..f7efa9bb0fce473f614c3b9d6235476053bbd4cc 100644
--- a/src/SRHD/setup.jl
+++ b/src/SRHD/setup.jl
@@ -9,7 +9,8 @@ function Project(env::Environment, prms, mesh::Mesh1d)
tci = TCI.make_TCI(env.mesh, prms["TCI"])
hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
- P = Project(equation, hrsc, tci)
+ fixedprms = (; av_regularization=:covariant, id_smooth=true)
+ P = Project(equation, hrsc, tci, fixedprms)
# register variables
# TODO Somehow replace _register_variables! with register_variables!
@@ -55,7 +56,7 @@ function Project(env::Environment, prms, mesh::Mesh2d)
# TODO Somehow replace _register_variables! with register_variables!
@unpack cache = env
_register_variables!(mesh, P)
- _register_variables!(mesh, hrsc)
+ _register_variables!(mesh, hrsc, nothing)
_register_variables!(mesh, bdrycond)
display(cache)
@@ -108,23 +109,23 @@ function _register_variables!(mesh::Mesh, P::Project)
cell_variablenames = (:cellmax_v,),
global_variablenames = (:t, :tm1) # time steps
)
- _register_variables!(mesh, P.hrsc)
+ _register_variables!(mesh, P.hrsc, P.prms)
_register_variables!(mesh, P.tci)
@unpack x = get_static_variables(mesh.cache)
# @. x = env.mesh.x[:]
end
-function _register_variables!(mesh, tci_or_hrsc::Nothing) end
+function _register_variables!(mesh, tci_or_hrsc::Nothing, args...) end
-_register_variables!(mesh, hrsc::AbstractHRSC) =
+_register_variables!(mesh, hrsc::AbstractHRSC, prms) =
error("_register_variables: Not implemented for hrsc of type '$(typeof(hrsc))'")
-_register_variables!(mesh, hrsc::HRSC.AbstractReconstruction) = nothing
+_register_variables!(mesh, hrsc::HRSC.AbstractReconstruction, prms) = nothing
-function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity)
+function _register_variables!(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),
@@ -133,17 +134,23 @@ function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity)
end
-function _register_variables!(mesh, hrsc::HRSC.SmoothedArtificialViscosity)
- _register_variables!(mesh, hrsc.av)
+function _register_variables!(mesh, hrsc::HRSC.SmoothedArtificialViscosity, prms)
+ _register_variables!(mesh, hrsc.av, prms)
register_variables!(mesh.cache,
static_variablenames = (:smoothed_mu,),
bdry_variablenames = (:bdry_ldg_D, :bdry_ldg_S, :bdry_ldg_tau,
:bdry_smoothed_mu),
)
+ reg = prms.av_regularization
+ if reg == :covariant
+ register_variables!(mesh.cache,
+ bdry_variablenames = (:bdry_rhs_D, :bdry_rhs_S, :bdry_rhs_tau)
+ )
+ end
end
-function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity)
+function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity, prms)
register_variables!(mesh.cache,
static_variablenames = (:E, :Em1, :flx_E_x, :flx_E_y, :flxm1_E_x, :flxm1_E_y,
:ldg_D_x, :ldg_D_y, :ldg_Sx_x, :ldg_Sx_y,
diff --git a/src/SRHD/types.jl b/src/SRHD/types.jl
index ea57c51bb3e068b6ab317d0c7bc99dee9e573d66..c406928205294408511cd2ec531692adb39885df 100644
--- a/src/SRHD/types.jl
+++ b/src/SRHD/types.jl
@@ -21,11 +21,13 @@ end
struct Project{T_HRSC <:Maybe{HRSC.AbstractHRSC},
- T_TCI <:Maybe{TCI.AbstractTCI}}
+ T_TCI <:Maybe{TCI.AbstractTCI},
+ T_Prms}
equation::Equation
hrsc::T_HRSC
tci::T_TCI
+ prms::T_Prms
end
diff --git a/test/IntegrationTests/refs/srhd_simple_wave_avmda/output.h5 b/test/IntegrationTests/refs/srhd_simple_wave_avmda/output.h5
index 645018934bcea95dd6bfe1c5bb38653df410ae60..d862f71838d8f971e66e0000249a4884267eda71 100644
Binary files a/test/IntegrationTests/refs/srhd_simple_wave_avmda/output.h5 and b/test/IntegrationTests/refs/srhd_simple_wave_avmda/output.h5 differ
diff --git a/test/IntegrationTests/refs/srhd_simple_wave_fv/output.h5 b/test/IntegrationTests/refs/srhd_simple_wave_fv/output.h5
index 57e5df637099fbfb066c714db8f3d45f4c482e5b..88c8bb9c4de42a67ea15d6e1c8d33285f11658d3 100644
Binary files a/test/IntegrationTests/refs/srhd_simple_wave_fv/output.h5 and b/test/IntegrationTests/refs/srhd_simple_wave_fv/output.h5 differ