8.3 Summation by parts

Since numerical stability is, by definition, the discrete counterpart of well-posedness, one way to come up with schemes that are, by construction, numerically stable, is by designing them so that they satisfy the same properties used at the continuum when showing well-posedness through an energy estimate. As discussed in Section 3.2.3 one such property is integration by parts or the application of Gauss’ theorem, which leads to its numerical counterpart: SBP [265, 266].

Consider a discrete grid consisting of points {xi}Ni=0 and uniform spacing Δx on some, possibly unbounded, domain [a,b].

Definition 17. A difference operator D approximating ∂ ∕∂x is said to satisfy SBP on the domain [a,b] with respect to a positive definite scalar product Σ,

( ) Σ = Δx ( σij ) , (8.20 )
∑N ⟨u,v ⟩Σ := Δx uivjσij, (8.21 ) i,j=0
if the property
⟨u,Dv ⟩Σ + ⟨Du, v⟩Σ = u(b)v(b) − u(a)v(a) (8.22 )
holds for all grid functions u and v.

This is the discrete counterpart of integration by parts for the ddx operator,

-d- -d- ⟨f,dx g⟩ + ⟨dx f,g⟩ = f(b)g(b) − f(a)g(a), (8.23 )
for all continuously-differentiable functions f and g and the scalar product
∫ b ⟨f, g⟩ := f(x)g(x )dx. (8.24 ) a
Similar definitions for SBP can be introduced for higher-dimensional domains.

If the interval is infinite, say (− ∞, b) or (− ∞, ∞ ), certain fall-off conditions are required and Eq. (8.22View Equation) replaced by dropping the corresponding boundary term(s).

Example 50. Standard centered differences as defined by Eq. (8.19View Equation) in the domain (− ∞, ∞ ) or for periodic domains and functions satisfy SBP with respect to the trivial scalar product (σij = δij),

∑ ⟨u, D0v ⟩Σ + ⟨D0u, v⟩Σ = 0, ⟨u,v⟩Σ := Δx uivi. (8.25 ) i∈ ℤ

The scalar product or associated norm are said to be diagonal if

σij = σiiδij, (8.26 )
that is, if Σ is diagonal. It is called restricted full if
σibj = σibibδibj, (8.27 )
where ib denote boundary point indices, ib = 0 or ib = N:
( ) σ00 0 ⋅⋅⋅ 0 | | || 0. . || Σ = Δx | .. Σinterior .. | . (8.28 ) |( 0 |) 0 ⋅⋅⋅ 0 σNN
In the case of bounded, non-periodic domains, one possibility is to use centered differences and the trivial scalar product in the interior and modify both of them at and near boundaries. That is, the scalar product has the form
( ) Σl Σ = Δx ( I ) , (8.29 ) Σr
with Σl and Σr blocks of size independent of Δx.

Accuracy and Efficiency.
As mentioned, in the absence of boundaries, standard centered FDs (which have even order of accuracy 2p) satisfy SBP with respect to the trivial (Σ = ΔxI) scalar product. In their presence the operators can be modified at and near boundaries so as to satisfy SBP, examples are given below. It can be seen that the accuracy at those points drops to p in the diagonal case and to 2p − 1 in the restricted full one. Therefore, the latter is more desirable from an accuracy perspective, but less so from a stability one, as we will discuss at the end of this subsection. Depending on the system, numerical dissipation might be enough to stabilize the discretization in the restricted full case. This is discussed below in Section 8.5.

When constructing SBP operators, the discrete scalar product cannot be arbitrarily fixed and afterward the difference operator solved for so that it satisfies the SBP property (8.22View Equation) – in general this leads to no solutions. The coefficients of Σ and those of D have to be simultaneously solved for. The resulting systems of equations lead to SBP operators being in general not unique, with increasing freedom with the accuracy order. In the diagonal case the resulting norm is automatically positive definite but not so in the full-restricted case.

We label the operators by their order of accuracy in the interior and near boundary points. For diagonal norms and restricted full ones this would be D2p −p and D2p −(p− 1), respectively.

Example 51. D2 −1: For the simplest case, p = 1, the SBP operator and scalar product are unique:

( | (− u0 + u1) for i = 0, 1 |{ ( 1 1 ) D2 −1ui = ---- -ui+1 − --ui− 1 for i = 1 ...(N − 1), (8.30 ) Δx ||( 2 2 (un − un−1) for i = N ,
with σ00 = σNN = 1∕2 and σij = δij otherwise. That is, the SBP scalar product (8.21View Equation) is the simple trapezoidal rule for integration.

The operator D4 −2 and its associated scalar product are also unique in the diagonal norm case:

Example 52. D4 −2:

( ( ) | 24 59 4 3 |||| − ---u0 + --u1 − --u2 − ---u3 for i = 0, ||| ( 17 34) 17 34 ||| − 1-u0 + 1u2 for i = 1, |||| ( 2 2 ) ||| 4-- 59- 59- 4-- ||| 43u0 − 86 u1 + 86u3 − 43u4 for i = 2, |||| ( 3 59 32 4 ) ||| --u0 − ---u2 + --u4 − --u5 for i = 3, ||| 98 98 49 49 |||{ ( ) -1-- 1-- 2- 2- 1-- D4−2ui = Δx | 12ui−2 − 3ui−1 + 3 ui+1 − 12ui+2 for i = 4...(N − 4), (8.31 ) |||| ||| ( ) ||| − 3-uN − 59uN −2 + 32-uN −4 − 4-uN −5 for i = N − 3, |||| ( 98 98 49 49 ) ||| 4-- 59- 59- 4-- ||| − 43uN − 86uN −1 + 86 uN −3 − 43uN −4 for i = N − 2, |||| ( 1 1 ) ||| − − --uN + -uN −2 for i = N − 1, ||| ( 2 2 ) |||( − − 24-u + 59u − -4-u − 3-u for i = N , 17 N 34 N−1 17 N− 2 34 N −3
with scalar product
{ } 17 59 43 49 49 43 59 17 Σ = Δx diag 48-,48, 48,48-,1,1,...,1,1,48-,48-,48, 48- . (8.32 )

On the other hand, the operators D6 −3,D8 −4,D10−5 have one, three and ten free parameters, respectively. Up to D 10−5 their associated scalar products are unique, while for D 10−5 one of the free parameters enters in Σ. For the full-restricted case, D4 −3,D6 −5,D8 −7 have three, four and five free parameters, respectively, all of which appear in the corresponding scalar products.

A possibility [396] is to use the non-uniqueness of SBP operators to minimize the boundary stencil size s. If the difference operator in the interior is a standard centered difference with accuracy-order 2p then there are b points at and near each boundary, where the accuracy is of order q (with q = p in the diagonal case and q = 2p − 1 in the full restricted one). The integer b can be referred to as the boundary width. The boundary stencil size s is the number of gridpoints that the difference operator uses to evaluate its approximation at those b boundary points.

However, minimizing such size, as well as any naive or arbitrary choice of the free parameters, easily leads to a large spectral radius and as a consequence restrictive CFL (see Section 7) limit in the case of explicit evolutions. Sometimes it also leads to rather large boundary truncation errors. Thus, an alternative is to numerically compute the spectral radius for these multi-parameter families of SBP operators and find in each case the parameter choice that leads to a minimum [399, 281Jump To The Next Citation Point]. It turns out that in this way the order of accuracy can be increased from the very low one of D2 −1 to higher-order ones such as D10− 5 or D8 −7 with a very small change in the CFL limit. It involves some work, but since the SBP property (8.22View Equation) is independent of the system of equations one wants to solve, it only needs to be done once. In the full-restricted case, when marching through parameter space and minimizing the spectral radius, this minimization has to be constrained with the condition that the resulting norm is actually positive definite.

The non-uniqueness of high-order SBP operators can be further used to minimize a combination of the average of the boundary truncation error (ABTE), defined below, without a significant increase in the spectral radius. For definiteness consider a left boundary. If a Taylor expansion of the FD operator is written as

| | du | q dq+1u | Du |xi = ---|| + ci(Δx ) ---q+1-|| for i = 0,1,...,b, (8.33 ) dx xi dx xi
then
( ) ∑ b 1∕2 ABTE := 1- c2i . (8.34 ) b i=0
Table 6 illustrates the results of this optimization procedure for the D10 −5 operator; see [141Jump To The Next Citation Point] for more details.

The coefficients for the SBP operators

D2 −1,D4 −2,D4 −3,D6− 3,D6− 5,D8 −4,D8 −7,D10 −5, (8.35 )
and, in particular, for their optimized versions, are available, along with their associated dissipation operators described below in Section 8.5, from the arXiv in [141Jump To The Next Citation Point] and also as complete source code from the Einstein Toolkit [151].


Table 6: Comparison, for the D10−5 operator, of both the spectral radius and average boundary truncation error (ABTE) when minimizing the bandwidth or a combination of the spectral radius and ABTE. For comparison, the spectral radius and ABTE for the lowest-accuracy operator, D1 −2, (which is unique) are 1.414 and 0.25, respectively. Note: the ABTE, as defined, is larger for this operator, but its convergence rate is faster.
Operator Min. bandwidth Min. ABTE and spectral radius
Spectral radius 995.9 2.240
ABTE 20.534 0.7661

Remarks:


  Go to previous page Go up Go to next page