8.5 Numerical dissipation

The use of numerical dissipation consistently with the underlying system of equations is a standard way of filtering unresolved modes, stabilizing a scheme, or both, without spoiling the convergence order of the scheme. As an example of unresolved modes, for centered differences, the mode with highest frequency for any given resolution does not propagate at all, while the semi-discrete group velocity with highest frequency is exactly in the opposite direction to the one at the continuum. In addition, the speed increases with the order of the scheme. See, for example, [281Jump To The Next Citation Point], for more details.

Some schemes, such as those with upwind FDs, are intrinsically dissipative, with a fixed “amount” of dissipation for a given resolution. Another approach is to add to the discretization a dissipative operator Q d with a tunable strength factor 𝜖 ≥ 0,

u˙= (...) → ˙u = (...) + 𝜖Qdu. (8.51 )
The operator Qd has derivatives of higher order than the ones in the principal part of the equation, mimicking dissipative physical systems and/or parabolic equations but in such a way that ∥Qd ∥ → 0 as Δx → 0. Furthermore, Qd is usually chosen so that Qd scales with the gridspacing as the highest FD approximation (so that the amplification factor depends only on Δt ∕Δx). For example, for first-order-in-space systems, FDs scale as
-1-- D ∼ Δx (8.52 )
and Qd is usually chosen to scale in the same way. More precisely, in the absence of boundaries the standard way to add numerical dissipation to a first-order-in-space system is through Kreiss–Oliger dissipation
r−1 2r− 1 (r) (r) Qd = (− 1) (Δx ) D + D − , (8.53 )
where D (r) denotes the application of D r times and D ,D + − denote forward and backward one-sided FDs, respectively, as defined in Eq. (7.70View Equation). Thus, Qd scales with the gridspacing as −1 (Δx ), like Eq. (8.52View Equation). If the accuracy order of the scheme is not higher than (2r − 1 ) in the absence of dissipation, it is not decreased by the addition of numerical dissipation of the form (8.53View Equation).

The main property that sometimes allows numerical dissipation to stabilize otherwise unstable schemes is when they strictly carry away energy (as in the energy definitions involved in well-posedness or numerical-stability analysis) from the system. For example, the operators (8.53View Equation) are semi-negative definite

⟨u,Qdu ⟩Σ ≤ 0, (8.54 )
with respect to the trivial scalar product Σ = Δx 𝕀, under which centered differences satisfy SBP.

In the presence of boundaries, it is standard to simply set the operators (8.53View Equation) to zero near them. The result is, in general, not semi-negative definite as in (8.54View Equation), which cannot only not help resolve instabilities but also trigger them. Many times this is not the case in practice if the boundary is an outer one, where the solution is weak, but not for inter-domain boundaries (see Section 10). For example, for a discretization of the standard wave equation on a multi-domain, curvilinear grid setting, using the D6− 5 SBP operator with Kreiss–Oliger dissipation set to zero near interpatch boundaries does not lead to stability while the more elaborate construction below does [141Jump To The Next Citation Point].

For SBP-based schemes, adding artificial dissipation may lead to an unstable scheme unless the dissipation operator is semi-negative under the SBP scalar product. In addition, the dissipation operator should ideally be non-vanishing all the way up to the boundary and preserve the accuracy of the scheme everywhere (which is more difficult in the SBP case, as it is non-uniform). In [303Jump To The Next Citation Point], a prescription for operators satisfying both conditions for arbitrary–high-order SBP scalar products, is presented. A compatible dissipation operator is constructed as

Qd = − (Δx )2pΣ −1DTp BpDp, (8.55 )
where Σ is the SBP scalar product, Dp is a consistent approximation of p p d ∕dx with minimal bandwidth (other choices are presumably possible), and Bp is called the boundary operator. The latter has to be positive semi-definite and its role is to allow boundary points to be treated differently from interior points. Bp cannot be chosen freely, but has to follow certain restrictions (which become somewhat involved in the non-diagonal SBP case) based on preserving the accuracy of the schemes near and at boundaries; see [303] for more details.
  Go to previous page Go up Go to next page