compressible Euler

intsharp/solvers/euler_1d.py, euler_5eq_1d.py, euler_5eq_2d.py, euler_dg_1d.py | intsharp/flux_ausm.py, flux_hllc.py | intsharp/eos.py, exact_riemann_1d.py | intsharp/limiters.py

Overview

The compressible Euler solver implements inviscid compressible-flow schemes in 1D and 2D with finite-volume (FV), plus a 1D discontinuous Galerkin (DG) path.

FV single-phase (3 equations): Mass, momentum, and energy conservation for a single fluid.
FV two-phase (5 equations): Allaire–Clerc–Kokh / Kapila et al. model with pressure equilibrium, tracking partial densities, mixture momentum/energy, and volume fraction.
FV two-phase 2D RTI: Directionally split 2D extension of the 5-equation model with modular x/y boundary conditions and optional gravity source term.
DG single-phase (3 equations): Element-local polynomial solution (P1/P2/P3) with interface Riemann fluxes and positivity-preserving controls.

All variants use the stiffened gas equation of state. FV supports AUSM+UP and HLLC (with two-phase v_riem form), while DG currently uses HLLE. FV supports optional MUSCL reconstruction with Barth–Jespersen limiting. DG supports polynomial orders P1/P2/P3 and SSP-RK3 time stepping.

PDE System

Single-phase (3 equations)

Conservation of mass, momentum, and energy:

\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} = 0

\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} = 0

\frac{\partial E}{\partial t} + \frac{\partial \bigl((E+p)u\bigr)}{\partial x} = 0

Conserved variables: \(\mathbf{U} = (\rho,\, \rho u,\, E)\).

Total energy per unit volume: \(E = \rho e + \tfrac{1}{2}\rho u^2\), where \(e\) is specific internal energy.

Two-phase 5-equation model

Allaire–Clerc–Kokh / Kapila et al. model with pressure equilibrium (\(p_1 = p_2 = p\)):

\frac{\partial (\alpha_1\rho_1)}{\partial t} + \frac{\partial (\alpha_1\rho_1 u)}{\partial x} = 0

\frac{\partial (\alpha_2\rho_2)}{\partial t} + \frac{\partial (\alpha_2\rho_2 u)}{\partial x} = 0

\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} = 0

\frac{\partial E}{\partial t} + \frac{\partial \bigl((E+p)u\bigr)}{\partial x} = 0

\frac{\partial \alpha_1}{\partial t} + u\,\frac{\partial \alpha_1}{\partial x} = 0 \quad \text{(non-conservative)}

Conserved variables: \((\alpha_1\rho_1,\, \alpha_2\rho_2,\, \rho u,\, E,\, \alpha_1)\).

Mixture density: \(\rho = \alpha_1\rho_1 + \alpha_2\rho_2\), with \(\alpha_2 = 1 - \alpha_1\).

Equation of State (Stiffened Gas)

Both phases use the stiffened gas EOS:

p = (\gamma - 1)\,\rho e - \gamma p_\infty

Specific internal energy:

e = \frac{p + \gamma p_\infty}{(\gamma - 1)\rho}

Sound speed:

c = \sqrt{\frac{\gamma(p + p_\infty)}{\rho}}

Total energy from primitives:

E = \rho e + \tfrac{1}{2}\rho u^2 = \frac{p + \gamma p_\infty}{\gamma - 1} + \tfrac{1}{2}\rho u^2

Specific total enthalpy: \(H = (E + p) / \rho\).

Typical parameters:

Ideal gas (air): \(\gamma = 1.4\), \(p_\infty = 0\).
Water: \(\gamma \approx 4.4\), \(p_\infty \approx 6\times 10^8\) Pa.

Decoding (Conserved → Primitive)

Single-phase

From \(\mathbf{U} = (\rho,\, \rho u,\, E)\):

u = \frac{\rho u}{\rho}, \qquad \rho e = E - \tfrac{1}{2}\rho u^2, \qquad p = (\gamma-1)\,\rho e - \gamma p_\infty

Enthalpy for flux: \(H = (E + p) / \rho\).

Two-phase (5-eq)

From \((\alpha_1\rho_1,\, \alpha_2\rho_2,\, \rho u,\, E,\, \alpha_1)\):

\rho = \alpha_1\rho_1 + \alpha_2\rho_2, \quad u = \frac{\rho u}{\rho}, \quad \rho_1 = \frac{\alpha_1\rho_1}{\alpha_1}, \quad \rho_2 = \frac{\alpha_2\rho_2}{\alpha_2}

Pressure from mixture internal energy (pressure equilibrium):

E_{\text{int}} = E - \tfrac{1}{2}\rho u^2 = \alpha_1\,\frac{p + \gamma_1 p_{\infty,1}}{\gamma_1 - 1} + \alpha_2\,\frac{p + \gamma_2 p_{\infty,2}}{\gamma_2 - 1}

Solve for \(p\): \(E_{\text{int}} = p\,A + B\) with

A = \frac{\alpha_1}{\gamma_1 - 1} + \frac{\alpha_2}{\gamma_2 - 1}, \qquad B = \frac{\alpha_1\gamma_1 p_{\infty,1}}{\gamma_1 - 1} + \frac{\alpha_2\gamma_2 p_{\infty,2}}{\gamma_2 - 1}

Hence \(p = (E_{\text{int}} - B) / A\).

Mixture sound speed (Wood's formula):

\frac{1}{\rho c^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2} \quad \Rightarrow \quad c = \sqrt{\frac{1}{\rho\left(\frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}\right)}}

Pure phases: \(c \to c_1\) when \(\alpha_1 \approx 1\), \(c \to c_2\) when \(\alpha_1 \approx 0\).

Flux Calculation (AUSM+UP / HLLC / HLLE)

The FV solver provides modular interface fluxes. AUSM+UP (below) is the default all-speed split form for single-phase and two-phase FV; HLLC is also available for robust shock/contact capturing in FV; DG uses HLLE in the current 1D implementation.

1. Interface sound speed

a_{1/2} = \tfrac{1}{2}(c_L + c_R)

2. Mach number splitting (M₄±)

\(M_L = u_L/a_{1/2}\), \(M_R = u_R/a_{1/2}\). Fourth-order polynomials:

M_4^+(M) = \begin{cases} M_1^+ = \tfrac{1}{2}(M + |M|) & |M| \ge 1 \\[4pt] M_2^+(1 - 16\beta M_2^-) & |M| < 1 \end{cases}

M_4^-(M) = \begin{cases} M_1^- = \tfrac{1}{2}(M - |M|) & |M| \ge 1 \\[4pt] M_2^-(1 + 16\beta M_2^+) & |M| < 1 \end{cases}

where \(M_2^+ = \tfrac{1}{4}(M+1)^2\), \(M_2^- = -\tfrac{1}{4}(M-1)^2\), \(\beta = 1/8\).

Pressure diffusion (low-speed):

\bar{M}^2 = \tfrac{1}{2}(M_L^2 + M_R^2), \qquad M_p = -K_p\,\max(1 - \bar{M}^2,\, 0)\,\frac{p_R - p_L}{\rho_{\text{avg}}\,a_{1/2}^2}

M_{1/2} = M_4^+(M_L) + M_4^-(M_R) + M_p

3. Pressure splitting (P₅±)

P_5^+(M) = \begin{cases} M_1^+/M & |M| \ge 1 \\[4pt] M_2^+\bigl((2-M) - 16\alpha M M_2^-\bigr) & |M| < 1 \end{cases}

P_5^-(M) = \begin{cases} M_1^-/M & |M| \ge 1 \\[4pt] M_2^-\bigl((-2-M) + 16\alpha M M_2^+\bigr) & |M| < 1 \end{cases}

\(\alpha = 3/16\).

Velocity diffusion:

p_u = -K_u\,P_5^+(M_L)\,P_5^-(M_R)\cdot 2\rho_{\text{avg}}\,a_{1/2}\,(u_R - u_L)

p_{1/2} = P_5^+(M_L)\,p_L + P_5^-(M_R)\,p_R + p_u

Default: \(K_p = 0.25\), \(K_u = 0.75\).

4. Mass and convective fluxes

\rho_{\text{upwind}} = \begin{cases} \rho_L & M_{1/2} \ge 0 \\ \rho_R & M_{1/2} < 0 \end{cases}, \qquad \dot{m} = a_{1/2}\,M_{1/2}\,\rho_{\text{upwind}}

\mathbf{F} = \dot{m}\begin{bmatrix} 1 \\ u_{\text{upwind}} \\ H_{\text{upwind}} \end{bmatrix} + \begin{bmatrix} 0 \\ p_{1/2} \\ 0 \end{bmatrix}

\(H = (E + p)/\rho\).

HLLC/HLLE note: HLL-family methods estimate left/right wave speeds and construct a conservative face flux from wave-fan states. HLLC restores the contact wave (sharper contacts); HLLE is more dissipative but robust and is used by the current DG module.

DG Formulation (P1/P2/P3)

The 1D DG solver evolves element-local polynomial coefficients for the single-phase Euler equations using interface Riemann fluxes and SSP-RK3 time integration. For polynomial order \(p\), each conserved variable has \(p+1\) local degrees of freedom.

\vec{q}_h(x,t)\big|_{K_e}=\sum_{j=0}^{p}\vec{Q}_{e,j}(t)\,\ell_j(\xi), \quad \xi\in[-1,1]

M_e\frac{d\vec{Q}_e}{dt} =\vec{R}^{\mathrm{vol}}_e(\vec{Q}_e)+\vec{R}^{\mathrm{surf}}_e(\vec{Q}_{e-1},\vec{Q}_{e},\vec{Q}_{e+1})

P1, P2, and P3 are supported via YAML configuration. Positivity safeguards are applied to keep \(\rho>0\) and \(p>0\) in shock-dominated regions.

Slope Limiter (MUSCL + Barth–Jespersen)

Optional second-order reconstruction. Cell-centered \(q_i\) with ghost cells.

Centered gradient

\Delta q_i = \tfrac{1}{2}(q_{i+1} - q_{i-1})

Barth–Jespersen limiter

Compute \(\phi_i \in [0, 1]\) such that \(q_i \pm \tfrac{1}{2}\phi_i\,\Delta q_i\) stays within \([q_{\min},\, q_{\max}]\) of cell \(i\) and its neighbors.

q_{\min} = \min(q_{i-1},\, q_i,\, q_{i+1}), \qquad q_{\max} = \max(q_{i-1},\, q_i,\, q_{i+1})

For each face \(q_{\text{face}} = q_i \pm \tfrac{1}{2}\Delta q_i\): if \(q_{\text{face}}\) exceeds bounds, reduce \(\phi\) so the extrapolated value lies within bounds. Take the minimum over both faces.

Reconstructed interface states

q_{L,i+1/2} = q_i + \tfrac{1}{2}\,\phi_i\,\Delta q_i, \qquad q_{R,i+1/2} = q_{i+1} - \tfrac{1}{2}\,\phi_{i+1}\,\Delta q_{i+1}

Positivity enforced for \(\rho\), \(p\), \(c\). First-order: \(q_L = q_i\), \(q_R = q_{i+1}\).

Time Integration

FV path: uses configured timestepper (euler or rk4) with explicit conservative update, e.g. Forward Euler:

\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x}\bigl(\mathbf{F}_{i+1/2} - \mathbf{F}_{i-1/2}\bigr)

DG path: uses an internal SSP-RK3 integrator in euler_dg_1d.py (see ssp_rk3).

CFL condition: \(\text{CFL} = \max(|u| + c)\,\Delta t/\Delta x < 1\) (DG typically requires tighter effective limits with increasing order).

Boundary Conditions

Transmissive: Zero-gradient (copy edge values to ghost cells).
Reflective: Mirror state with reversed velocity (\(u \to -u\)).
Periodic: Wrap around (left ghost = right interior, right ghost = left interior).

In 2D two-phase runs, set boundary conditions per direction with physics.euler_bc_x and physics.euler_bc_y.

Two-Phase Extensions

1. Wood's formula at interfaces

Interface sound speeds \(c_L\), \(c_R\) use Wood's formula with phase densities \(\rho_1\), \(\rho_2\) and phase sound speeds \(c_1\), \(c_2\) from stiffened gas.

2. Flux vector splitting (v_riem)

Interface velocity from mass flux: \(v_{\text{riem}} = \dot{m}/\rho_{\text{upwind}}\).

Partial density fluxes:

F_{\alpha_1\rho_1} = (\alpha_1\rho_1)_{\text{upwind}}\,v_{\text{riem}}, \qquad F_{\alpha_2\rho_2} = (\alpha_2\rho_2)_{\text{upwind}}\,v_{\text{riem}}

Ensures \(F_{\alpha_1\rho_1} + F_{\alpha_2\rho_2} = F_\rho\).

3. Volume fraction flux

Non-conservative \(\alpha\) advection: \(F_{\alpha_1} = \alpha_{1,\text{upwind}}\,v_{\text{riem}}\).

Corrective term (dissipative, for stability at interfaces):

F_{\alpha_1} \mathrel{+}= K_\alpha\,(u_L - u_R)\,\bar{\alpha}(1 - \bar{\alpha}), \quad \bar{\alpha} = \tfrac{1}{2}(\alpha_{1,L} + \alpha_{1,R}), \quad K_\alpha = 0.1

4. Post-step α normalization

After each step: \(\alpha_1 \leftarrow (\alpha_1\rho_1) / \rho\) to enforce consistency with partial densities and \(\alpha_1 + \alpha_2 = 1\).

In 2D RTI mode, optional constant gravity is applied as an explicit source after split flux updates:

\partial_t (\rho u) = \cdots + \rho g_x,\qquad \partial_t (\rho v) = \cdots + \rho g_y,\qquad \partial_t E = \cdots + \rho (g_x u + g_y v)

Analytical Convergence Mode

The Sod setup supports analytical-reference convergence studies through dedicated YAML configuration. The workflow runs multiple resolutions and compares numerical density against the exact Riemann solution, reporting and plotting \(L_1\) errors on log-log axes for FV and DG variants (e.g., FV, DG-P1, DG-P2, DG-P3).

Implementation & File Locations

Component	File	Function / Class
Single-phase 1D solver	`intsharp/solvers/euler_1d.py`	`euler_step_1d`, `EulerState1D`
Single-phase 1D DG solver	`intsharp/solvers/euler_dg_1d.py`	`euler_dg_step_1d`, DG reconstruction/evolution utilities
5-equation two-phase 1D solver	`intsharp/solvers/euler_5eq_1d.py`	`euler_step_5eq_1d`, `FiveEqState1D`
5-equation two-phase 2D solver	`intsharp/solvers/euler_5eq_2d.py`	`euler_step_5eq_2d`, `FiveEqState2D`, `create_initial_state_rti_5eq_2d`
AUSM+UP flux	`intsharp/flux_ausm.py`	`ausm_plus_up_flux_1d`, `ausm_plus_up_flux_1d_with_v_riem`
HLLC/HLLE flux	`intsharp/flux_hllc.py`	`hllc_flux_1d`, `hlle_flux_1d`
Stiffened gas EOS	`intsharp/eos.py`	`sound_speed`, `pressure_from_energy`, `enthalpy`
Exact Sod reference	`intsharp/exact_riemann_1d.py`	Exact 1D Riemann solution helpers for convergence/plot overlays
Wood's mixture sound speed	`intsharp/eos.py`	`mixture_sound_speed_wood`
MUSCL + Barth–Jespersen	`intsharp/limiters.py`	`muscl_reconstruct_1d`, `barth_jespersen_1d`

References

Liou, M.S. (2006). A sequel to AUSM, Part II: AUSM+-up for all speeds. J. Comput. Phys. 214, 137–170.
Toro, E.F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer.
Allaire, G., Clerc, S., Kokh, S. (2002). A five-equation model for the simulation of interfaces between compressible fluids. J. Comput. Phys. 181, 577–616.
Kapila et al. (2001). Two-phase modeling of deflagration-to-detonation transition in granular materials.

Validation: sod_shock_tube_1d, water_air_shock_tube_1d, rayleigh_taylor_2d.