dragonfly_sim.utils.Euler_utils
Functions
|
|
|
|
|
|
|
HLL (Harten-Lax-van Leer) approximate Riemann solver, a drop-in |
|
HLLC approximate Riemann solver: HLL with the contact wave restored. |
|
Local Lax-Friedrichs (LLF) numerical flux, as implemented in dolfin_dg. |
|
|
|
|
|
|
|
|
|
The signed extremal acoustic wave speeds (u.n - c, u.n + c) of the |
|
|
|
|
|
Module Contents
- dragonfly_sim.utils.Euler_utils.boundary_flux(U, U_ext, n, gamma)
- dragonfly_sim.utils.Euler_utils.energy_density(p_in, rho_in, u_in, gamma)
- dragonfly_sim.utils.Euler_utils.flux(U, gamma)
- dragonfly_sim.utils.Euler_utils.hll_flux(U_p, U_m, n, gamma)
HLL (Harten-Lax-van Leer) approximate Riemann solver, a drop-in alternative to llf_flux for flows containing shocks.
The Riemann fan is bounded by two waves with Davis speed estimates
SL = min(un_p - c_p, un_m - c_m) SR = max(un_p + c_p, un_m + c_m)
and the flux is the corresponding single averaged intermediate state
SL >= 0 : F(U_p).n (fully upwind) SR <= 0 : F(U_m).n (fully upwind) else : (SR F(U_p).n - SL F(U_m).n + SL SR (U_m - U_p))
/ (SR - SL)
Written branchlessly below via SL_m = min(SL, 0), SR_p = max(SR, 0), which reproduces all three cases exactly: SL >= 0 sends SL_m to 0 and collapses the quotient to F(U_p).n, SR <= 0 sends SR_p to 0 and collapses it to F(U_m).n. This is not merely a tidier spelling. It needs only min_value/max_value rather than nested ufl.conditional nodes, so the expression is genuinely continuous across the branch boundaries instead of only piecewise defined – which is what keeps the assembled Jacobian consistent with the residual for a scheme whose branch boundaries (un = +/- c) sit exactly on the sonic lines of a transonic solution. That matters twice over here, since this model’s Jacobian also feeds the adjoint used for shape optimization.
Verified numerically against an explicit three-branch reference in test_algorithms/test_hll_flux.py – branchless vs branched agree to ~1e-16 with states swept across sub-, trans- and supersonic so that every branch is exercised, in both 2D and 3D.
- dragonfly_sim.utils.Euler_utils.hllc_flux(U_p, U_m, n, gamma)
HLLC approximate Riemann solver: HLL with the contact wave restored.
hll_flux collapses the entire Riemann fan into a single averaged intermediate state, which throws the contact wave away entirely. HLLC reinstates it as a third wave of speed S*, so the two intermediate states either side of the contact are resolved separately. Wave speeds SL/SR come from the same signed_wave_speeds/Davis estimate hll_flux uses – deliberately, so the two schemes differ only in the contact treatment – and the contact speed is
- S* = (p_m - p_p + rho_p un_p (SL - un_p)
rho_m un_m (SR - un_m))
/ (rho_p (SL - un_p) - rho_m (SR - un_m))
with the star state on side K (Toro eq. 10.39; E is the volumetric total energy, as primitives() returns it)
- U*_K = rho_K (S_K - un_K)/(S_K - S*) *
- [ 1,
u_K + (S* - un_K) n, E_K/rho_K + (S* - un_K)(S* + p_K/(rho_K (S_K - un_K))) ]
- dragonfly_sim.utils.Euler_utils.llf_flux(U_p, U_m, n, gamma)
Local Lax-Friedrichs (LLF) numerical flux, as implemented in dolfin_dg. Mathematically, this evaluates to the same expression as the Rusanov flux for the compressible Euler equations.
- dragonfly_sim.utils.Euler_utils.llf_penalty_term(U_p, U_m, n, gamma)
- dragonfly_sim.utils.Euler_utils.normal_flux(U, n_normal, gamma)
- dragonfly_sim.utils.Euler_utils.pressure(U, gamma)
- dragonfly_sim.utils.Euler_utils.primitives(U)
- dragonfly_sim.utils.Euler_utils.signed_wave_speeds(U, n, gamma)
The signed extremal acoustic wave speeds (u.n - c, u.n + c) of the normal flux Jacobian.
The Lax-Friedrichs family only ever needs the largest wave speed in magnitude, |u.n| + c, whereas an approximate Riemann solver needs to know which direction the fan is travelling in – that is precisely what lets it become fully one-sided in supersonic flow instead of always dissipating both ways. So this returns the signed quantity that the magnitude throws away.
Floors rho/p at 1e-8 before c = sqrt(gamma*p/rho). Worth being explicit about why this floor is not merely cosmetic here: sqrt of a negative argument returns NaN, not just an implausibly large number, so a single transiently negative pointwise pressure – which a DG discretization can readily produce in early, far-from-converged Newton iterations – would otherwise poison the entire assembled residual rather than just the facet it occurred on.
Works for facet-restricted traces (U(‘+’)/U(‘-‘)) exactly like primitives()/pressure().
- dragonfly_sim.utils.Euler_utils.slip_wall_state(U, n)
- dragonfly_sim.utils.Euler_utils.subsonic_inflow_state(U, rho_in, u_in, gamma)
- dragonfly_sim.utils.Euler_utils.subsonic_outflow_state(U, p_out, gamma)