Source-linked AI summary
Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method
Andreas Klöckner, Tim Warburton, Jan S. Hesthaven
TL;DR
Discontinuous solutions in nonlinear DG can produce Gibbs phenomena that the method must detect and mitigate. The paper develops a detector to guide artificial-viscosity shock capturing, reporting stabilization on nonlinear benchmark problems without excessive smoothing.
Problem
The paper seeks reliable detection of Gibbs phenomena that arise when nonlinear DG solutions develop discontinuities.
Source: §1 · p. 2Method
A shock detector supplies guidance for selective artificial-viscosity smoothing in DG, with rescaling allowing use in other shock-detection settings.
Source: §1 · p. 2Results
The method stabilized the shock-wave-interaction example and other important benchmarks without introducing excessive smoothness.
Source: §7 · p. 23Takeaways & Limitations
With appropriate rescaling, the detector may suit applications beyond guiding artificial viscosity that require shock detection.
Source: §1 · p. 2The explicit viscosity time-step constraint includes a potentially damaging diffusion factor N^4∥ν∥L∞/h^2 at high polynomial degrees.
Source: §2 · p. 5
Abstract
from arXiv · showhide
We present a novel, cell-local shock detector for use with discontinuous Galerkin (DG) methods. The output of this detector is a reliably scaled, element-wise smoothness estimate which is suited as a control input to a shock capture mechanism. Using an artificial viscosity in the latter role, we obtain a DG scheme for the numerical solution of nonlinear systems of conservation laws. Building on work by Persson and Peraire, we thoroughly justify the detector's design and analyze its performance on a number of benchmark problems. We further explain the scaling and smoothing steps necessary to turn the output of the detector into a local, artificial viscosity. We close by providing an extensive array of numerical tests of the detector in use.
1 Introduction
DG methods provide high-order approximations for hyperbolic conservation laws, but nonlinear solutions that develop discontinuities can produce Gibbs oscillations. This paper seeks reliable shock detection to guide selective artificial-viscosity smoothing within DG.
1 Introduction
- The paper targets reliable detection of Gibbs phenomena in discontinuous DG solutions, then uses artificial viscosity to mitigate them.Its detector is intended to guide selective smoothing and, with rescaling, may serve other shock-detection applications.Source: §1 · p. 2
- Unlike limiting approaches that modify solutions and can reach low orders of accuracy, artificial viscosity smooths discontinuities; a detector can guide that smoothing selectively.The authors position their detector as an improvement on existing detector-guided DG viscosity schemes.Source: §1 · p. 2
- Mesh adaptation can add resolution around discontinuities, but it too depends on a detector that reliably identifies where refinement is needed.The paper notes that artificial viscosity may also help when mathematical models contain discontinuities.Source: §1 · p. 2
1.1 The Discontinuous Galerkin Method
- DG approximates each hyperbolic system with local polynomials, using the same polynomial space for the numerical solution and test functions.The method rewrites the conservation law in a weak form on each element.Source: §1.1 · p. 3
- The lifting matrix combines facial contributions and lifts them into volume contributions.Its construction is illustrated in Figure 1.Source: §1.1 · pp. 3–4, Figure 1 · p. 4
- Applying the inverse mass matrix is element-wise and requires no global communication, enabling explicit time stepping and simplifying parallel implementation.The paper identifies this local operation as a distinguishing feature of DG relative to other finite-element methods.Source: §1.1 · p. 4
2 Basic Design Considerations
The design is motivated by GPU-oriented, matrix-free DG and the trade-offs of explicit integration. The section highlights viscosity-driven time-step constraints, the value of temporal adaptivity, and the limited influence of time-integrator choice in the reported setting.
2 Basic Design Considerations
- The work develops shock detection and capturing for DG as a step toward applying GPU-DG to nonlinear problems.The motivation comes from prior work mapping DG operators efficiently onto throughput-oriented GPUs.Source: §2 · p. 4
- Matrix-free methods suit GPUs’ memory and access constraints, but make implicit integration harder by requiring preconditioners and communication-intensive global reductions.The authors therefore examine explicit time stepping, whose small time-step restriction is its principal disadvantage.Source: §2 · pp. 4–5
- The viscosity-related diffusion time scale contains the factor N^4∥ν∥L∞/h^2, which can be damaging at high polynomial degrees.The section notes that high viscosity is expected to be localized in space and time, and that time adaptivity can address temporal localization.Source: §2 · p. 5
- In initial experiments, spatially generated oscillations dominated temporal ones at the encountered time steps, so changing the time integrator did not appreciably affect reported results.The discussion concerns oscillations originating in polynomial expansions of the spatial discretization.Source: §2 · p. 5
3 Applications and Equations
The paper applies artificial viscosity to advection, a first-order wave system, and Euler’s equations. Across the examples, the formulations use conservative diffusion, while the Euler sensor acts on density and a physical-viscosity form is insufficient to control its oscillations.
3.1 Advection Equation
- For advection, the paper adds viscosity in conservative form and discretizes diffusion with a first-order dual interior penalty method.The DG formulation uses a conventional upwind flux, with the gradient in strong form and divergence in weak form.Source: §3.1 · p. 6
3.2 Second-Order Wave Equation
- The second-order wave equation is rewritten as a first-order conservation-law system, letting the tests examine effects of two coupled characteristics.The sensor operates on scalar u, and the diffusive terms are discretized analogously to the advection case.Source: §3.2 · p. 6
3.3 Euler’s Equations of Gas Dynamics
- For Euler’s equations, one artificial viscosity applies to all components while the sensor operates on density ρ; a Navier-Stokes-like physical viscosity provides insufficient control of density oscillations.The system is intended for compressible, inviscid flow problems.Source: §3.3 · p. 7
4 A Smoothness-Estimating Detector
The detector builds on modal smoothness estimates but focuses on making the signal meaningfully scaled and robust to misleading modal coefficients. Its refinements include skyline pessimization and baseline modal decay, with a practical limitation when a discontinuity shifts toward an element boundary.
4 A Smoothness-Estimating Detector
- The detector addresses the scaling ambiguity of smoothness indicators by refining the Persson–Peraire approach, whose scaled value is intended to represent modal smoothness.The paper also notes a likely typographical error in Persson and Peraire’s mapping description: the comparable reference should scale with log(1/N^4).Source: §4 · pp. 7–8
4.1 Estimating Solution Smoothness
- The intended degree range is moderate to high polynomial order: low-order methods have flux-limiting alternatives, while degrees N ≳ 5 provide sufficient smoothness information but may face restrictive time steps.The authors hope the modal expansion is informative at lower degrees as well.Source: §4.1 · p. 8
- A least-squares fit to modal coefficient decay estimates an exponent, expected to be about 1 for a jump and 2 for a C0 kink under the Fourier-decay analogy.For the Heaviside interpolant, the fitted exponent is exactly 1.Source: §4.1 · pp. 9–10
- Skyline pessimization raises each coefficient to the largest higher-numbered coefficient, suppressing non-monotone decay that otherwise makes the kink appear spuriously smooth.The corrected kink estimate is 1.67, reasonably close to the expected value of 2.Source: §4.1 · pp. 11–12
- For the C1 but non-C2 truncated polynomial, skyline pessimization changes a raw estimate of about 13 to nearly the expected exponent of 3.The authors take the jump, kink, and truncated-polynomial tests as support for estimating smoothness from Legendre coefficients.Source: §4.1 · p. 12
- The highest representable mode in a length-10 expansion receives exponent 0 after skyline pessimization, while two smooth test functions receive values around 4.Source: §4.1 · p. 12
- A constant plus small white noise can receive a near-zero smoothness estimate because removing the constant mode loses scale; baseline modal decay scaled by the element norm is proposed to suppress that noise.Introducing baseline decay does not unduly modify the earlier experimental estimates.Source: §4.1 · pp. 12–13
- For a Heaviside jump shifted to an element edge, the final estimate is 0.57 rather than 1.05, and the paper reports no entirely practical remedy for this location sensitivity.The decay-fit indicator is described as more robust than a single-mode indicator, though slightly more expensive.Source: §4.1 · p. 14
5 From Smoothness to Viscosity
The section maps an estimated modal-decay exponent to an element-wise viscosity, then smooths that viscosity because elementwise constants are discontinuous. Its scaling balances diffusion and characteristic propagation, while its smoothing procedure avoids requiring higher-order continuity.
5.1 Scaling the Viscosity
- The activation targets ν(1) = ν0, ν(2) = ν0/2, and ν(3) = 0, corresponding roughly to discontinuous, continuous, and C1 solutions.The ramp reverses the earlier convention so that lower smoothness estimates produce greater viscosity.Source: §5.1 · p. 14
- Matching diffusion width to characteristic travel distance gives the viscosity scale; σ = 3h/N proved more realistic than σ = h/N.The authors say this reproduces Barter and Darmofal’s value while offering more insight into its meaning.Source: §5.1 · p. 15
- Using λmax is simple but sets ν0 higher than necessary for Euler contact discontinuities, and the scaling is imperfect across test problems.Most subsonic Euler problems benefit from an additional factor of 1/2, while Burgers problems work well without it in a picture-norm sense.Source: §5.1 · p. 15
5.2 Smoothing the Viscosity
- Because viscosity is constant per element, it is discontinuous; smoothing is preferred to avoid issues associated with discontinuous viscosity inside a conservative derivative.Earlier work reported issues from localized discontinuous viscosities and a marked decrease in H1 error for smooth viscosities.Source: §5.2 · pp. 15–16
- Diffusion-based post-processing can make viscosity smooth in space and time, but sufficient smoothing requires large diffusivity and yields a stiff ODE system.Experiments found no apparent advantage to requiring ν ∈ Ck for k > 0.Source: §5.2 · p. 16
- The chosen smoother takes the maximum adjacent-element viscosity at each vertex, propagates it back to adjoining elements, and extends it with a linear P1 interpolant.The authors report that this approach is cheap and reasonably straightforward to implement, including on GPUs.Source: §5.2 · p. 16
6 Experience with and Evaluation of the Scheme
Across advection, wave, and Euler tests, artificial viscosity suppresses visible advection overshoots and stabilizes classical Euler cases, but interacting regions show reduced convergence. The Euler Sod study finds about first-order convergence in h and less than first-order convergence in polynomial degree.
6.1 Advection: Basic Functionality, Interaction with Time Discretization
- With viscosity, the advection jump is resolved within eight points—less than one element—without visible overshoots, and its steep profile remains stable over repeated round-trips.The comparison is shown against the unmodified DG solution, which exhibits a strong Gibbs-type overshoot.Source: §6.1 · p. 17, Figure 6 · p. 17
- Viscosity spikes recur as discontinuities cross element boundaries; spatial sensitivity delays a profile that no longer triggers the detector until about t⪆12.The final profile is expected to be twice smoother than necessary absent the detector’s oversensitivity.Source: §6.1 · pp. 17–18, Figure 7 · p. 18
- Explicit artificial-viscosity shock capture requires adaptive time stepping to remain competitive; chosen steps show an initial settling period and later reductions.Some reductions coincide with discontinuity transitions, though the pattern is less clear than for viscosity activations.Source: §6.1 · p. 18, Figure 7 · p. 18
6.2 Waves: Shock Spreading and Spurious Coupling
- Without viscosity, full convergence holds away from discontinuities, while a slowly growing non-convergence zone remains near them.The study notes that this behavior matches predictions for the inviscid DG scheme.Source: §6.2 · p. 19, Figure 8 · p. 19
- With artificial viscosity, wave-equation convergence remains roughly unchanged outside interacting regions but does not recover full order within them, consistent with cross-coupling from the added viscosity.The reported loss of accuracy also affects the smooth rarefaction-wave profile in the Euler case.Source: §6.2 · p. 19, Figure 8 · p. 19, §6.3 · p. 20
6.3 Euler’s Equations
- Sod solutions show element-scale wrinkles on constant states inside the shock fan; their magnitude appears to grow with resolution, while their cause and link to reduced convergence remain unknown.The oscillations do not appear to contribute significantly to the error at the examined resolutions.Source: §6.3 · pp. 21–22, Figure 10 · p. 21
- Sod tests support about first-order convergence in h but less than first order in N, making added elemental resolution more effective than higher polynomial degree.The study reports that increasing N does not improve convergence, while computational complexity grows superlinearly in N.Source: §6.3 · p. 21, Table 1 · p. 22
- Euler tests produce a correct Lax solution despite a more-smeared contact, and a stable, meaningful Shu–Osher solution without excessive smearing, though high-order convergence is not established there.For sufficiently resolved smooth isentropic vortices, the detector does not activate.Source: §6.3 · p. 22, Figure 11 · p. 23
7 Conclusions and Future Work
The method focuses on reliably scaling a local shock detector for explicit, GPU-suited DG calculations. Nonlinear benchmark results are encouraging, while simple-problem behavior and multidimensional generalization remain areas for further investigation.
7 Conclusions and Future Work
- Reliable scaling, explicit local GPU-suited calculation, and new detector and viscosity-smoothing concepts distinguish the method.These include skyline pessimization, baseline decay, and P 1 viscosity smoothing.Source: §7 · p. 23
- Simple linear-wave and transport problems revealed behavior the authors say should be investigated in other shock-capturing methods.Source: §7 · p. 23
- The method stabilizes shock-wave interaction and other important nonlinear benchmarks without introducing excessive smoothness.The authors characterize these results as encouraging.Source: §7 · p. 23
- Further testing on established shock-capturing benchmarks and investigation of multidimensional generalization remain future work.Initial multidimensional tests showed promising results, which the authors planned to report in a future article.Source: §7 · p. 23