From 2dda124f244ad7a4815a56d1357fbd4978580a65 Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Sun, 3 Sep 2023 20:01:43 +0000
Subject: [PATCH] Make EulerEq2d run with Kelvin Helmholtz test (dg/dg1d.jl!54)
* add TODO on bdryconds not being implemented for AV 2d
* fixup dg_rhs in special case when one flux component is nonzero
* cleanup initialdata
* make EulerEq2d run with Kelvin Helmholtz test (not stable yet)
also substitute all max with absmax in equations.jl
* export broadcast_bdry_2!
* clean up EulerEq/equation.jl; remove local LLF definitions
* move out numerical utils in separate file
---
src/EulerEq/equation.jl | 238 +++++++++++-------
src/EulerEq/initialdata.jl | 24 +-
src/EulerEq/rhs.jl | 15 +-
src/EulerEq/setup.jl | 13 +-
src/SRHD/equation.jl | 3 -
src/ScalarEq/equation.jl | 3 -
src/dg1d.jl | 17 +-
src/dg_rhs.jl | 30 ++-
src/numutils.jl | 106 ++++++++
src/utils.jl | 88 -------
.../euler2d_kelvin_helmholtz_entropyav.toml | 8 +-
11 files changed, 335 insertions(+), 210 deletions(-)
create mode 100644 src/numutils.jl
diff --git a/src/EulerEq/equation.jl b/src/EulerEq/equation.jl
index e68070cd..7e8a3f3c 100644
--- a/src/EulerEq/equation.jl
+++ b/src/EulerEq/equation.jl
@@ -43,7 +43,7 @@ end
@with_signature function entropy(equation::EulerEquation{EquationOfState.IdealGas})
- @accepts rho, q, p, x
+ @accepts rho, p, x
@unpack eos = equation
gamma = eos.gamma
S = rho / (gamma - 1) * log(p / rho^gamma)
@@ -64,7 +64,7 @@ end
@accepts Prefix(rhs), rho, flx_rho, q, flx_q, E, flx_E, v
@accepts nx
- vmax = max(lhs_v, rhs_v)
+ vmax = absmax(lhs_v, rhs_v)
lhs_nflx = nx*lhs_flx_rho
rhs_nflx = nx*rhs_flx_rho
@@ -108,7 +108,7 @@ end
rhs_flx_q = rhs_q^2 / rhs_rho + rhs_p
rhs_flx_E = rhs_q / rhs_rho * (rhs_E + rhs_p)
- vmax = max(lhs_v, rhs_init_v)
+ vmax = absmax(lhs_v, rhs_init_v)
lhs_nflx = nx*lhs_flx_rho
rhs_nflx = nx*rhs_flx_rho
@@ -156,17 +156,41 @@ end
end
+@with_signature function ldg_nflux(eq::EulerEquation)
+ @accepts Prefix(lhs), rho, q, E
+ @accepts Prefix(rhs), rho, q, E
+ @accepts nx
+
+ nflx_rho = 0.5 * nx * (lhs_rho + rhs_rho)
+ nflx_q = 0.5 * nx * (lhs_q + rhs_q)
+ nflx_E = 0.5 * nx * (lhs_E + rhs_E)
+
+ @returns nflx_rho, nflx_q, nflx_E
+end
+@with_signature function ldg_bdryllf(eq::EulerEquation)
+ @accepts Prefix(lhs), rho, q, E
+ @accepts Prefix(rhs), rho, q, E
+ @accepts nx
+
+ # set exterior state to 0 according to stability analysis
+ # see dg-notes/tex/evm.pdf for reasoning
+ nflx_rho = nx * lhs_rho
+ nflx_q = nx * lhs_q
+ nflx_E = nx * lhs_E
+
+ @returns nflx_rho, nflx_q, nflx_E
+end
+
+
#######################################################################
# 2D #
#######################################################################
-LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
-
-
@with_signature function internal_energy(::EulerEquation2d)
- @accepts rho, q, E
- eps = E / rho - q^2 / (2 * rho^2)
+ @accepts rho, qx, qy, E
+ q2 = qx^2 + qy^2
+ eps = E / rho - q2 / (2 * rho^2)
@returns eps
end
@@ -174,7 +198,8 @@ end
@with_signature function pressure(equation::EulerEquation2d)
@accepts rho, qx, qy, E
q = sqrt(qx^2+qy^2)
- eps, = internal_energy((rho,q,E), equation)
+ eps, = internal_energy((rho,qx,qy,E), equation)
+ @unpack eos = equation
p = eos(Pressure, Density, rho, InternalEnergy, eps)
@returns p
end
@@ -184,7 +209,7 @@ end
@accepts rho, qx, qy, E, p
@unpack eos = equation
q = sqrt(qx^2+qy^2)
- eps, = internal_energy((rho,q,E), equation)
+ eps, = internal_energy((rho,qx,qy,E), equation)
p = eos(Pressure, Density, rho, InternalEnergy, eps)
flx_rho_x = qx
flx_rho_y = qy
@@ -197,13 +222,14 @@ end
end
+# TODO Replace with llf(::EulerEquation2d) below
@with_signature function llf_rho(equation::EulerEquation2d)
@accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y, v
@accepts Prefix(rhs), rho, flx_rho_x, flx_rho_y, v
@accepts nx, ny
lhs_nflx = nx*lhs_flx_rho_x + ny*lhs_flx_rho_y
rhs_nflx = nx*rhs_flx_rho_x + ny*rhs_flx_rho_y
- vmax = max(lhs_v, rhs_v)
+ vmax = absmax(lhs_v, rhs_v)
nflx_rho = LLF(lhs_nflx, rhs_nflx, lhs_rho, rhs_rho, vmax)
@returns nflx_rho
end
@@ -213,7 +239,7 @@ end
@accepts nx, ny
lhs_nflx = nx*lhs_flx_qx_x + ny*lhs_flx_qx_y
rhs_nflx = nx*rhs_flx_qx_x + ny*rhs_flx_qx_y
- vmax = max(lhs_v, rhs_v)
+ vmax = absmax(lhs_v, rhs_v)
nflx_qx = LLF(lhs_nflx, rhs_nflx, lhs_qx, rhs_qx, vmax)
@returns nflx_qx
end
@@ -223,7 +249,7 @@ end
@accepts nx, ny
lhs_nflx = nx*lhs_flx_qy_x + ny*lhs_flx_qy_y
rhs_nflx = nx*rhs_flx_qy_x + ny*rhs_flx_qy_y
- vmax = max(lhs_v, rhs_v)
+ vmax = absmax(lhs_v, rhs_v)
nflx_qy = LLF(lhs_nflx, rhs_nflx, lhs_qy, rhs_qy, vmax)
@returns nflx_qy
end
@@ -233,18 +259,18 @@ end
@accepts nx, ny
lhs_nflx = nx*lhs_flx_E_x + ny*lhs_flx_E_y
rhs_nflx = nx*rhs_flx_E_x + ny*rhs_flx_E_y
- vmax = max(lhs_v, rhs_v)
+ vmax = absmax(lhs_v, rhs_v)
nflx_E = LLF(lhs_nflx, rhs_nflx, lhs_E, rhs_E, vmax)
@returns nflx_E
end
@with_signature function llf(equation::EulerEquation2d)
- @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_x, flx_qy_y, E, flx_E_x, flx_E_y, v
- @accepts Prefix(rhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_x, flx_qy_y, E, flx_E_x, flx_E_y, v
+ @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_y, E, flx_E_x, flx_E_y, v
+ @accepts Prefix(rhs), rho, flx_rho_x, flx_rho_y, qx, flx_qx_x, flx_qx_y, qy, flx_qy_y, E, flx_E_x, flx_E_y, v
@accepts nx, ny
- vmax = max(lhs_v, rhs_v)
+ vmax = absmax(lhs_v, rhs_v)
lhs_nflx = nx*lhs_flx_rho_x + ny*lhs_flx_rho_y
rhs_nflx = nx*rhs_flx_rho_x + ny*rhs_flx_rho_y
@@ -254,6 +280,8 @@ end
rhs_nflx = nx*rhs_flx_qx_x + ny*rhs_flx_qx_y
nflx_qx = LLF(lhs_nflx, rhs_nflx, lhs_qx, rhs_qx, vmax)
+ lhs_flx_qy_x = lhs_flx_qx_y
+ rhs_flx_qy_x = rhs_flx_qx_y
lhs_nflx = nx*lhs_flx_qy_x + ny*lhs_flx_qy_y
rhs_nflx = nx*rhs_flx_qy_x + ny*rhs_flx_qy_y
nflx_qy = LLF(lhs_nflx, rhs_nflx, lhs_qy, rhs_qy, vmax)
@@ -266,11 +294,73 @@ end
end
+@with_signature function bdryllf(equation::EulerEquation2d)
+ @accepts Prefix(lhs), rho, flx_rho_x, flx_rho_y,
+ qx, flx_qx_x, flx_qx_y,
+ qy, #=flx_qy_x,=# flx_qy_y,
+ E, flx_E_x, flx_E_y, v
+ @accepts init_rho, init_qx, init_qy, init_E, init_v
+ @accepts nx, ny
+
+ lhs_flx_qy_x = lhs_flx_qx_y
+
+ # enforce dirichlet/neumann conditions following stability analysis
+ # see dg-notes/tex/evm.pdf
+ # Hesthaven, Numerical Methods for Conservation Laws uses same logic
+
+ # Dirichlet conditions
+ rhs_rho = -lhs_rho + 2*init_rho
+ rhs_qx = -lhs_qx + 2*init_qx
+ rhs_qy = -lhs_qy + 2*init_qy
+ rhs_E = -lhs_E + 2*init_E
+
+ # Neumann conditions
+ # rhs_rho = lhs_rho
+ # rhs_qx = lhs_qx
+ # rhs_qy = lhs_qy
+ # rhs_E = lhs_E
+
+ error()
+ @unpack eos = equation
+ rhs_q = sqrt(rhs_qx^2 + rhs_qy^2)
+ rhs_eps, = internal_energy((rhs_rho,rhs_qx,rhs_qy,rhs_E), equation)
+ rhs_p = eos(Pressure, Density, rhs_rho, InternalEnergy, rhs_eps)
+ rhs_flx_rho_x = rhs_qx
+ rhs_flx_rho_y = rhs_qy
+ rhs_flx_qx_x = rhs_qx^2 / rhs_rho + rhs_p
+ rhs_flx_qx_y = rhs_qx*rhs_qy / rhs_rho
+ rhs_flx_qy_x = rhs_flx_qx_y
+ rhs_flx_qy_y = rhs_qy^2 / rhs_rho + rhs_p
+ rhs_flx_E_x = rhs_qx / rhs_rho * (rhs_E + rhs_p)
+ rhs_flx_E_y = rhs_qy / rhs_rho * (rhs_E + rhs_p)
+
+ vmax = absmax(lhs_v, init_v)
+
+ lhs_nflx = nx*lhs_flx_rho_x + ny*lhs_flx_rho_y
+ rhs_nflx = nx*rhs_flx_rho_x + ny*rhs_flx_rho_y
+ nflx_rho = LLF(lhs_nflx, rhs_nflx, lhs_rho, rhs_rho, vmax)
+
+ lhs_nflx = nx*lhs_flx_qx_x + ny*lhs_flx_qx_y
+ rhs_nflx = nx*rhs_flx_qx_x + ny*rhs_flx_qx_y
+ nflx_qx = LLF(lhs_nflx, rhs_nflx, lhs_qx, rhs_qx, vmax)
+
+ lhs_nflx = nx*lhs_flx_qy_x + ny*lhs_flx_qy_y
+ rhs_nflx = nx*rhs_flx_qy_x + ny*rhs_flx_qy_y
+ nflx_qy = LLF(lhs_nflx, rhs_nflx, lhs_qy, rhs_qy, vmax)
+
+ lhs_nflx = nx*lhs_flx_E_x + ny*lhs_flx_E_y
+ rhs_nflx = nx*rhs_flx_E_x + ny*rhs_flx_E_y
+ nflx_E = LLF(lhs_nflx, rhs_nflx, lhs_E, rhs_E, vmax)
+
+ @returns nflx_rho, nflx_qx, nflx_qy, nflx_E
+end
+
+
@with_signature function maxspeed(equation::EulerEquation2d)
@accepts rho, qx, qy, E
@unpack eos = equation
q = sqrt(qx^2+qy^2)
- eps, = internal_energy((rho,q,E), equation)
+ eps, = internal_energy((rho,qx,qy,E), equation)
p = eos(Pressure, Density, rho, InternalEnergy, eps)
v = abs(q / p) + sqrt(abs(eos.gamma * p / rho))
@returns v
@@ -294,82 +384,50 @@ end
@returns flx_S_x, flx_S_y
end
-#######################################################################
-# LDG implementation #
-#######################################################################
-
-
-@with_signature function bak_ldg_flux(equation::EulerEquation)
- @accepts State(ldg_rho, ldg_q, ldg_E), rho, q, E
- flx_ldg_rho, flx_ldg_q, flx_ldg_E = rho, q, E
- @returns flx_ldg_rho, flx_ldg_q, flx_ldg_E
-end
-@with_signature function ldg_nflux(eq::EulerEquation)
- @accepts Prefix(lhs), rho, q, E
- @accepts Prefix(rhs), rho, q, E
- @accepts nx
-
- nflx_rho = 0.5 * nx * (lhs_rho + rhs_rho)
- nflx_q = 0.5 * nx * (lhs_q + rhs_q)
- nflx_E = 0.5 * nx * (lhs_E + rhs_E)
-
- @returns nflx_rho, nflx_q, nflx_E
-end
-@with_signature function ldg_bdryllf(eq::EulerEquation)
- @accepts Prefix(lhs), rho, q, E
- @accepts Prefix(rhs), rho, q, E
- @accepts nx
-
- # set exterior state to 0 according to stability analysis
- # see dg-notes/tex/evm.pdf for reasoning
- nflx_rho = nx * lhs_rho
- nflx_q = nx * lhs_q
- nflx_E = nx * lhs_E
- @returns nflx_rho, nflx_q, nflx_E
-end
+@with_signature function ldg_nflux(eq::EulerEquation2d)
+ @accepts Prefix(lhs), rho, qx, qy, E
+ @accepts Prefix(rhs), rho, qx, qy, E
+ @accepts nx, ny
+ avg_rho = 0.5 * (lhs_rho + rhs_rho)
+ avg_qx = 0.5 * (lhs_qx + rhs_qx)
+ avg_qy = 0.5 * (lhs_qy + rhs_qy)
+ avg_E = 0.5 * (lhs_E + rhs_E)
-@with_signature function ldg_speed(equation::EulerEquation)
- @accepts State(ldg_rho, ldg_q, ldg_E), rho, q, E
- v_ldg_rho, v_ldg_q, v_ldg_E = 1, 1, 1
- @returns v_ldg_rho, v_ldg_q, v_ldg_E
-end
+ nflx_rho_x = nx * avg_rho
+ nflx_rho_y = ny * avg_rho
+ nflx_qx_x = nx * avg_qx
+ nflx_qx_y = ny * avg_qx
-# TODO Remove this together with ldg_flux, ldg_speed, av_speed when
-# moving everything to the new broadcast interface
-@with_signature function bak_av_flux(equation::EulerEquation)
- @accepts State(ldg_rho, ldg_q, ldg_E), smoothed_mu
- flx_ldg_rho, flx_ldg_q, flx_ldg_E =
- smoothed_mu * ldg_rho, smoothed_mu * ldg_q, smoothed_mu * ldg_E
- @returns flx_ldg_rho, flx_ldg_q, flx_ldg_E
-end
+ nflx_qy_x = nx * avg_qy
+ nflx_qy_y = ny * avg_qy
+ nflx_E_x = nx * avg_E
+ nflx_E_y = ny * avg_E
-@with_signature function av_speed(equation::EulerEquation)
- @accepts State(ldg_rho, ldg_q, ldg_E), smoothed_mu
- v_av_rho, v_av_q, v_av_E = smoothed_mu, smoothed_mu, smoothed_mu
- @returns v_av_rho, v_av_q, v_av_E
+ @returns nflx_rho_x, nflx_rho_y, nflx_qx_x, nflx_qx_y,
+ nflx_qy_x, nflx_qy_y, nflx_E_x, nflx_E_y
end
-
-
-@with_signature function ldg_nflux(eq::EulerEquation2d)
- @accepts Prefix(lhs), rho, qx, qy, E
- @accepts Prefix(rhs), rho, qx, qy, E
+@with_signature function ldg_bdryllf(eq::EulerEquation2d)
+ @accepts Prefix(in), rho, qx, qy, E
+ @accepts Prefix(out), rho, qx, qy, E
@accepts nx, ny
- nflx_rho_x = 0.5 * nx * (lhs_rho + rhs_rho)
- nflx_rho_y = 0.5 * ny * (lhs_rho + rhs_rho)
+ # set exterior state to 0 according to stability analysis
+ # see dg-notes/tex/evm.pdf for reasoning
+ nflx_rho_x = nx * in_rho
+ nflx_rho_y = ny * in_rho
- nflx_qx_x = 0.5 * nx * (lhs_qx + rhs_qx)
- nflx_qx_y = 0.5 * ny * (lhs_qx + rhs_qx)
+ nflx_qx_x = nx * in_qx
+ nflx_qx_y = ny * in_qx
- nflx_qy_x = 0.5 * nx * (lhs_qy + rhs_qy)
- nflx_qy_y = 0.5 * ny * (lhs_qy + rhs_qy)
+ nflx_qy_x = nx * in_qy
+ nflx_qy_y = ny * in_qy
- nflx_E_x = 0.5 * nx * (lhs_E + rhs_E)
- nflx_E_y = 0.5 * ny * (lhs_E + rhs_E)
+ nflx_E_x = nx * in_E
+ nflx_E_y = ny * in_E
@returns nflx_rho_x, nflx_rho_y, nflx_qx_x, nflx_qx_y,
nflx_qy_x, nflx_qy_y, nflx_E_x, nflx_E_y
@@ -379,25 +437,29 @@ end
@with_signature function av_flux_2d(eq::EulerEquation2d)
@accepts flx_rho_x, flx_rho_y, ldg_rho_x, ldg_rho_y,
flx_qx_x, flx_qx_y, ldg_qx_x, ldg_qx_y,
- flx_qy_x, flx_qy_y, ldg_qy_x, ldg_qy_y,
+ #=flx_qy_x,=# flx_qy_y, ldg_qy_x, ldg_qy_y,
flx_E_x, flx_E_y, ldg_E_x, ldg_E_y,
smoothed_mu
+ flx_qy_x = flx_qx_y
+
flx_rho_x += smoothed_mu * ldg_rho_x
flx_rho_y += smoothed_mu * ldg_rho_y
- flx_qx_x += smoothed_mu * ldg_qx_x
- flx_qx_y += smoothed_mu * ldg_qx_y
+ flx_qx_x += smoothed_mu * ldg_qx_x
+ flx_qx_y += smoothed_mu * ldg_qx_y
- flx_qy_x += smoothed_mu * ldg_qy_x
- flx_qy_y += smoothed_mu * ldg_qy_y
+ flx_qy_x += smoothed_mu * ldg_qy_x
+ flx_qy_y += smoothed_mu * ldg_qy_y
- flx_E_x += smoothed_mu * ldg_E_x
- flx_E_y += smoothed_mu * ldg_E_y
+ flx_E_x += smoothed_mu * ldg_E_x
+ flx_E_y += smoothed_mu * ldg_E_y
@returns flx_rho_x, flx_rho_y, flx_qx_x, flx_qx_y,
flx_qy_x, flx_qy_y, flx_E_x, flx_E_y
end
+
+
@with_signature function av_nflux_2d(eq::EulerEquation2d)
@accepts Prefix(lhs), ldg_rho_x, ldg_rho_y, ldg_qx_x, ldg_qx_y, ldg_qy_x, ldg_qy_y, ldg_E_x, ldg_E_y, smoothed_mu
@accepts Prefix(rhs), ldg_rho_x, ldg_rho_y, ldg_qx_x, ldg_qx_y, ldg_qy_x, ldg_qy_y, ldg_E_x, ldg_E_y, smoothed_mu
diff --git a/src/EulerEq/initialdata.jl b/src/EulerEq/initialdata.jl
index 1b34bccc..7ba9dd2a 100644
--- a/src/EulerEq/initialdata.jl
+++ b/src/EulerEq/initialdata.jl
@@ -136,10 +136,8 @@ function initialdata_equation_of_state!(env, P, eos::EquationOfState.IdealGas, p
Ï0, vx0, vy0, p0 = [ similar(x) for _ = 1:4 ]
CI = CartesianIndices((Npts*Kx,Npts*Ky))
for (i,(xi,yi)) in enumerate(zip(x,y))
- ix, _ = Tuple(CI[i])
ylow = 0.25 + Y(xi)
yup = 0.75 + Y(xi)
- # if yi >= ylow[ix] && yi <= yup[ix]
if yi >= ylow && yi <= yup
Ï0[i] = 2.0
vx0[i] = -0.5
@@ -190,13 +188,20 @@ end
function initialdata_hrsc!(env, P, hrsc::Maybe{HRSC.AbstractHRSC}, prms, mesh) end
-function initialdata_hrsc!(env, P, hrsc::HRSC.SmoothedArtificialViscosity, prms, mesh)
+function initialdata_hrsc!(env, P, hrsc::HRSC.SmoothedArtificialViscosity, prms, mesh::Mesh1d)
@unpack cache = env
@unpack mu = get_cell_variables(mesh.cache)
@unpack smoothed_mu = get_static_variables(mesh.cache)
initialdata_hrsc!(env, P, hrsc.av, prms, mesh)
hrsc.smoother(smoothed_mu, mu, env.mesh, false)
end
+function initialdata_hrsc!(env, P, hrsc::HRSC.SmoothedArtificialViscosity, prms, mesh::Mesh2d)
+ @unpack cache = env
+ @unpack mu = get_cell_variables(mesh.cache)
+ @unpack smoothed_mu = get_static_variables(mesh.cache)
+ initialdata_hrsc!(env, P, hrsc.av, prms, mesh)
+ hrsc.smoother(smoothed_mu, mu, env.mesh, true, true)
+end
function initialdata_hrsc!(env, P::Project, hrsc::HRSC.EntropyViscosity, prms, ::Mesh1d)
@unpack mu = get_cell_variables(env.cache)
@@ -265,7 +270,7 @@ function initialdata_tci!(env, P::Project, tci::TCI.AbstractTCI, prms)
end
-function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh)
+function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh::Mesh1d)
@unpack rho, q, E = get_dynamic_variables(mesh.cache)
@unpack flx_rho, flx_q, flx_E, v = get_static_variables(mesh.cache)
@unpack init_flx_rho, init_flx_q, init_flx_E, init_v = get_static_variables(mesh.cache)
@@ -280,6 +285,17 @@ function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mes
init_flx_E .= flx_E
init_v .= v
end
+function initialdata_bdrycond!(env, P::Project, bdrycond::dg1d.DirichletBC2, mesh::Mesh2d)
+ @unpack rho, qx, qy, E = get_dynamic_variables(mesh.cache)
+ @unpack init_rho, init_qx, init_qy, init_E, v, init_v = get_static_variables(mesh.cache)
+ broadcast_volume_2!(speed, P.equation, mesh)
+ broadcast_volume_2!(flux, P.equation, mesh)
+ init_rho .= rho
+ init_qx .= qx
+ init_qy .= qy
+ init_E .= E
+ init_v .= v
+end
function initialdata_smooth!(env, P, mesh)
diff --git a/src/EulerEq/rhs.jl b/src/EulerEq/rhs.jl
index a46e85c7..d2c5808b 100644
--- a/src/EulerEq/rhs.jl
+++ b/src/EulerEq/rhs.jl
@@ -179,7 +179,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
## solve auxiliary equation: q + ∂x u = 0
broadcast_faces_2!(ldg_nflux, equation, mesh)
- dg1d.broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
+ broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
compute_rhs_weak_form!(ldg_rho, rho, nflx_rho, mesh)
compute_rhs_weak_form!(ldg_q, q, nflx_q, mesh)
compute_rhs_weak_form!(ldg_E, E, nflx_E, mesh)
@@ -189,7 +189,7 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
broadcast_faces_2!(llf, equation, mesh)
broadcast_volume_2!(av_flux, equation, mesh)
broadcast_faces_2!(av_nflux, equation, mesh)
- dg1d.broadcast_bdry_2!(bdryllf, equation, P.bdrycond, mesh)
+ broadcast_bdry_2!(bdryllf, equation, P.bdrycond, mesh)
# TODO don't we also need av_bdrylff here?
compute_rhs_weak_form!(rhs_rho, flx_rho, nflx_rho, mesh)
@@ -219,7 +219,8 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
nflx_qx_x, nflx_qx_y,
nflx_qy_x, nflx_qy_y,
nflx_E_x, nflx_E_y = get_bdry_variables(mesh.cache)
- @unpack t = get_global_variables(mesh.cache)
+
+ broadcast_volume_2!(pressure, equation, mesh)
## TODO Maybe implement a macro with this syntax
# @weakform(mesh) ∂t ldg_rho_x + ∂x rho = [ nflx_rho_x ]
@@ -234,7 +235,8 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
## solve auxiliary equation: q + ∂x u = 0
broadcast_faces_2!(ldg_nflux, equation, mesh)
- # broadcast_boundaryconditions!(central_flux, ldg_bdryconds, cache, mesh, 0.0)
+ # TODO Need to revise Bdry conditions
+ # broadcast_bdry_2!(ldg_bdryllf, equation, P.bdrycond, mesh)
compute_rhs_weak_form!(ldg_rho_x, rho, 0, nflx_rho_x, mesh)
compute_rhs_weak_form!(ldg_rho_y, 0, rho, nflx_rho_y, mesh)
compute_rhs_weak_form!(ldg_qx_x, qx, 0, nflx_qx_x, mesh)
@@ -247,11 +249,12 @@ function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::M
## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
broadcast_volume_2!(flux, equation, mesh)
broadcast_faces_2!(llf, equation, mesh)
+ # TODO Need to revise Bdry conditions
+ # broadcast_bdry_2!(bdryllf, equation, mesh)
- # broadcast_boundaryconditions!(lax_friedrich_flux, bdryconds, cache, mesh, 0.0)
broadcast_volume_2!(av_flux_2d, equation, mesh)
broadcast_faces_2!(av_nflux_2d, equation, mesh)
- # broadcast_boundaryconditions!(mean_flux, av_bdryconds, cache, mesh, 0.0)
+ # TODO don't we also need av_bdrylff here?
compute_rhs_weak_form!(rhs_rho, flx_rho_x, flx_rho_y, nflx_rho, mesh)
compute_rhs_weak_form!(rhs_qx, flx_qx_x, flx_qx_y, nflx_qx, mesh)
diff --git a/src/EulerEq/setup.jl b/src/EulerEq/setup.jl
index 85da1d33..297ba3c5 100644
--- a/src/EulerEq/setup.jl
+++ b/src/EulerEq/setup.jl
@@ -25,7 +25,7 @@ function Project(env::Environment, mesh::Mesh1d, prms)
# TODO Somehow replace _register_variables! with register_variables!
@unpack cache = env
_register_variables!(mesh, P)
- _register_variables!(mesh.cache, bdrycond)
+ _register_variables!(mesh, bdrycond)
register_variables!(mesh, rsolver, dont_register=true)
display(cache)
@@ -65,6 +65,7 @@ function Project(env::Environment, mesh::Mesh2d, prms)
equation = EulerEquation2d(eos)
ldg_rsolver, av_rsolver = nothing, nothing
+ bdrycond = dg1d.DirichletBC2()
P = Project2d(equation, nothing, nothing,
nothing, nothing, hrsc, nothing, nothing, nothing)
@@ -72,6 +73,7 @@ function Project(env::Environment, mesh::Mesh2d, prms)
# TODO Somehow replace _register_variables! with register_variables!
@unpack cache = env
_register_variables!(mesh, P, mesh)
+ _register_variables!(mesh, bdrycond)
_register_variables!(mesh, hrsc)
display(cache)
@@ -187,9 +189,14 @@ end
# TODO Already overwritten above
# function _register_variables!(cache, bdryconds::Nothing) end
-function _register_variables!(cache, bdryconds::dg1d.AbstractBC)
- register_variables!(cache,
+function _register_variables!(mesh::Mesh1d, bdryconds::dg1d.AbstractBC)
+ register_variables!(mesh.cache,
static_variablenames = (:init_rho, :init_q, :init_E, :init_v,
:init_flx_rho, :init_flx_q, :init_flx_E)
)
end
+function _register_variables!(mesh::Mesh2d, bdryconds::dg1d.AbstractBC)
+ register_variables!(mesh.cache,
+ static_variablenames = (:init_rho, :init_qx, :init_qy, :init_E, :init_v)
+ )
+end
diff --git a/src/SRHD/equation.jl b/src/SRHD/equation.jl
index 6c7740c7..400d54c3 100644
--- a/src/SRHD/equation.jl
+++ b/src/SRHD/equation.jl
@@ -312,9 +312,6 @@ end
end
-@inline LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
-
-
# TODO Remove llf_D, llf_Sx, llf_Sy, llf_tau in favor of lld below
@with_signature function llf_D(equation::Equation2d)
@accepts Prefix(lhs), D, flx_D_x, flx_D_y, max_v
diff --git a/src/ScalarEq/equation.jl b/src/ScalarEq/equation.jl
index dd56db27..652c4a2f 100644
--- a/src/ScalarEq/equation.jl
+++ b/src/ScalarEq/equation.jl
@@ -112,9 +112,6 @@ end
#######################################################################
-LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
-
-
@with_signature function flux(eq::Advection2d)
@accepts State(u)
flx_x, flx_y = eq.v_x * u, eq.v_y * u
diff --git a/src/dg1d.jl b/src/dg1d.jl
index 05e4851b..ed481f15 100644
--- a/src/dg1d.jl
+++ b/src/dg1d.jl
@@ -34,9 +34,6 @@ export # parameter handling
query, splitallext, trim_keys, match_keys, prepend_keys,
# type utils
concretetypes,
- # numeric utilities
- periodize, periodic_index, # deprecated
- estimate_convergence_rate,
# parameter files and file system operations
mk_outputdir, dirname_path, flatten, get_parameter_file, get_parameter_file_index,
select_parameter_files, gather_output_directories, gather_parameter_files,
@@ -44,10 +41,18 @@ export # parameter handling
find_duplicated_parameter_files,
run_projects, analyze_projects,
TODO, @TODO, @stringify_names,
- indent,
- absmax
+ indent
include("utils.jl")
+
+export LLF,
+ absmax,
+ periodize,
+ estimate_convergence_rate,
+ # deprecated
+ periodic_index
+include("numutils.jl")
+
export Output, read_variable, read_variable!, read_mesh
include("output.jl")
@@ -86,7 +91,7 @@ export Cache, register_variables!,
get_bdry_variables, get_rhs_variables, get_static_variables,
get_variable, wrap_dynamic_variables!, wrap_rhs_variables!,
broadcast_volume!, variablenames, variablegroups,
- broadcast_volume_2!, broadcast_faces_2!
+ broadcast_volume_2!, broadcast_faces_2!, broadcast_bdry_2!
include("cache.jl")
export Mesh, Mesh1d, Mesh2d, MeshInterpolator, layout, n_points, n_cells, dimension,
diff --git a/src/dg_rhs.jl b/src/dg_rhs.jl
index 727330c0..26841c1e 100644
--- a/src/dg_rhs.jl
+++ b/src/dg_rhs.jl
@@ -73,29 +73,37 @@ function compute_rhs_weak_form!(rhs, fx, fy::Real, nf, mesh::Mesh2d)
@unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
Ml_down, Ml_up = Ml_lhs, Ml_rhs
Nx, Ny = Npts, Npts
- @inbounds for k = 1:n_cells(mesh.tree)
+ @inbounds for kk = 1:n_cells(mesh.tree)
# bulk
- idxs_bulk = cellindices(mesh, k)
+ idxs_bulk = cellindices(mesh, kk)
v_rhs = view(rhs, idxs_bulk)
v_fx = view(fx, idxs_bulk)
v_dxdX = view(dxdX, idxs_bulk)
+ v_dydY = view(dydY, idxs_bulk)
# faces
- idxs_lhs, idxs_rhs, _, _ = faceindices(mesh, k)
+ idxs_lhs, idxs_rhs, idxs_down, idxs_up = faceindices(mesh, kk)
v_nf_lhs = view(nf, idxs_lhs)
v_nf_rhs = view(nf, idxs_rhs)
+ v_nf_down = view(nf, idxs_down)
+ v_nf_up = view(nf, idxs_up)
# TODO Benchmark and see whether LoopVectorization can used here.
@inbounds for j=1:Ny, i=1:Nx
rhsij = 0
dxdXij = v_dxdX[i,j]
+ dydYij = v_dydY[i,j]
for k = 1:Nx
rhsij += MDM[i,k] * v_fx[k,j]
end
v_rhs[i,j] = rhsij * dxdXij
v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j] * Ml_lhs[i]) * dxdXij
+ # TODO Can we really neglect num fluxes in orthogonal direction in general?
+ # Or we need to adjust for triangular meshes?
+ # v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j] * Ml_lhs[i]) * dxdXij +
+ # (v_nf_up[i] * Ml_up[j] + v_nf_down[i] * Ml_down[j]) * dydYij
end
end
@@ -107,29 +115,37 @@ function compute_rhs_weak_form!(rhs, fx::Real, fy, nf, mesh::Mesh2d)
@unpack invM, MDM, Ml_lhs, Ml_rhs, Npts = element
Ml_down, Ml_up = Ml_lhs, Ml_rhs
Nx, Ny = Npts, Npts
- @inbounds for k = 1:n_cells(mesh.tree)
+ @inbounds for kk = 1:n_cells(mesh.tree)
# bulk
- idxs_bulk = cellindices(mesh, k)
+ idxs_bulk = cellindices(mesh, kk)
v_rhs = view(rhs, idxs_bulk)
v_fy = view(fy, idxs_bulk)
+ v_dxdX = view(dxdX, idxs_bulk)
v_dydY = view(dydY, idxs_bulk)
# faces
- _, _, idxs_down, idxs_up = faceindices(mesh, k)
+ idxs_lhs, idxs_rhs, idxs_down, idxs_up = faceindices(mesh, kk)
+ v_nf_lhs = view(nf, idxs_lhs)
+ v_nf_rhs = view(nf, idxs_rhs)
v_nf_down = view(nf, idxs_down)
v_nf_up = view(nf, idxs_up)
# TODO Benchmark and see whether LoopVectorization can used here.
@inbounds for j=1:Ny, i=1:Nx
rhsij = 0
+ dxdXij = v_dxdX[i,j]
dydYij = v_dydY[i,j]
for k = 1:Ny
rhsij += MDM[j,k] * v_fy[i,k]
end
- v_rhs[i,j] += rhsij * dydYij
+ v_rhs[i,j] = rhsij * dydYij
v_rhs[i,j] -= (v_nf_up[i] * Ml_up[j] + v_nf_down[i] * Ml_down[j]) * dydYij
+ # TODO Can we really neglect num fluxes in orthogonal direction in general?
+ # Or we need to adjust for triangular meshes?
+ # v_rhs[i,j] -= (v_nf_rhs[j] * Ml_rhs[i] + v_nf_lhs[j] * Ml_lhs[i]) * dxdXij +
+ # (v_nf_up[i] * Ml_up[j] + v_nf_down[i] * Ml_down[j]) * dydYij
end
end
diff --git a/src/numutils.jl b/src/numutils.jl
new file mode 100644
index 00000000..d803d82d
--- /dev/null
+++ b/src/numutils.jl
@@ -0,0 +1,106 @@
+"""
+ LLF(nfl,nfr,ul,ur,v)
+
+Local Lax-Friedrich flux:
+
+ nf* = 1/2 ( f(ul) + f(ur) - v * (ur - ul) )
+"""
+@inline LLF(nfl,nfr,ul,ur,v) = 1/2 * (nfl + nfr - v * (ur - ul))
+
+
+"""
+ estimate_convergence_rate(n1, n2, f1, f2)
+
+Experimental rate of convergence based on resolutions `n1, n2` and data `f1, f2`.
+"""
+function estimate_convergence_rate(n1, n2, f1, f2)
+ return log10(f1/f2) / log10(n2/n1)
+end
+
+
+"""
+ purge_zeros!(M, thrs=1e-13)
+
+Scan `M` for elements smaller than `thrs` and set them to 0.
+"""
+purge_zeros!(M, thrs=1e-13) = M[abs.(M).<thrs] .= 0.0
+
+
+"""
+ least_sqr_linear(x, y)
+
+Compute a linear least squares fit for data `x, y`.
+Returns fit function `f, a, b`, where `f(x) = a + b * x` and `a,b` are the fitting coefficients.
+Alogrithm taken from https://mathworld.wolfram.com/LeastSquaresFitting.html
+"""
+function least_sqr_linear(x, y)
+
+ n = length(x)
+ @assert n == length(y) "Mismatch in number of data points x, y!"
+
+ yx = dot(y, x)
+ xx = dot(x, x)
+ mean_y = sum(y) / n
+ mean_x = sum(x) / n
+
+ ss_xx = xx - mean_x^2 * n
+ ss_xy = yx - mean_x * mean_y * n
+
+ b = ss_xy / ss_xx
+ a = mean_y - b * mean_x
+
+ f(x) = a + b * x
+ return f, a, b
+end
+
+
+stepdata(x, x0, fL, fR) = [ (xi < x0) ? fL : fR for xi in x ]
+
+
+"""
+ periodize(i, N)
+
+Wrap `i` around period `N` if `i > N`, otherwise return `i`.
+Behaves similar when called with `Int`s or `Real`s.
+"""
+@inline function periodize(i::T, N::S) where {T<:Int,S<:Int}
+ mod = i % N
+ return mod <= 0 ? N+mod : mod
+end
+@inline function periodize(x::T, L::S) where {T<:Real,S<:Real}
+ xround = floor(x/L) * L
+ return x <= xround ? x : x - xround
+end
+
+
+@deprecate periodic_index(i, N) periodize(i, N)
+
+
+@inline absolute_maximum(x) = norm(x, Inf)
+
+
+"""
+ absolute_maximum_index(x::AbstractArray)
+
+Return index and value of the absolute maximum of array `x`.
+"""
+function absolute_maximum_index(x::AbstractArray)
+ max_i = 1
+ max_x = abs(x[1])
+ for (i,xi) in enumerate(x)
+ abs_x = abs(xi)
+ if abs_x > max_x
+ max_i = i
+ max_x = abs_x
+ end
+ end
+ return max_i, max_x
+end
+
+
+"""
+ absmax(vs...)
+
+Return absolute maximum of `vs...`.
+"""
+@inline absmax(vs...) = mapreduce(abs, max, vs)
diff --git a/src/utils.jl b/src/utils.jl
index 64792915..088b1485 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,21 +1,3 @@
-"""
- estimate_convergence_rate(n1, n2, f1, f2)
-
-Experimental rate of convergence based on resolutions `n1, n2` and data `f1, f2`.
-"""
-function estimate_convergence_rate(n1, n2, f1, f2)
- return log10(f1/f2) / log10(n2/n1)
-end
-
-
-"""
- purge_zeros!(M, thrs=1e-13)
-
-Scan `M` for elements smaller than `thrs` and set them to 0.
-"""
-purge_zeros!(M, thrs=1e-13) = M[abs.(M).<thrs] .= 0.0
-
-
function get_outputdir()
thisdir = dirname(@__FILE__)
outputdir = "$thisdir/../outputs/"
@@ -30,37 +12,6 @@ end
@deprecate get_output_dir() get_outputdir()
-"""
- least_sqr_linear(x, y)
-
-Compute a linear least squares fit for data `x, y`.
-Returns fit function `f, a, b`, where `f(x) = a + b * x` and `a,b` are the fitting coefficients.
-Alogrithm taken from https://mathworld.wolfram.com/LeastSquaresFitting.html
-"""
-function least_sqr_linear(x, y)
-
- n = length(x)
- @assert n == length(y) "Mismatch in number of data points x, y!"
-
- yx = dot(y, x)
- xx = dot(x, x)
- mean_y = sum(y) / n
- mean_x = sum(x) / n
-
- ss_xx = xx - mean_x^2 * n
- ss_xy = yx - mean_x * mean_y * n
-
- b = ss_xy / ss_xx
- a = mean_y - b * mean_x
-
- f(x) = a + b * x
- return f, a, b
-end
-
-
-stepdata(x, x0, fL, fR) = [ (xi < x0) ? fL : fR for xi in x ]
-
-
# stolen from https://discourse.julialang.org/t/convert-type-to-dictionary/6982/7
typedict(x::T) where T = Dict(fn=>getfield(x,fn) for fn ∈ fieldnames(T))
@@ -78,45 +29,6 @@ end
-"""
- periodize(i, N)
-
-Wrap `i` around period `N` if `i > N`, otherwise return `i`.
-Behaves similar when called with `Int`s or `Real`s.
-"""
-@inline function periodize(i::T, N::S) where {T<:Int,S<:Int}
- mod = i % N
- return mod <= 0 ? N+mod : mod
-end
-@inline function periodize(x::T, L::S) where {T<:Real,S<:Real}
- xround = floor(x/L) * L
- return x <= xround ? x : x - xround
-end
-
-
-@deprecate periodic_index(i, N) periodize(i, N)
-
-
-@inline absolute_maximum(x) = norm(x, Inf)
-
-
-function absolute_maximum_index(x::AbstractArray)
- max_i = 1
- max_x = abs(x[1])
- for (i,xi) in enumerate(x)
- abs_x = abs(xi)
- if abs_x > max_x
- max_i = i
- max_x = abs_x
- end
- end
- return max_i, max_x
-end
-
-
-@inline absmax(vs...) = mapreduce(abs, max, vs)
-
-
function query_key(prms, key, default)
value = get(prms, key, missing)
if ismissing(value)
diff --git a/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml b/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml
index c844f086..860254bc 100644
--- a/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml
+++ b/test/IntegrationTests/refs/euler2d_kelvin_helmholtz_entropyav.toml
@@ -1,5 +1,6 @@
[EulerEq]
id = "kelvin_helmholtz"
+bc = "periodic"
[EquationOfState]
eos = "idealgas"
@@ -8,12 +9,13 @@ idealgas_gamma = 1.4
[Mesh]
range = [ 0.0,1.0, 0.0,1.0 ]
cfl = 0.1
+# n = [4,4]
n = [1,1]
-k = [256,256]
+k = [64,64]
basis = "lgl"
+periodic = [true, true]
[Output]
-every_iteration = 1
every_dt = 0.02
variables = [ "rho", "qx", "qy", "E", "p", "smoothed_mu" ]
enable1d = false
@@ -26,3 +28,5 @@ tend = 2.0
[HRSC]
method = "av"
av_method = "entropy"
+entropy_cmax = 0.1
+entropy_ce = 0.1
--
GitLab