PaperlayerView exact revision on arXiv

Source-linked AI summary

Stiff source terms inside ADER-DG method

I. S. Popov, A. P. Popova

arXiv:2610.30852v1physics.comp-phcs.NAmath-phmath.MPmath.NAphysics.app-phphysics.flu-dyn

TL;DR

ADER-DG with an LST-DG predictor faces stability questions for hyperbolic equations with stiff source terms, where established Courant-number estimates leave uncertainty. The paper develops rigorous approximation and stability analyses and evaluates the method on linear and reacting-flow problems. It finds that the maximum Courant number rises with stiffness and reports accurate solutions in the reacting-flow tests.

  • Problem

    The paper addresses stability and approximation of ADER-DG with an LST-DG predictor for hyperbolic equations with stiff source terms, a setting not covered by the cited source-free stability analysis.

    Source: Introduction · p. 2

  • Method

    The paper constructs rigorous approximation and von Neumann stability analyses for ADER-DG with an LST-DG predictor, deriving Courant-number conditions from polynomial equations.

    Source: Introduction · p. 1, §3.2 · p. 22, Conclusion · p. 87

  • Results

    For |z| ≫ 1, CFLmax(N, z) increases practically linearly with |z|, while for |z| ≪ 1 it reaches the source-free CFLmax(N, 0) plateau.

    Source: Conclusion · p. 88

  • Takeaways & Limitations

    In the tested reacting-flow problems, the method reproduces the detonation front, chemical spike, contact discontinuity, and small rarefaction wave accurately on fairly coarse grids.

    Source: Conclusion · p. 89

    The stability theory is based on a linear advection-reaction equation, while quasilinear systems can develop discontinuities even from smooth initial conditions.

    Source: §4.2 · p. 72, Conclusion · p. 87

Abstract

from arXiv · show

ADER-DG numerical method with the LST-DG predictor, when used to solve hyperbolic equations with stiff source terms, exhibits properties characteristic of IMEX computational schemes. With increasing stiffness, characterized by the parameter $z = \nu{\Delta t}^{n} < 0$, the maximum of the Courant number $\mathrm{CFL}_{\rm max}$ of the ADER-DG numerical method increases. In the region of strong stiffness $|z| \gg 1$, a linear increase in the Courant number $\mathrm{CFL}_{\rm max}(z) \sim |z|$ occurs. In the region of small source terms $|z| \ll 1$, the results $\mathrm{CFL}_{\rm max}(z)$ transform into known results for the problem $\mathrm{CFL}_{\rm max}(0)$ without source terms. In the transition region $|z| \sim 1$, a range of insignificant decrease in the Courant number $\mathrm{CFL}_{\rm max}(z)$ is discovered. In the linear case, instability is violated precisely when one of the matrix eigenvalues reaches $\lambda = -1$, regardless of the phase $\theta$ of the numerical solution. A rigorous stability theory of the ADER-DG numerical method for stiff problems for arbitrary degree $N$ of basis polynomials is constructed. The stability condition is simplified to solving a polynomial equation for the Courant number $\mathrm{CFL}$. Highly accurate dependencies of the Courant number $\mathrm{CFL}_{\rm max}(N, z)$ of the ADER-DG method on the parameter $z$ are obtained. A rigorous direct proof of the approximation of the ADER-DG method is presented. Approximation orders $p = N+1$ for arbitrary degrees $N$ are rigorously derived. A set of numerical experiments is carried out to apply the ADER-DG method to solving both a linear advection-reaction equation and the system of Euler equations for a two-component reacting medium. The results obtained in these calculations confirm the theoretical results well.

Introduction

The introduction places ADER-DG with an LST-DG predictor among high-order methods used across many PDE applications, including stiff reacting flows. It motivates a rigorous study of stability and approximation for stiff source terms because existing Courant-number estimates and computations leave important uncertainty.

Introduction

  • ADER-DG methods have been developed from the ADER paradigm and are used in applications ranging from hydrodynamics to reacting flows.The LST-DG predictor uses a discrete space-time solution and was reported to improve accuracy over the original finite-volume ADER version, particularly for stiff problems.Source: Introduction · pp. 1–2
  • Strict time-step limits and small maximum Courant numbers, especially at high polynomial degrees, constrain ADER-DG calculations.The method's small CFL limits are relevant to subcell limiting, where finite-volume limiters can use less restrictive time steps.Source: Introduction · p. 2
  • Existing high-degree CFL estimates differ, and the widely used Cockburn–Shu estimate can overestimate ADER-DG limits, sometimes by more than twofold.Prior estimates relied on spectral-radius calculations over a continuous range of phases, with eigenvalue robustness, matrix operations, and finite precision affecting the computation.Source: Introduction · p. 2
  • Prior rigorous stability analysis for the source-free case motivates investigating whether stiff source terms raise the ADER-DG maximum Courant number.The introduction contrasts that analysis with detonation calculations using an effective Courant number of 0.90 even at degrees N ≥ 5, suggesting a higher CFL limit under stiffness.Source: Introduction · p. 3
  • The paper develops rigorous approximation and stability analyses for stiff-source problems, then tests the method on linear advection-reaction and reacting Euler equations.Its roadmap includes approximation orders, stability bounds, and numerical experiments for a two-component reacting medium.Source: Introduction · p. 3

1 General description of the numerical method

The section formulates ADER-DG with an LST-DG predictor for a linear advection-reaction equation. It specifies the spatial and space-time discretizations used to derive the predictor and the single-step update.

1 General description of the numerical method

  • The target equation combines advection at velocity a with a decay term νu, where ν ≤ 0, and serves as a PDE analogue of the Dahlquist test equation.The analysis assumes sufficiently smooth initial data; otherwise the method's approximation order decreases.Source: §1 · pp. 3–4
  • The spatial and temporal domains are discretized into coordinate-time elements and mapped to a reference element for the weak-form scheme.The mesh uses constant spatial steps for simplicity; the time discretization uses steps Δt^n.Source: §1 · p. 4
  • The method uses a nodal Lagrange basis at shifted Legendre roots, with Gauss–Legendre quadrature and pointwise expansion coefficients.The quadrature is exact for P2N+1, and the coefficient representation is exact for functions in PN.Source: §1 · p. 4
  • The resulting ADER-DG scheme is a completely single-step finite-element update using the local predictor and a classical Riemann solver for interface fluxes.The Riemann solver has an exact closed form for the advection equation.Source: §1 · p. 5
  • The LST-DG predictor computes a local discrete space-time solution q on each coordinate-time element to supply local and flux terms.The predictor represents q using a tensor-product basis and imposes the initial condition in weak form, allowing a discontinuity at the initial time.Source: §1 · p. 5
  • The predictor's algebraic system uses matrices assembled from the basis, including the mass and stiffness matrices, and can be expressed with Kronecker products.For the chosen basis, the mass matrix is diagonal through L2 orthogonality.Source: §1 · p. 6

2 Approximation analysis

Section 2 rigorously analyzes approximation of the ADER-DG method with the LST-DG predictor, examining its local discrete space-time solution and corrected time-step solution. It estimates local and interface contributions to the update, establishes approximation order p = N + 1, and gives order 2N + 1 for the exponential-integrator flux approximation.

2.1 Discrete space-time solution

  • The LST-DG predictor system is well-defined because its matrix κ is invertible, as shown by proving that κg = 0 implies g = 0.The proof uses the linear independence of the basis functions and a choice of test functions to establish the trivial kernel.Source: §2.1 · p. 7
  • The local discrete space-time solution has approximation order p_L = N + 1 within the element-local region.The estimate applies locally within element Ξ_i^n.Source: §2.1 · p. 12

2.2 Corrected solution at time steps

  • The corrected time-step analysis rewrites the one-step ADER-DG update and uses the ODE stability function and integration by parts to isolate terms that compensate for fluxes.The ODE stability function is introduced for the time-DG predictor, and integration by parts is used to isolate the flux-compensation terms.Source: §2.2 · pp. 12–14
  • The exponential-integrator approximation of the numerical total flux has order 2N + 1, above the DG solution representation's order N + 1.The flux estimate is used to analyze the interface contribution to the update.Source: §2.2 · p. 18
  • For sufficiently smooth solutions, the local error construction includes terms estimated as O((Δt^n)^(N+1)) and O((Δt^n)^(2N+1)).These local estimates are combined with interface estimates to obtain the final approximation order.Source: §2.2 · p. 20
  • The ADER-DG method with LST-DG predictor has approximation order p = N + 1.The final error estimate is derived by combining local and interface contributions.Source: §2.2 · p. 20

3 Stability analysis

This section develops qualitative and rigorous quantitative stability analyses for the ADER-DG method, characterizing its stability boundary through spectral analysis and a phase-independent polynomial-root procedure. It examines how the boundary varies with source-term stiffness and polynomial degree, including its occurrence when an eigenvalue reaches −1.

3.1 Qualitative considerations

  • As stiffness increases, the ADER-DG/LST-DG stability region expands, consistent with IMEX-like behavior; the high-degree asymptotic CFL estimate is useful for preliminary cost estimates but overestimates exact values by 21–47%.The estimate is asymptotically relevant for large polynomial degrees N, not a reliable quantitative prediction of CFLmax.Source: §3.1 · pp. 21–22

3.2 Statement of the problem of quantitative considerations

  • The Von Neumann analysis represents one-step evolution by a matrix R(CFL, z, θ), and the method rigorously reproduces spatially constant solutions’ reaction decay independently of CFL.For the corresponding ODE, the decay is given by a Padé approximation with approximation order 2N + 1.Source: §3.2 · pp. 22, 27

3.3 Empirical considerations

  • Because arbitrary-degree stability boundaries are not analytically tractable, the analysis combines rigorous algebraic equations for CFLmax with empirically motivated assumptions about where eigenvalue boundaries occur.Stability requires the spectral radius of R(CFL, z, θ) to be at most one for every phase θ.Source: §3.3 · p. 27, §3.2 · p. 27
  • Figure 1 plots eigenvalue magnitudes across phase for selected CFL values, source stiffnesses, and polynomial degrees, showing stable-boundary and nearby stable and unstable cases.The figure samples 1,000 phase values for z = −0.1, −1.0, −10.0, and −100.0, with N = 1, 2, and 3.Source: Figure 1 · p. 28
  • As |z| increases for z < 0, all but one eigenvalue decrease in magnitude; the largest-magnitude eigenvalue occurs at θ = 0 for odd N and θ = π for even N.For N ≥ 4, the compressed spectrum does not decrease entirely monotonically with increasing |z|.Source: §3.3 · pp. 30–31
  • Advection-reaction stability is lost strictly when an eigenvalue reaches -1, unlike advection, where tiny 10^-5–10^-9 excesses near +1 occur at CFL values 1.5–4.0 times below CFLmax.The exception is limited to |z| ≤ 10^-5–10^-6 and is described as having no significant interest.Source: §3.3 · p. 33
  • The eigenvalue set changes from one closed curve to two as |z| increases, with the transition at z = -10.0 for N = 1–3 and z = -1.0 for N = 4–12.One curve contains the eigenvalue near -1; the other contains the remaining eigenvalues closer to 0.Source: §3.3 · pp. 33–34

3.4 Rigorous quantitative considerations

  • The quantitative analysis transforms the characteristic equation to remove its phase factor and derives a phase-independent relation for the Courant stability boundary, where an eigenvalue reaches -1.The derivation uses the characteristic equation of the evolution operator and transformations that eliminate the phase θ dependence.Source: §3.4 · pp. 34, 37
  • CFLmax is the infimum of Courant numbers for which the phase-independent boundary equation has a solution for some complex λ with |λ| ≥ 1.The criterion combines the cases a > 0 and a < 0.Source: §3.4 · p. 40
  • For a > 0 and a < 0, the boundary equations yield the same F values after accounting for the sign change CFL ↔ −|CFL|.Source: §3.4 · p. 40
  • Across z = −0.1, −1.0, −10.0, and −100.0 and N = 1,…,12, the stability boundary occurs at ϕ = π (λ = −1), with no false roots near ϕ = 0 or 2π.The reported roots confirm that λ = −1 marks the stability limit.Source: §3.4 · p. 40
  • Setting λ = −1 reduces the boundary condition to F(π; CFL, z) = 0, whose solution determines CFLmax; symbolic A and B expressions are tabulated for N = 1,…,6.The functions A(CFL, z) and B(CFL, z) are defined from H+(1; CFL, z, −1) and H+(0; CFL, z, −1), respectively.Source: §3.4 · p. 40, Table 1 · p. 48, Table 2 · p. 49
  • The expressions were computed symbolically for N=1–12, with those for N=1–5 presented; 5000-digit precision was used to address cancellation when finding roots.For N=12, the highest-CFL coefficient is a polynomial in z of degree 1170.Source: §3.4 · pp. 52, 54
  • The stability boundary CFLmax is determined by the minimum positive real root of Φ+(CFL; N, z)=0 or Φ−(CFL; N, z)=0, polynomials of degrees N+1 and N, respectively.The expressions are obtained by clearing the common denominator d(z).Source: §3.4 · p. 52
  • For N=1–12, CFLmax approaches CFLmax(N, 0) as |z|≪1, dips slightly near |z|∼1, then grows linearly for |z|≫1; the dip is at most 6% for N=12.N=1 has no intermediate decrease; |zmin| increases from 0.202 at N=2 to 0.796 at N=12.Source: Figure 15 · p. 55, §3.4 · pp. 55–56, Table 5 · p. 56
  • For |z|≫1, CFLmax(N, z)≈r(N)|z|, and r(N) decreases with N, primarily alongside CFLmax(N, 0)∝1/(N+1)^2.The qualitative estimate nearly matches for N=1; its differences stop systematically increasing from N=8 and remain approximately 26% or less.Source: §3.4 · p. 56, Table 5 · p. 56

4 Applications

The applications evaluate the ADER-DG method with the LST-DG predictor on linear advection–reaction and stiff reacting-flow problems, including ZND detonation waves. They examine its behavior across source stiffness and polynomial degrees, including high-Courant solutions and the resolution of reaction-dependent structures.

4.1 Linear advection-reaction equation

  • The compensating source s(t, x) = −νuinit(x − at) makes the exact solution a pure advection of the initial profile, and its added term does not change the homogeneous stability condition.The inhomogeneous term enters the scheme as an additional right-hand-side contribution proportional to Δt.Source: §4.1 · pp. 57–58
  • For strong stiffness, the tests demonstrate stable computations when CFLreact ≫ CFL ≫ CFLmax(N, 0), with the gap reaching at least 10^3–10^5.Here z = |ν|Δt is interpreted as an effective reaction Courant number; for high polynomial degrees, good agreement with the exact solution is reported even on 10–20-cell meshes.Source: §4.1 · pp. 60–61
  • Empirical convergence orders generally agree with p = N + 1, with modest reductions at N > 9 for CFL ≥ 1.0.The reductions are attributed to accumulated computational error over the long integration time; the method still obtains high convergence orders for strongly stiff problems at CFL values far above the source-free stability limit.Source: Table 7 · p. 63, §4.1 · pp. 63–66
  • At CFL = 100.0 and 1000.0, the solution can develop complex artifacts, and very large time steps produce significant dispersion-law deviations.The authors describe this behavior as expected and note that convergence may persist pointwise but becomes significantly non-uniform.Source: §4.1 · pp. 66–68
  • At CFL = 10.0, the stability boundary is almost 100 times the source-free limit for N = 3 and almost 1000 times that limit for N = 12.The reported boundary values increase with stiffness; the tests also show good qualitative agreement with the exact solution for CFL up to 10.0, including meshes of 10–20 cells at N ≥ 8–9.Source: §4.1 · pp. 60, 68
  • The method retains p = N + 1 convergence and is stable at CFL = 100.0 and 1000.0, though significant, clearly visible numerical artifacts arise.These results are reported for the ADER-DG method with the LST-DG predictor.Source: §4.1 · p. 69

4.2 ZND-detonation waves

  • In reacting systems, source Jacobians generally mix characteristic modes; strong negative source eigenvalues are expected to expand stability, but not with the linear problem’s exact CFL limits.The source term is not generally diagonal in the flux-characteristic basis, and a linearized system can have additional source eigenvalues.Source: §4.2 · p. 70
  • For quasilinear systems, solution discontinuities can arise from smooth initial data, so the linear ADER-DG method requires a posteriori correction to handle general cases.The correction uses PAD and NAD admissibility checks and a robust finite-volume limiter in troubled cells.Source: §4.2 · pp. 72–73
  • The nonlinear application is a two-component A → B detonation test using slow and fast reaction kinetics, and it omits physical-process splitting.The fast case uses τ0 = 4 · 10^-3 and the slow case τ0 = 10^-1; the study disables global or local time-step corrections for stiffness. It uses CJ initial conditions and a 5000-cell ADER-WENO2 reference solution.Source: §4.2 · pp. 73–75
  • For slow reactions, the extended burnout zone requires resolving kinetics as well as hydrodynamics; N=2–3 reproduce reference features at C=0.90–1.15.N=3 achieves this on 200 cells; N=4–6 produce high-quality solutions on 100–50 cells, with maximum effective C=1.10 for N=4 and N=6.Source: §4.2 · pp. 75, 77, 78, 79
  • Fast reactions create a very narrow burnout region, so inaccurate space-time resolution can trigger a spurious weak detonation ahead of the shock and prevent the stationary wave from forming.Source: §4.2 · pp. 83–84
  • For fast reactions, N≥7 improves solution quality, and N≥9 resolves all reference features on 40 cells with the detonation front smoothed by no more than one cell.For N≥7, the maximum effective C for a correct solution is 1.00.Source: §4.2 · p. 86
  • In both weak- and strong-stiffness cases, correct solutions occur above the source-free CFLmax(N, 0), with the boundary at least CFLCS(N).Source: §4.2 · pp. 76, 85, 87
  • CFLmax reached ≈1.35 CFLNSmax(1) for N=1 in both stiffness cases, ≈1.35 CFLNSmax(2) for weak-stiff N=2, and ≈1.46 CFLNSmax(2) for strong-stiff N=2.For weak-stiff N=3, CFLmax reached ≈1.58 CFLNSmax(3).Source: §4.2 · p. 87
  • At N = 5, the trend reverses: CFLmax is ≈2.00 CFLNS for weak stiffness and ≈1.90 CFLNS for strong stiffness.For degrees up to N = 4, the strong-stiffness boundary was not lower than the weak-stiffness boundary.Source: §4.2 · p. 87
  • For N ≥ 12, CFLmax with stiff sources is 3.0 or more times greater than CFLNS, the boundary Courant number without stiff source terms.The increase is less dramatic than for the linear advection-reaction equation, but the method is reported to suit strongly stiff ZND problems.Source: §4.2 · p. 87
  • For N ≥ 10, the stiff-source boundary CFLmax is well described by the Cockburn-Shu formula (1).At N = 11 and N = 12, CFLmax is ≈3.19 CFLNS in both weak- and strong-stiffness cases.Source: §4.2 · p. 87
  • The results explain the use of effective Courant number C = 0.90 in ZND detonation simulations with this method.The passage attributes the more limited increase in quasilinear systems to their different source-term stiffness structure.Source: §4.2 · p. 87

Conclusion

The work develops stability and approximation theory for ADER-DG with an LST-DG predictor on stiff-source hyperbolic problems. It finds that stiffness raises the allowable Courant number, while numerical experiments support the theory for linear and reacting-flow problems.

Conclusion

  • For the linear advection-reaction equation, instability occurs when an evolution-matrix eigenvalue reaches λ = -1, and high-precision calculations found no earlier small unstable eigenvalues of practical concern.The phase-eliminated stability function has roots at phase φ = π.Source: Conclusion · p. 88
  • For arbitrary polynomial degree N, the stability boundary can be computed by solving two degree-(N + 1) polynomial equations and selecting the smallest positive real root.The coefficients are polynomials in z, avoiding a detailed phase search and potentially unstable eigenvalue calculations.Source: Conclusion · p. 88
  • The ADER-DG method with the LST-DG predictor has rigorously established approximation order p = N + 1 for arbitrary N, and a constant-initial-condition decay property of order p = 2N + 1.The decay property is described as similar in essence to the constant-preserving property for linear advection.Source: Conclusion · p. 88
  • The stability boundary CFLmax grows almost linearly with |z| for strong stiffness, while remaining much larger than the source-free boundary CFLmax(N, 0) but much smaller than the reaction Courant scale.For |z| ≫ 1, the reported relation is CFLreact ≫ CFLmax(N, z) ≫ CFLmax(N, 0).Source: Conclusion · pp. 87–88
  • For |z| ≪ 1, CFLmax returns to the source-free value, whereas a narrow dip near |z| ∼ 1 stays within about 5–6% even at N = 12.In the strong-stiffness regime, the growth slope decreases with N in association with the decreasing source-free CFLmax(N, 0).Source: Conclusion · p. 88
  • In linear advection-reaction tests, empirical convergence orders agreed with p = N + 1 under strong conditions, while extremely large CFL values produced dispersion artifacts despite preserved pointwise convergence.The tests used ν = −10^2 and −10^3; at CFL = 100 and 1000, artifacts arose when a time step traversed the domain.Source: Conclusion · p. 88
  • In two-component reacting Euler tests, including fast-reaction detonation cases, accurate solutions were obtained on coarse grids even at CFL > 3 CFLmax(N, 0), supporting applicability to quasilinear stiff-source systems.The solutions reproduced the detonation front, chemical spike, contact discontinuity, and a small rarefaction wave.Source: Conclusion · pp. 88–89

Related papers

Loading 2610.30852v1…