Runge-Kutta Discontinuous Galerkin Schemes

• Runge-Kutta Discontinuous Galerkin schemes are high-order finite element methods that solve time-dependent PDEs using discontinuous representations and explicit Runge-Kutta integration.
• They incorporate δ-mapping at discontinuous flux interfaces to unify left/right states, ensuring steady-state preservation and non-oscillatory shock resolution.
• Practical applications include nonlinear elasticity in layered media and multi-class traffic flow, demonstrating entropy consistency and high accuracy.

A Runge-Kutta Discontinuous Galerkin (RKDG) scheme is a high-order finite element method for numerically solving time-dependent partial differential equations, notably hyperbolic conservation laws. The RKDG framework applies a spatial discretization by discontinuous Galerkin (DG) finite elements, combined with strong-stability-preserving (SSP) explicit Runge-Kutta time integration. When the flux function contains spatially discontinuous coefficients—i.e., $u_t + f\big(u, \theta(x)\big)_x = 0$ where $\theta(x)$ is a piecewise constant parameter—the DG scheme requires special treatment at interfaces of discontinuity to retain steady-state preservation, high accuracy, and non-oscillatory shock resolution. The hybrid RKDG scheme augmented with $\delta$-mapping achieves these goals by unifying left/right states at discontinuous interfaces and adjusting limiters to maintain entropy and supply-demand consistency (Qiao et al., 2015).

1. Problem Formulation and Flux Discontinuity

Consider the initial-value problem for a hyperbolic conservation law with spatially variable flux parameters,

$$u_t + \partial_x f(u, \theta(x)) = 0, \qquad x \in [0, L],~t > 0,$$

where $u \in \mathbb{R}^r$ and $\theta(x) \in \mathbb{R}^q$ is piecewise constant, possibly discontinuous at mesh interfaces $x_{j+1/2}$. At such interfaces, denote the left and right limits by $\theta^-(x_{j+1/2}) = \theta_j$ and $\theta^+(x_{j+1/2}) = \theta_{j+1}$, respectively.

Classical steady or stationary-shock solutions across jumps satisfy the local Rankine–Hugoniot condition,

$$[f(u, \theta_j)] = 0 \quad\text{across each jump of }\theta,$$

ensuring conservation is retained in the presence of spatially discontinuous flux.

2. Spatial Discretization via Discontinuous Galerkin

Partition $[0, L]$ into cells $I_j = (x_{j-1/2}, x_{j+1/2})$. On each cell, approximate $u$ by polynomials of degree $\leq k$,

$$u_h(x, t)\big|_{I_j} = \sum_{l=0}^k u_j^l(t)\, \varphi_j^l(x),$$

where $\{\varphi_j^l(x)\}$ is a local basis (e.g., scaled Legendre polynomials).

The weak DG form, for each basis index $l$, reads

$$\int_{I_j} (u_h)_t\, \varphi_j^l\,dx - \int_{I_j} f(u_h, \theta)\, \partial_x \varphi_j^l\,dx + \widehat{f}_{j+1/2}\, \varphi_j^l(x_{j+1/2}^-) - \widehat{f}_{j-1/2}\, \varphi_j^l(x_{j-1/2}^+) = 0.$$

After expansion, the method yields a finite-dimensional ODE system for each coefficient,

$$\frac{d}{dt} u_j^l(t) = \frac{2l + 1}{\Delta_j} \left[ \int_{I_j} f(u_h, \theta)\, \partial_x \varphi_j^l\,dx - \left\{ (-1)^l \widehat{f}_{j-1/2} + \widehat{f}_{j+1/2} \right\} \right].$$

3. Numerical Flux Definition via δ-Mapping

Direct use of standard numerical fluxes at interfaces where $\theta$ jumps fails to preserve steady states. Instead, the $\delta$-mapping algorithm maps left/right DG traces $(u^-_{j+1/2}, u^+_{j+1/2})$ to unified states $(\delta u^-_{j+1/2}, \delta u^+_{j+1/2})$ evaluated at an intermediate parameter $\bar\theta$, defined to satisfy:

• Supply-Demand Equalization: maximize $\gamma \leq 1$ subject to

$$f(\delta u^\mp, \theta_{j+1/2}) = \gamma f(u^\mp, \theta_{j+1/2\mp 1/2}),$$

• Entropy and Eigenvalue Compatibility:

$$\lambda_l(\delta u^\mp, \theta_{j+1/2})\, \lambda_l(u^\mp, \theta_{j+1/2\mp 1/2}) \geq 0,$$

ensuring no sign change for characteristic speeds.

The mapped states are used in a classical Riemann solver,

$$\widehat{f}_{j+1/2} = \widetilde{f}(\delta_{j+1/2} u^-_{j+1/2}, \delta_{j+1/2} u^+_{j+1/2}, \theta_{j+1/2}).$$

4. Slope Limiting with δ-Mapping

To control nonphysical oscillations near shocks and discontinuities, a TVB-type minmod slope limiter is applied before each RK stage. For $k=1$,

$$u_h|_{I_j}(x) = \bar{u}_j + \frac{2u_j^1}{\Delta_j}(x - x_j),$$

with slope replaced by a δ-mapped minmod,

$$\Lambda\Pi_h^1 u_h = \bar{u}_j + \frac{2}{\Delta_j} m(u_j^1, \delta_j \bar{u}_{j+1} - \bar{u}_j, \bar{u}_j - \delta_j \bar{u}_{j-1})(x - x_j),$$

where $\delta_j$ maps neighboring cell averages to $\theta_j$, guaranteeing steady-state preservation.

5. SSP Runge–Kutta Time Integration

Let the semi-discrete DG operator be $L_h(u_h)$. Time integration is performed by an SSP RK scheme of order $k+1$,

$$u_h^{(0)} = u_h^n, \quad u_h^{(i)} = \sum_{m=0}^{i-1} \alpha_{im}\, u_h^{(m)} + \Delta t \sum_{m=0}^{i-1} \beta_{im}\, L_h(u_h^{(m)}),$$

setting $u_h^{n+1} = u_h^{(k)}$ at the completion of the step.

The coefficients $(\alpha_{im}, \beta_{im})$ are taken from the TVD RK (Cockburn–Shu) Butcher tableau.

CFL stability constraint is

$$\Delta t \leq C\, \min_j \frac{\Delta_j}{\max_i |\lambda_i|}, \quad C = \frac{1}{2k+1},$$

with characteristic speeds evaluated on δ-mapped states. The formal accuracy is $\mathcal{O}(h^{k+1})$.

6. Algorithmic Workflow

The full algorithm proceeds as:

1. Decompose the domain into cells $I_j$ and project $u$ onto $P^k(I_j)$.
2. At each discontinuous interface, compute δ-mapped states by supply-demand and entropy constraints.
3. Construct the DG spatial operator $L_h$ using the unified flux.
4. Apply the δ-mapped slope limiter to suppress oscillations near discontinuities.
5. Advance in time via SSP RK$(k+1)$ under the CFL constraint.

This process achieves high-order accuracy, preserves nontrivial steady states exactly, and suppresses non-physical oscillations arising from flux discontinuities.

7. Representative Applications

The RKDG + δ-mapping methodology has been successfully applied to:

• Nonlinear Elasticity in Layered Media: For systems

$$\varepsilon_t - v_x = 0, \quad (\rho v)_t - \sigma(\varepsilon, K)_x = 0, \quad \sigma(\varepsilon, K) = K \varepsilon + \beta K^2 \varepsilon^2,$$

with periodically varying $(\rho, K)$, incident pulses split into amplitude-dependent solitary waves; the scheme robustly captures wave speed and structure in agreement with analytic theory.

• Multi-Class Inhomogeneous Traffic Flow: For multiple vehicle classes $u_l = a(x) \rho_l$, with flux $f_l = b_l(x)\, u_l\, v(\rho)$, the δ-mapping equalizes supply/demand at interface jumps (lane number/speed limit), cleanly resolving up to $m+2$ nonlinear waves (shocks, rarefactions, contact discontinuities).

8. Generalization and Extensions

The δ-mapping concept can be integrated into any classical high-order shock-capturing scheme (WENO, TVD finite difference, central schemes), by replacing cell values $u_i$ by mapped values at discontinuous interfaces. This approach extends naturally to complex network problems, higher-order model couplings, and multidimensional domains with junction or interface conditions, providing a systematic route to high-order, entropy-consistent numerical simulation in the presence of spatially discontinuous coefficients.

---

For further technical details, including explicit construction of δ-mapped states, entropy constraints, and specific numerical benchmarks, see (Qiao et al., 2015).