PaperlayerView exact revision on arXiv

Source-linked AI summary

Second-order invariant domain preserving approximation of the Euler equations using convex limiting

Jean-Luc Guermond, Murtazo Nazarov, Bojan Popov, Ignacio Tomas

arXiv:1710.00417v1math.NAcs.NA

TL;DR

The paper addresses how to approximate the Euler equations while preserving density, internal-energy, and specific-entropy constraints. It combines a first-order invariant-domain-preserving scheme with high-order approximation and convex limiting; the intermediate states satisfy the stated physical bounds, while reported test cases include first-order convergence rates.

  • Problem

    Standard FCT is not well suited to enforcing the quasiconcave constraints of positive internal energy and a minimum principle on specific entropy.

    Source: §1 · p. 2

  • Method

    An algebraic, convexity-based limiting process uses invariant-domain-preserving states to bound and limit a high-order approximation for general concave Euler constraints.

    Source: §1 · p. 2, §4 · p. 11, §4.2 · p. 12

  • Results

    Intermediate states have positive density and internal energy, satisfy the local specific-entropy minimum, and are discretely entropy satisfying; the rarefaction test reports O(h) when entropy limitation is unrelaxed, and the Leblanc shocktube reports first-order convergence for δ1(t).

    Source: §3.2 · p. 8, §5.3 · p. 21, §5.4 · p. 22

  • Takeaways & Limitations

    The algebraic formulation is described as applicable beyond continuous finite elements, including finite-volume and discontinuous Galerkin discretizations.

    Source: §1 · p. 2, §3.1 · p. 7

    Total-energy limitation was not found useful and was omitted from the reported numerical tests.

    Source: §4.4 · p. 15

Abstract

from arXiv · show

A new second-order method for approximating the compressible Euler equations is introduced. The method preserves all the known invariant domains of the Euler system: positivity of the density, positivity of the internal energy and the local minimum principle on the specific entropy. The technique combines a first-order, invariant domain preserving, Guaranteed Maximum Speed method using a Graph Viscosity (GMS-GV1) with an invariant domain violating, but entropy consistent, high-order method. Invariant domain preserving auxiliary states, naturally produced by the GMS-GV1 method, are used to define local bounds for the high-order method which is then made invariant domain preserving via a convex limiting process. Numerical tests confirm the second-order accuracy of the new GMS-GV2 method in the maximum norm, where 2 stands for second-order. The proposed convex limiting is generic and can be applied to other approximation techniques and other hyperbolic systems.

1. Introduction.

The introduction frames second-order Euler approximation with non-affine invariant constraints as a limitation of standard FCT, then motivates an algebraic, convexity-based limiting approach.

1. Introduction.

  • FCT methods can enforce pointwise bounds such as density positivity, but are not well suited to quasiconcave constraints on internal energy and specific entropy.The paper identifies convex combinations and concavity as a basis for extending limitation to these constraints.Source: §1 · pp. 1–2
  • The proposed approach treats limitation algebraically rather than limiting slopes, with the aim of enforcing general concave Euler constraints.The method is presented for continuous finite elements and is described as applicable to finite volumes and discontinuous Galerkin methods.Source: §1 · p. 2
  • The paper proceeds from Euler-equation and finite-element preliminaries to approximation methods, novel local bounds and limiting results, and numerical tests.The local bounds and Lemmas 4.3–4.4 are identified as the main results of Section 4.Source: §1 · p. 2, §2 · p. 2

2. Preliminaries.

The preliminaries define the Euler variables and entropy properties used by the method, then establish the continuous finite-element setting and its algebraic notation.

2.1. The Euler equations.

  • The Euler equations are written in conservative variables: density, momentum, and total energy, with velocity from momentum and specific internal energy derived from total energy.Pressure is determined by an equation of state; the exposition assumes periodic boundaries or compactly supported initial data before boundary influence arrives.Source: §2.1 · pp. 2–3

2.2. Intrinsic properties.

  • The set A_smin requires positive density and internal energy and a specific-entropy floor; it is convex and invariant, and generalized entropy inequalities describe entropy solutions.The paper's objective is an explicit second-order continuous finite-element method consistent with the entropy inequality and preserving A_smin.Source: §2.2 · p. 4

2.3. Finite element setting.

  • The finite-element setting uses continuous Lagrange spaces on mapped, shape-regular meshes, with nodal degrees of freedom and partition-of-unity shape functions.The text allows combinations of reference elements, including triangles and parallelograms in two dimensions.Source: §2.3 · pp. 4–5
  • Positive lumped masses and pairwise vectors c_ij provide the algebraic notation for conservation, while SSP Runge–Kutta time stepping is used in the numerical tests.The tests use SSP RK(3,3); the discussion develops the spatial techniques for forward Euler stepping.Source: §2.3 · p. 5, §3 · p. 6

3. Time and space approximation.

The approximation pairs a guaranteed-speed, first-order invariant-domain-preserving method with high-order updates, then constructs intermediate states and viscosity choices to support accuracy, entropy consistency, and limiting.

3.1. Low-order approximation (GMS-GV1).

  • GMS-GV1 uses a guaranteed upper bound on the Riemann-problem wave speed to define graph viscosity for the low-order update.The low-order method is described as the Guaranteed Maximum Speed method with first-order Graph Viscosity.Source: §3.1 · pp. 6–7

3.2. The intermediate limiting states.

  • Under the appropriate CFL condition, intermediate states are Riemann-solution averages, and the low-order update has positive density and internal energy, a local specific-entropy minimum principle, and a discrete entropy inequality.The entropy inequality is stated for generalized entropy pairs.Source: §3.2 · pp. 7–8

3.3. Smoothness-based approximation.

  • The smoothness-based viscosity uses a local smoothness indicator and is inherently only second-order accurate in space.The indicator is linearity-preserving for suitable positive parameters and is of order h in smooth regions without local extrema.Source: §3.3 · p. 8

3.4. Entropy viscosity commutator.

  • The entropy-viscosity method is formally high-order for any polynomial degree and aims to make the provisional update entropy consistent and close to invariant-domain preserving.The authors report that it is more robust than the smoothness-based viscosity, albeit slightly more oscillatory.Source: §3.4 · pp. 9–10
  • The entropy residual viscosity is at least one mesh-size order smaller than the low-order graph viscosity under the stated nonzero-entropy-derivative assumption.For polytropic gases, the authors use p^(1/γ) to construct the entropy residual in their reported tests because they found it minimizes dissipation across contact discontinuities.Source: §3.4 · pp. 10–11

4. Quasiconcavity-based limitation.

The high-order update may violate physical bounds, motivating a convexity-based limiting technique that enforces constraints without several common auxiliary strategies.

4. Quasiconcavity-based limitation.

  • The proposed convexity-based limiting technique avoids linearization, worst-case estimates, a posteriori fixes, and auxiliary discontinuous spaces.The authors trace its roots to Khobalatte and Perthame, Perthame and Qiu, and Perthame and Shu.Source: §4 · p. 11
  • The high-order update is not guaranteed to be oscillation-free or to preserve physical bounds, so it requires limitation.For the Euler equations, the text identifies the specific-entropy minimum principle as the exception to the absence of known maximum or minimum principles.Source: §4.1 · p. 11
  • Some existing limiters enforce properties not necessarily satisfied by the Euler equations, or only satisfied by a chosen first-order method.The passage notes that many such approaches use a Lax-Friedrichs-like first-order scheme.Source: §4.1 · p. 11

4.1. Bounds and quasiconcavity.

The bounds are based on properties naturally satisfied by the low-order solution, with quasiconcavity providing a unified framework for the constraints.

4.1. Bounds and quasiconcavity.

  • The strategy enforces bounds naturally satisfied by the low-order solution, alongside the local minimum principle on specific entropy.These bounds are used when modifying the graph viscosity for the high-order update.Source: §4.1 · p. 11
  • When the specific-entropy equation is unavailable, limiting the concave internal energy ρe offers an alternative discrete minimum principle.The passage states that ρe is concave in the conservative variables irrespective of the equation of state.Source: §4.1 · p. 12
  • Quasiconcavity unifies density, internal-energy, and specific-entropy constraints through convex upper level sets.The paper gives density, positive internal energy, and a specific-entropy lower bound as the three settings.Source: §4.1 · p. 12
  • The equivalent concave formulations use ρ, ρe, and ρ(Φ(U) − r), which the authors say are easier to constrain than quasiconcave functions.The passage notes that Newton's method on a concave function is guaranteed to converge independently of the initial guess, unlike the quasiconcave case.Source: §4.1 · p. 12

4.2. An abstract limiting scheme.

The scheme extends FCT-style limiting to general quasiconcave constraints while retaining conservative updates. Symmetric limiting parameters and convexity of constraint upper-level sets underpin the construction, which also extends to other hyperbolic systems.

4.2. An abstract limiting scheme.

  • The method targets general convex constraints, including the specific-entropy minimum principle, that Zalesak-style FCT cannot easily handle without losing second-order accuracy.The limiting process is convexity-based.Source: §4.2 · p. 12
  • The provisional high-order solution is expressed as a correction to the low-order solution, and symmetric limiting parameters scale that correction while preserving conservation.The passage states that the skew-symmetry of A^n and symmetry of the limiter imply conservation.Source: §4.2 · pp. 12–13
  • For a quasiconcave constraint Ψ, admissible pairwise limited states yield an admissible combined update because the constraint’s upper-level set is convex.The limiter parameters are symmetric, ℓ_ij = ℓ_ji.Source: §4.2 · pp. 13–14
  • The limiting lemmas also apply to arbitrary hyperbolic systems when their invariant domains are described by quasiconcave constraints.The section presents these lemmas as the basis for enforcing the constraints.Source: §4.2 · p. 14, §4.3 · p. 14

4.3. Application to the Euler equations.

The section applies the limiting lemmas to the Euler constraints, imposing bounds on density, internal energy, and specific entropy. It also notes that extra total-energy limiting is not used in the reported tests and that internal-energy limiting can over-dissipate contact discontinuities.

4.3. Application to the Euler equations.

  • The Euler application limits density, internal energy, and specific entropy using the quasiconcave constraints described for each quantity.The section identifies these as the constraints imposed in the algorithm.Source: §4.3 · p. 14
  • The numerical simulations use k_max = 1 for the iterative limiting algorithm.The passage specifies this setting for the reported simulations.Source: §4.3 · p. 14
  • Additional limiting based on other quasiconcave constraints, such as bounds on total energy, did not improve performance; internal-energy limiting tends to over-dissipate contact discontinuities.The reported tests limit density, internal energy, and specific entropy.Source: §4.3 · p. 15
  • The section states that the described method is independent of the equation of state.Source: §4.3 · p. 15

4.4. Limitation on the density.

For a covolume equation of state, the stated invariant domain includes positivity of density and internal energy, a specific-entropy lower bound, and bρ < 1. The method can enforce the additional affine constraint 1 − bρ > 0; total-energy limiting is described as unhelpful and omitted from the tests.

4.4. Limitation on the density.

  • For the covolume equation of state p(1−bρ) = (γ−1)ρe, the stated invariant domain requires ρ > 0, e > 0, s ≥ s_min, and bρ < 1.Source: §4.4 · p. 15
  • The method can enforce the additional affine constraint 1 − bρ > 0 for the covolume equation of state.Source: §4.4 · p. 15
  • Total-energy limiting can be performed as for density, but the authors found it unhelpful and omitted it from the reported numerical tests.Source: §4.4 · p. 15

4.5. Limitation on the internal energy ρe.

The internal-energy limiter bounds the correction from a low-order state using a quadratic auxiliary function. Its admissible step is capped by the smallest positive root, and the proposed limitation is independent of the equation of state.

4.5. Limitation on the internal energy ρe.

  • The internal-energy limiter estimates a step length so that the quasiconcave constraint remains nonnegative along the correction from the low-order state.The section defines an auxiliary function for this constraint.Source: §4.5 · pp. 15–16
  • The auxiliary constraint evaluated along the correction is quadratic, and the limiter uses its smallest positive root as the cap on the step.If there is no positive root, the convention is to set the root value to 1.Source: §4.5 · p. 16
  • The root-based cap keeps the auxiliary constraint nonnegative through the selected interval of the correction.The passage justifies this interval using the smallest-positive-root property.Source: §4.5 · p. 16
  • The proposed limitation on ρe is independent of the equation of state.Source: §4.5 · p. 16
  • Khobalatte and Perthame is identified as the earliest work the authors know of that performs this type of specific-entropy limiting.Source: §4.5 · p. 16

4.6. Minimum principle on the specific entropy.

The specific-entropy minimum is enforced through a concave constraint and a line search along each correction segment. When the direct formulation is unavailable, the method uses a density-weighted entropy constraint; strict enforcement can reduce convergence to first order, motivating relaxation.

4.6. Minimum principle on the specific entropy.

  • The line search selects the largest admissible step along a correction segment, using the boundary of the convex constraint set to determine the limiting value.When the endpoint violates the constraint, the method finds a root; otherwise it takes the full step.Source: §4.6 · p. 16
  • When g satisfies ∂ρρg ≥ 0, the entropy bound is rewritten as ρe(U) − g(smin, ρ) ≥ 0, yielding a concave constraint for the line search.If the full step violates the bound, a combination of secant and Newton methods finds the boundary point.Source: §4.6 · p. 17
  • For general g, the method instead enforces ρΦ(U) − sminρ ≥ 0, whose concavity follows from the convexity of −ρΦ(U).The line search again takes the full step when feasible, or computes the boundary root otherwise.Source: §4.6 · p. 17
  • Strictly enforcing the specific-entropy minimum has been observed to degrade convergence to first order, so the method relaxes the constraint to make its violation second-order accurate.The cited discussion also notes that limiting can reduce accuracy near solution maxima and minima.Source: §4.7 · p. 17

4.7. Relaxation.

The section introduces relaxation variants to retain accuracy while controlling invariant-domain bounds. For density and internal energy, averaging relaxation—not minmod—achieves second-order accuracy in the maximum norm; entropy relaxation is O(h^2) and is unnecessary near shocks.

4.7.1. Relaxation on the density and the internal energy.

  • Only averaging relaxation gives second-order L∞ accuracy; minmod is more restrictive, while without relaxation the method is second-order in L1 for density and internal energy.Both relaxation variants are reported robust.Source: §4.7.1 · p. 18
  • The density and internal-energy relaxation never exceeds 1% of the legitimate lower bound, preserving positivity of both quantities.The definitions use a 0.99ϱmin threshold.Source: §4.7.1 · p. 18

4.7.2. Relaxation on the specific entropy.

  • For specific entropy, relaxation is O(h^2); it is unnecessary near shocks, where first-order viscosity preserves the minimum principle, and smooth isentropic regions have constant entropy.In smooth regions, the local variation measures curvature induced by the nonlinearity of entropy.Source: §4.7.2 · p. 18
  • The reported internal-energy positivity threshold is not reached when the mesh is reasonably fine.The passage states that this threshold guarantees positivity of internal energy.Source: §4.7.2 · p. 18

5. Numerical illustrations.

The illustrations test smooth, rarefaction, shock, vortex, and step flows across three implementations. Reported convergence ranges from second-order accuracy in the maximum norm to first-order shock results.

5.1. Technical details.

  • Three implementations use continuous finite elements and SSP(3,3) time stepping; Code 2's smoothness-based viscosity adds extrema diffusion that degrades L∞ accuracy.The three codes differ in libraries, element families, and high-order viscosity choices.Source: §5.1 · p. 19

5.2. 1D smooth wave.

  • Code 1's δΠ∞ shows O(h^3) superconvergence, while Code 2's error decays as O(h^1.5) under smoothness-indicator clipping.These results are from the one-dimensional smooth-wave test.Source: §5.2 · p. 20, Table 1 · p. 20

5.3. Rarefaction wave.

  • Without relaxing the specific-entropy minimum, the rarefaction convergence rate falls to O(h); Code 1's added dissipation is of truncation-error order.For the maximum norm, the reported rates are O(h) for Code 1 and O(h^0.75) for Code 2.Source: §5.3 · p. 21
  • With density and specific-entropy limitation, the proposed method converges at least O(h^1.5); the low-order method reaches rate 0.96 on 12,800 grid points.Convergence tests for this rarefaction problem are reported for Code 1 and Code 2.Source: §5.3 · p. 21, Table 3 · p. 21

5.4. Leblanc shocktube.

  • For Leblanc shocktube simulations with x0 = 0.33 through t = 2/3 using Code 2 and Code 3, consolidated error indicator δ1(t) converges clearly at first order, optimal for this problem.Table 5’s caption specifies P1 meshes, Code 1 and Code 2, and CFL = 0.25.Source: §5.4 · p. 22, Table 5 · p. 22

5.5. Sod, Lax, Blast wave.

  • For Sod and Lax shocktubes, δ1 convergence rates lie between O(h^0.9) and O(h), described as near optimal.Figure 1 shows density profiles for Sod, Lax, and the Woodward-Collela blast wave.Source: §5.5 · p. 22, Figure 1 · p. 23

5.6. 2D isentropic vortex.

  • Code 3 achieves second-order accuracy in the maximum norm for the isentropic vortex.The reported error indicator is δ∞ at t = 2.Source: §5.6 · p. 23, Table 6 · p. 23

5.7. Mach 3 step.

  • The Mach 3 step density snapshot shows the Kelvin–Helmholtz instability of the contact discontinuity without smoothing the step's top-left corner.The snapshot is from Code 2 at t = 4.Source: §5.7 · p. 23, Figure 2 · p. 24
Loading 1710.00417v1…