presents a comparison, for the FD case, of numerical boundary conditions through injection, orthogonal projections, and penalty terms. One additional advantage of the penalty approach is that, for advection-diffusion problems, in time it damps away violations of the compatibility conditions between the initial and boundary data; see Section 5. Figure 7 shows the error as a function of time for an advection-diffusion equation where the initial data is perturbed so that an inconsistency with the boundary condition is present, and the error is defined as the difference between the unperturbed and perturbed solutions. See  for more details. One expects this not to be the case for hyperbolic problems, though. For example, at the continuum an incompatibility between the initial data and initial boundary condition for the advection would propagate, not dissipate in time, and a consistent and convergent numerical scheme should reflect this behavior.
Many systems in numerical relativity, starting with Einstein’s equations themselves, are numerically solved by reducing them to first-order systems, because there is a large pool of advanced and mature analytical and numerical techniques for them. However, this is at the expense of enlarging the system and, in particular, introducing extra constraints (though this seems to be less of a concern; see, for example, [286, 82]). It seems more natural to solve such equations directly in second-order (at least in space) form. It turns out, though, that it is considerably more complicated to ensure stability for such systems. Trying to “integrate back” an algorithm for a first-order reduction to its original second-order form is, in general, not possible; see, for example, [101, 413] for a discussion of this apparent paradox and the difficulties associated with constructing boundary closures for second-order systems such that an energy estimate follows.
The SAT approach was generalized to a class of wave equations in second-order form in space in [300, 302, 299, 121]. Part of the difficulty in obtaining an energy estimate for second-order variable coefficient systems in these approaches is related to the property that the FD operators now depend on the equations being solved. This also complicates their generalization to arbitrary systems. For example,  deals with systems of the form for generalizations, which include shift-type terms.
In  the projection method and high-order SBP operators approximating second derivatives were used to provide interface boundary conditions for a wave equation (directly in second-order-in-space form) in discontinuous media, while guaranteeing an energy estimate. The domain in this work is a rectangular multi-block domain with piecewise constant coefficients, where the interfaces coincide with the location of the discontinuities. This approach was generalized, using the SAT for the jump discontinuities instead of projection, to variable coefficients and complex geometries in .
The difficulty is not related to FDs: in  a penalty multi-domain method was derived for second-order systems. For the Legendre and constant coefficient case the method guarantees an energy estimate but difficulties are reported in guaranteeing one in the variable coefficient case. Nevertheless, it appears to work well in practice in the variable coefficient case as well. An interesting aspect of the approach of  is that an energy estimate is obtained by applying the penalty in the whole domain (as an analogy we recall the above discussion about the Legendre–Chebyshev penalty method in Section 10.1.3).
A recent generalization, valid both for FDs and – at least Legendre – collocation methods (as discussed, the underlying tool is the same: SBP), to more general penalty couplings, where the penalty terms are not scalar but matrices (i.e., there is coupling between the penalty for different characteristic variables) can be found in .
Energy estimates are in general lost when the different grids are not conforming (different types of domain decompositions are discussed in the next Section 11), and interpolation is needed. This is the case when using overlapping patches with curvilinear coordinates but also mesh refinement with Cartesian, nested boxes (see, for example,  and references therein). A recent promising development has been the introduction of a procedure for systematically constructing interpolation operators preserving the SBP property for arbitrary high-order cases; see . Numerical tests are presented with a 2:1 refinement ratio, where the design convergence rate is observed. It is not clear whether reflections at refinement boundaries such as those reported in  would still need to be addressed or if they would be taken care of by the high-order accuracy.
Finally, we mention some results in numerical relativity concerning absorbing artificial boundaries. In , boundary conditions based on the work of , which are perfectly absorbing for quadrupolar solutions of the flat-wave equation, were numerically implemented via spectral methods, and proposed to be used in a constrained evolution scheme of Einstein’s field equations . For a different method, which provides exact, nonlocal outer boundary conditions for linearized gravitational waves propagating on a Schwarzschild background; see [271, 270, 272]. A numerical implementation of the well-posed harmonic IBVP with Sommerfeld-type boundary conditions given in  was worked out in , where the accuracy of the boundary conditions was also tested.
In , various boundary treatments for the Einstein equations were compared to each other using the test problem of a Schwarzschild black hole perturbed by an outgoing gravitational wave. The solutions from different boundary algorithms were compared to a reference numerical solution obtained by placing the outer boundary at a distance large enough to be causally disconnected from the interior spacetime region where the comparison was performed. The numerical implementation in  was based on the harmonic formulation described in .
In Figure 8, a comparison is shown between (a) simple boundary conditions, which just freeze the incoming characteristic fields to their initial value, and (b) constraint-preserving boundary conditions controlling the complex Weyl scalar at the boundary. The boundary surface is an approximate metric sphere of areal radius , with the mass of the black hole. The left side of the figure demonstrates that case (a) leads to significantly larger reflections than case (b). The difference with the reference solution after the first reflection at the boundary is not only large in case (a), but also it does not decrease with increasing resolution. Furthermore, the violations of the constraints shown in the right side of the figure do not converge away in case (a), indicating that one does not obtain a solution to Einstein’s field equations in the continuum limit. In contrast to this, the difference with the reference solution and the constraint violations both decrease with increasing resolution in case (b).
A similar comparison was performed for alternative boundary conditions, including spatial compactifications and sponge layers. The errors in the gravitational waves were also estimated in , by computing the complex Weyl scalar for the different boundary treatments; see Figure 10 in that reference.
Based on the construction of exact in- and outgoing solutions to the linearized Bianchi identities on a Minkowksi background (cf. Example 32), the reflection coefficient for spurious reflections at a spherical boundary of areal radius , which sets the Weyl scalar to zero, was estimated to be 
For higher-order absorbing boundary conditions, which involve derivatives of the Weyl scalar ; see  and  for their numerical implementation.
Living Rev. Relativity 15, (2012), 9
This work is licensed under a Creative Commons License.