### 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 and uniform spacing on some, possibly unbounded, domain .

Definition 17. A difference operator approximating is said to satisfy SBP on the domain with respect to a positive definite scalar product ,

if the property
holds for all grid functions and .

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

for all continuously-differentiable functions and and the scalar product
Similar definitions for SBP can be introduced for higher-dimensional domains.

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

Example 50. Standard centered differences as defined by Eq. (8.19) in the domain or for periodic domains and functions satisfy SBP with respect to the trivial scalar product (),

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

that is, if is diagonal. It is called restricted full if
where denote boundary point indices, or :
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
with and blocks of size independent of .

##### Accuracy and Efficiency.
As mentioned, in the absence of boundaries, standard centered FDs (which have even order of accuracy ) satisfy SBP with respect to the trivial () 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 in the diagonal case and to 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.22) – in general this leads to no solutions. The coefficients of and those of 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 and , respectively.

Example 51. : For the simplest case, , the SBP operator and scalar product are unique:

with and otherwise. That is, the SBP scalar product (8.21) is the simple trapezoidal rule for integration.

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

Example 52. :

with scalar product

On the other hand, the operators have one, three and ten free parameters, respectively. Up to their associated scalar products are unique, while for one of the free parameters enters in . For the full-restricted case, 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 . If the difference operator in the interior is a standard centered difference with accuracy-order then there are points at and near each boundary, where the accuracy is of order (with in the diagonal case and in the full restricted one). The integer can be referred to as the boundary width. The boundary stencil size is the number of gridpoints that the difference operator uses to evaluate its approximation at those 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, 281]. It turns out that in this way the order of accuracy can be increased from the very low one of to higher-order ones such as or with a very small change in the CFL limit. It involves some work, but since the SBP property (8.22) 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

then
Table 6 illustrates the results of this optimization procedure for the operator; see [141] for more details.

The coefficients for the SBP operators

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 [141] and also as complete source code from the Einstein Toolkit [151].

Table 6: Comparison, for the 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, , (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:

• The requirement of uniform spacing is not an actual restriction, since a coordinate transformation can always be used so that the computational grid is uniformly spaced even though the physical distance between gridpoints varies. In fact, this is routinely done in the context of multiple-domains or curvilinear coordinates (see Section 11). In that case, though, stability needs to be guaranteed for systems with variable coefficients, since they appear due to the coordinate transformation(s) even if the original system had constant coefficients. This has relevance in terms of the distinction between diagonal and block-diagonal SBP norms, as mentioned below.
• A similar concept of SBP holds for discrete expansions into Legendre polynomials using Gauss-type quadratures, as discussed in Section 9.4.
• The definition of SBP depends only on the computational domain, not on the system of equations being solved. This allows the construction of optimized SBP operators once and for all.
• Difference operators satisfying SBP, which are genuinely multi-dimensional, can be explicitly constructed (see, for example, [103, 102]). However, they become rather complicated even for simple geometries as higher-order accuracy is sought. An easier approach, for the case in which the domain is the cross product of one dimensional ones (say, topologically a cube in three dimensions), which is usually the case in many domain-decompositions for complex geometries (Section 11) is to simply apply a one-dimensional operator satisfying SBP in each direction, and this is the approach that we will discuss from now on. The question then is whether SBP holds in several dimensions; the answer is affirmative in the case of diagonal norms but not necessarily otherwise.