PaperlayerView exact revision on arXiv

Source-linked AI summary

Targeted Entropy-Stable Nonlinear Artificial Dissipation for High-Order Spectral-Element Summation-By-Parts Discretizations

Zelalem Arega Worku, Alex Bercik, David C. Del Rey Fernández, David W. Zingg

arXiv:2610.17567v1math.NAcs.NA

TL;DR

Under-resolved high-order methods need added damping, but excessive dissipation can harm accuracy and suppress relevant structures. This paper introduces sensor-activated nonlinear dissipation for SBP discretizations; its experiments report improved robustness while retaining high-order accuracy for smooth solutions.

  • Problem

    Artificial dissipation can improve robustness for under-resolved flows, but excessive dissipation can degrade accuracy and suppress physically relevant structures.

    Source: §1 · p. 2

  • Method

    The paper adds an element-local nonlinear dissipation term with a conservative-variable distribution, a separately specified entropy-viscosity budget, and an operator-normalized smoothness sensor.

    Source: §1 · p. 3

  • Results

    All C-SBP and D-SBP Kelvin–Helmholtz computations with p ≤7 reach final time, whereas every D-SBP computation with LLF interface dissipation alone fails; smooth-solution tests retain high-order convergence.

    Source: §7 · p. 36

  • Takeaways & Limitations

    The added term preserves conservation and semidiscrete entropy stability whenever the base method already has those properties.

    Source: §7 · p. 35

    The dissipation does not guarantee positivity: runs with p ≥8 fail in Kelvin–Helmholtz tests, and the tested entropy-stable discretizations fail the two-dimensional Riemann problem because of positivity violations.

    Source: §7 · p. 36

Abstract

from arXiv · show

Artificial dissipation improves the robustness of high-order methods for under-resolved flows, but can compromise accuracy and suppress physically relevant structures if its strength and spatial distribution are not appropriately controlled. We present a sensor-activated nonlinear artificial-dissipation framework, inspired by the residual correction method, for continuous and discontinuous summation-by-parts discretizations. In contrast to many existing spectral-element dissipation operators that either reduce the formal order of accuracy or target only the highest polynomial modes, the proposed framework permits design- and higher-order dissipation with explicit control over its magnitude, spatial distribution, activation, and convergence rate under mesh refinement. This modular construction also recovers several existing methods through appropriate choices of its components. The dissipation preserves discrete conservation and semidiscrete entropy stability whenever the base method satisfies these properties. It is dimensionally consistent, element local, and discretization agnostic, as the added term is constructed independently of the base spatial residual. For smooth solutions under suitable mesh assumptions, we establish parameter choices that retain the formal accuracy of the base discretization. Numerical experiments with summation-by-parts discretizations of the Burgers and Euler equations demonstrate substantially improved robustness for under-resolved flows and flows with shock waves while retaining high-order accuracy for smooth solutions, using the same parameter choices across several test problems. Overall, these results show that targeted dissipation can improve the robustness of high-order methods without excessive damping of resolved flow structures.

1. Introduction

The introduction motivates controlled artificial dissipation as a way to stabilize under-resolved high-order methods without sacrificing accuracy or physically relevant structures. It presents an element-local nonlinear framework and outlines tests of its accuracy, stability, and robustness.

1. Introduction

  • Artificial dissipation can damp under-resolved modes and improve robustness, but excessive strength can degrade accuracy, suppress relevant structures, and reduce the stable time step.The introduction frames spatial targeting, magnitude, and mesh-refinement decay as controls for balancing these effects.Source: §1 · p. 2
  • Existing spectral-element approaches can require extra nodes or relax polynomial exactness, while the proposed construction supplies smooth-state decay without high-order annihilation of its distribution operator.On minimally nodal elements, the cited Bercik et al. construction acts only on the highest-degree mode and sacrifices one order of accuracy.Source: §1 · pp. 2–3
  • The added term has zero conservative moment and nonpositive entropy contribution, preserving conservation and semidiscrete entropy stability when the base method has those properties.It does not generally make a non-entropy-stable base method entropy stable.Source: §1 · p. 3
  • The proposed element-local dissipation combines a conservative-variable distribution, a nonnegative entropy-viscosity budget, and an operator-normalized smoothness sensor.Its nonlinear activation provides adjustable high-order decay under mesh refinement, and its formulation is independent of the base spatial residual.Source: §1 · p. 3
  • Experiments assess smooth-solution accuracy, under-resolved vortical flows, shocks, and kinetic-energy evolution for Burgers and Euler discretizations using the same parameter-selection rules.The test set includes isentropic vortices, Kelvin–Helmholtz instabilities, Burgers shock problems, Sod, Shu–Osher, two-dimensional Riemann, and Taylor–Green tests.Source: §1 · p. 3

2. Preliminaries

The preliminaries define the hyperbolic conservation-law and entropy setting, then establish the mesh, notation, and SBP-operator conventions used in the paper.

2.1. Hyperbolic PDE and entropy identities

  • A strictly convex mathematical entropy defines entropy variables, while the entropy flux and entropy identity characterize the continuous entropy balance for smooth solutions.For entropy solutions with shocks, the entropy equality becomes an inequality; periodic boundaries cancel the boundary contributions.Source: §2.1 · p. 4

2.2. Mesh and SBP operators

  • The mesh is conforming, quasi-uniform, and shape regular, and the element length scale h is the maximum element diameter.The paper uses physical SBP scaling on these meshes and specifies the dimensions of the mass, derivative, and boundary operators.Source: §2.2 · pp. 5–6
  • The diagonal-norm SBP operators satisfy the SBP property and use collocated volume and facet nodes.The preliminaries also define componentwise identity notation and distinguish assembled shared degrees of freedom for C-SBP from independently coupled D-SBP element vectors.Source: §2.2 · pp. 5–6

3. Base spatial discretization

This section establishes the base residuals and their conservation and entropy-stability requirements. It describes the Hadamard-form two-point-flux SBP setting and the different coupling mechanisms for discontinuous and continuous discretizations.

3. Base spatial discretization

  • The analysis focuses on conservative entropy-stable base residuals, especially Hadamard-form two-point-flux SBP discretizations, while also applying the dissipation to other Burgers and Euler discretizations.The experiments use these Hadamard-form discretizations as the default Euler bases.Source: §3 · p. 6
  • For D-SBP, element-local two-point-flux residuals couple across facets, whereas C-SBP assembles local volume contributions whose interior boundary terms cancel pairwise.A symmetric consistent two-point flux satisfies Tadmor’s entropy-consistency condition for entropy conservation.Source: §3 · pp. 6–7
  • D-SBP can add interface dissipation as a conservative SAT, but C-SBP lacks double-valued interface states and instead uses volume or edge stabilization.The cited C-SBP alternatives can require additional nodes, relax polynomial exactness, or use coefficient scaling; the proposed sensor-based dissipation is presented as targeting a wider range of modes without compromising accuracy in the numerical comparisons considered.Source: §3 · pp. 7–8

4. Nonlinear high-order artificial dissipation

This section develops an element-local, sensor-activated dissipation framework that adds a prescribed entropy budget while preserving conservation and entropy stability, with components controlling its spatial action, activation, and mesh-refinement decay. It also presents the cheaper entropy-viscosity budget used in the numerical results, based on elementwise scaling, a cell-average entropy Jacobian, and reference-element derivatives.

4.1. A general algebraic structure for entropy-stable artificial dissipation

  • The general construction adds an element-local dissipation term with zero conservative moment and a prescribed nonnegative entropy-production budget.A distribution vector controls spatial placement, while a scalar budget sets the amount of dissipation; the budget combines a smoothness sensor, an entropy-viscosity factor, and an optional magnitude parameter.Source: §4.1 · pp. 9–10

4.2. Distribution vector

  • The derivative-based distribution is built from conservative variables; a graph-Laplacian safeguard is optional but unused because its robustness-versus-resolution tradeoff depends on the problem.Exploratory comparisons found improved Kelvin–Helmholtz robustness but stronger attenuation in Shu–Osher and Taylor–Green; the study retains the derivative-based construction.Source: §4.2 · pp. 10–11

4.3. Solution smoothness sensor

  • The novel bounded sensor measures high-order derivative content in dimensionlessly scaled conservative variables, vanishes on constant states, and is O(h^sθ) in smooth regions but can be O(1) for under-resolved data.Scaling avoids unit-dependent component weights without changing the polynomial degree detected by the reference derivative operator.Source: §4.3 · pp. 12–13

4.3.1. Operator-normalization

  • Operator normalization divides the sensor’s Rayleigh quotient by a generalized-eigenvalue upper bound, yielding a value below one and removing operator-dependent scaling.The bound is computed once for each reference operator, polynomial degree, and sensor order.Source: §4.3.1 · p. 13

4.3.2. Removing the constant state

  • The mean-subtracted sensor compares local oscillatory variation with grid resolution, ignores fluctuations below a resolution-based threshold, and decays as O(h^2⌈p/2⌉) in smooth regions for p ≥ 2.It responds to local variation rather than classifying shocks, shear layers, physical structures, or numerical oscillations separately; for p = 1, the normalized operator sensor is used.Source: §4.3.2 · pp. 13–14

4.3.3. Adaptive dissipation order

  • A Hill-function indicator adapts the sensor exponent so weak activations decay at high order while strong activations retain low-order decay; the study sets pL = p and pH = 2p.The mapping distinguishes shocks and large oscillations, small under-resolved structures, and very smooth oscillations requiring little dissipation.Source: §4.3.3 · p. 15

4.4. Entropy-viscosity budget

  • The physical reference entropy budget approximates scalar viscosity using directional nodal viscosity coefficients and a positive-definite approximation to the entropy Jacobian; the continuous viscosity has dimension L2/T.More sophisticated matrix viscosities are deferred to future work.Source: §4 · p. 9, §4.4 · p. 16

4.4.1. Cheaper entropy-viscosity budget

  • Because the exact continuous viscosity magnitude is unknown, the cheaper budget simplifies the physical construction by replacing nodal viscosity values with a constant coefficient.The constant coefficient follows the viscosity scaling relation and is intended to retain a similar budget structure.Source: §4.4.1 · p. 16
  • The approximation uses λ_k h_eff,k as a nodal viscosity scale and a constant entropy-Jacobian matrix based on a cell-average state, with derivatives computed on the reference element.The scaling preserves units and mesh-size scaling but not directional geometric information.Source: §4.4.1 · p. 16
  • The implemented budget is nonnegative, vanishes for constant entropy variables, and has the same physical dimensions as the entropy-production budget.The paper uses this approximation for all numerical results; implementation applies reference derivative matrices and the cell-average entropy-Jacobian matrix, then reduces over nodes and coordinate directions.Source: §4.4.1 · pp. 16–17
  • For highly skewed elements, metric terms at each node can provide local scaling instead of explicit h_eff,k and μ_k factors.This alternative approximates directional viscosity with the same units and mesh-size scaling as the implemented budget.Source: §4.4.1 · p. 17

5. Properties of the dissipation operator

The dissipation has the same dimensions as the time derivative, retains high-order accuracy for smooth solutions under stated mesh and solution assumptions, and is conservative and entropy dissipative. Its low- and high-order parameters let designers trade stronger active-element dissipation against smooth-region accuracy without changing conservation or entropy stability.

5.1. Dimensional consistency

  • The artificial-dissipation term has the same physical dimension as ∂U/∂t.The analysis establishes this for the overall dissipation operator.Source: §5.1 · p. 17

5.2. High-order accuracy in smooth regions

  • Under shape-regular, quasi-uniform meshes with affine mappings and smooth admissible solutions, parameter choices recover high-order dissipation scaling as the smoothness indicator vanishes.The result distinguishes low-order consistency for any indicator value from high-order behavior as the indicator vanishes.Source: §5.2 · pp. 17, 18, 20
  • The low-order parameter pL controls worst-case consistency and may be set below p to strengthen dissipation in active elements, while pH governs resolved smooth elements.Design-order consistency follows from pL ≥ p or, with the stated indicator, pH ≥ p; these choices do not affect conservation or entropy stability.Source: §5.2 · p. 20

5.3. Conservation and entropy stability

  • The assembled dissipation is conservative and contributes nonpositively to entropy, so it preserves semidiscrete entropy stability when the base residual satisfies the corresponding entropy condition.The local contribution is also conservative and entropy dissipative.Source: §5.3 · p. 20

6. Numerical experiments

The experiments assess accuracy, conservation, entropy dissipation, shock performance, and robustness for under-resolved vortical flows using continuous and discontinuous SBP discretizations. They use common dissipation parameters across test cases, with specified base discretizations and time-integration settings.

6.1.1. Accuracy

  • For smooth isentropic vortices, C-SBP errors with and without artificial dissipation agree closely, while solution-error rates approach the expected p- to (p + 1)-order range in both dimensions.The maximum elementwise entropy budget decreases at least at order 2p + d + 1 under refinement.Source: §6.1.1 · p. 23, Table 1 · p. 23, Table 2 · p. 24

6.1.2. Conservation, entropy stability, and long-time error behavior

  • The added dissipation preserves conservation to machine precision and yields nonpositive entropy changes with RRK4, while smooth-solution long-time errors do not increase significantly.For D-SBP, adding artificial dissipation to LLF produces essentially no further change in entropy or long-time error for these smooth tests.Source: §6.1.2 · pp. 24–25, Figure 1 · p. 24, Figure 3 · p. 25

6.2. One-dimensional Burgers shock problem

  • Both C-SBP and D-SBP capture the Burgers shock at the correct location, with the entropy budget sharply localized at the shock despite localized oscillations near the discontinuity.The dissipation is not designed as a rigorous shock-capturing mechanism.Source: §6.2 · pp. 26–27, Figure 5 · p. 27

6.3. Sod shock tube

  • For the Sod shock tube, both formulations capture the rarefaction, contact, and shock with minimal oscillations and remain close to the exact Riemann solution.The budget concentrates near the shock; for D-SBP, adding the proposed dissipation reduces discontinuity-adjacent oscillations relative to LLF alone.Source: §6.3 · pp. 27–28, Figure 6 · p. 28

6.4. Shu–Osher shock tube

  • For Shu–Osher, both formulations with the proposed dissipation reach tf = 1.8, while D-SBP with LLF alone develops much larger oscillations at p = 4 and terminates early at p = 8.At fixed nedge = 100, p = 8 better matches the reference post-shock density oscillations than p = 4.Source: §6.4 · p. 29, Figure 8 · p. 30

6.5. Kelvin–Helmholtz instability

  • 21 of 27 Kelvin–Helmholtz cases reach tf = 15 for each formulation with the proposed dissipation; all split-form D-SBP runs with LLF alone fail by t = 5.34.The remaining failures with the proposed dissipation are confined to p ≥8; representative density fields retain characteristic roll-up and vortical structures.Source: §6.5 · p. 31, Table 3 · p. 31

6.6. Two-dimensional Riemann problem

  • Both C-SBP and D-SBP reach the final time in the two-dimensional Riemann test and produce closely agreeing density fields, with entropy budgets concentrated near shocks.The test uses a direct divergence-form base discretization because the entropy-stable schemes did not complete the simulation due to positivity requirements.Source: §6.6 · pp. 31–32, Figure 9 · p. 32

6.7. Three-dimensional inviscid Taylor–Green vortex

  • In the Taylor–Green vortex, D-SBP histories with and without artificial dissipation initially nearly coincide, but later the artificial-dissipation histories show larger peak rates and greater total kinetic-energy decay than LLF alone.For both p = 3 and p = 4, C-SBP with artificial dissipation has a lower initial dissipation rate than D-SBP with LLF for t ≲3.Source: §6.7 · p. 35, Figure 12 · p. 35

7. Conclusion

The paper develops element-local dissipation with independently controlled design features and conservation and entropy-stability guarantees. Tests show improved robustness with limited damping of resolved structures, while positivity remains an important boundary.

7. Conclusion

  • The element-local dissipation has separately controlled spatial distribution, magnitude, activation, and mesh-refinement convergence rate, independent of the base residual.Its added term has zero conservative moment and a nonpositive entropy contribution, preserving conservation and semidiscrete entropy stability when the base method has those properties.Source: §7 · p. 35
  • For smooth admissible solutions on affine, shape-regular, quasi-uniform meshes, the dissipation is O(h^pL), improving to O(h^pH) when p ≥2 and pL ≤pH.Suitable order choices retain the base discretization’s formal accuracy.Source: §7 · p. 36
  • The tests show high-order smooth-solution accuracy and substantially improved under-resolved-flow robustness without excessive damping of resolved structures.All C-SBP and D-SBP Kelvin–Helmholtz computations with p ≤7 reach final time, while every D-SBP run with LLF interface dissipation alone fails; the same parameter-selection rules are used across tests.Source: §7 · p. 36
  • The dissipation does not guarantee positivity: Kelvin–Helmholtz runs with p ≥8 and entropy-stable schemes in the two-dimensional Riemann problem still fail from positivity violations.The paper identifies parameter adaptation and combination with positivity-preserving limiting as future directions.Source: §7 · p. 36

Appendix A. Relation to Existing Volume Dissipation Schemes

The appendix recalls the general volume-dissipation structure and identifies the elementwise dissipation budget and positive quantities used in it.

Appendix A. Relation to Existing Volume Dissipation Schemes

  • The appendix begins by recalling the general structure of the volume dissipation.The cited passage identifies that structure as equations (26) and (27).Source: Appendix A · p. 36
  • The elementwise budget bν,k specifies the total amount of volume dissipation in each element.Source: Appendix A · p. 36
  • The appendix states that ak and ϵmach are positive.Source: Appendix A · p. 36

Appendix A.1. Dissipation of Bercik, Craig Penner, and Zingg

The appendix relates the proposed construction to the Bercik, Craig Penner, and Zingg scheme by comparing distribution order, smooth-state decay control, and variable choice.

Appendix A.1. Dissipation of Bercik, Craig Penner, and Zingg

  • The referenced derivative operator is equivalent up to a constant to the reference-element difference operator used here, with reference-to-physical scaling handled implicitly by metric terms in Ci.Setting θk = 1 is part of placing that operator in the framework.Source: Appendix A.1 · p. 37
  • The referenced scheme’s distribution order sa = p restricts dissipation to the highest-degree mode on minimally nodal spectral elements and relaxes polynomial exactness by one degree.Source: Appendix A.1 · p. 37
  • The proposed method instead fixes sa = 1 and uses sensor order sθ with θβk to obtain prescribed smooth-state decay.Source: Appendix A.1 · p. 37
  • The proposed method uses uk in ak for accuracy and wk in bν,k to mimic an entropy viscosity, unlike the referenced scheme’s shared entropy-variable choice [6].Source: Appendix A.1 · p. 37

Appendix A.2. Local Projection Stabilization

The appendix describes local projection stabilization and explains how the proposed framework can retain entropy dissipation while using a different distribution and activation strategy.

Appendix A.2. Local Projection Stabilization

  • The LPS adaptation penalizes high-order polynomial degrees, with Π isolating entropy-variable content outside the degree-p polynomial space.Source: Appendix A.2 · p. 37
  • The framework need not apply the projection directly to entropy variables to ensure entropy stability, unlike the cited LPS approach.Source: Appendix A.2 · p. 37
  • A first-derivative distribution acting on lower-degree variations, paired with a nonlinear budget, avoids the zero-projection restriction on minimally nodal elements while supplying smooth-state decay.Source: Appendix A.2 · p. 37

Appendix A.3. Residual-based artificial viscosity

The appendix recasts residual-based viscosity as a smoothness sensor within the framework, with an entropy budget and distribution that preserve conservative, entropy-dissipative structure. It also relates this construction to scalar residual-viscosity diffusion and to the more specialized viscosity of.

Appendix A.3. Residual-based artificial viscosity

  • The compared residual-viscosity formulations separate first-order and high-order viscosities, using wavespeed, mesh spacing, residual, and normalization in the coefficient.The viscosity coefficient multiplies an appropriate viscous flux.Source: Appendix A.3 · p. 38
  • The residual-viscosity coefficient is interpreted primarily as a smoothness sensor, with its first-order viscosity used to compute the entropy budget.The construction uses θk directly as activation when βk = 1.Source: Appendix A.3 · p. 38
  • The distribution uses an upwind SBP operator and a boundary modification moved inside the operator to ensure conservation.The boundary modification facilitates boundary-condition imposition.Source: Appendix A.3 · p. 38
  • For scalar problems, the unregularized construction reduces to scalar residual-viscosity diffusion; for nonlinear systems, entropy normalization and a positive entropy budget modify monolithic conservative-variable residual viscosity.The physics-based viscosity in treats Navier–Stokes components separately and would require further generalization.Source: Appendix A.3 · pp. 38–39
  • Nodal residual activation can be incorporated through the distribution and budget, retaining conservative, entropy-dissipative structure and the spatial localization of nodal residual-viscosity methods.The elementwise scalar activation controls the overall amount of residual activation.Source: Appendix A.3 · p. 39

Appendix B. Implementation Summary

The implementation summary links core quantities to their related equations and records the construction of the distribution vector.

Appendix B. Implementation Summary

  • The summary associates d_k = −α_k,ε ε_k with equations (27) and (29), respectively.It also lists b_k with equation (28).Source: Appendix B · p. 39
  • The implementation summary identifies V_k as related to equation (9), alongside a_k from equation (31).The summary points to equation (30) for the matrix construction.Source: Appendix B · p. 39
  • Equation (31) defines a_k = M_k u_k, with M_k constructed from μ_k and the matrices shown in the implementation summary.The summary cross-references equations (30) and (32) for the related quantities.Source: Appendix B · p. 39

Artificial-dissipation parameter choices

The appendix lists default artificial-dissipation parameter choices, author contributions, and the authors’ declaration about generative AI assistance.

Artificial-dissipation parameter choices

CRediT authorship contribution statement

  • Worku led conceptualization, methodology, analysis, software, investigation, validation, visualization, and writing; Bercik contributed to methodology, analysis, investigation, software, validation, and writing.Del Rey Fernández and Zingg are credited with supervision, validation, and review and editing.Source: CRediT authorship contribution statement · p. 40

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

Loading 2610.17567v1…