10.3 Going further, applications in numerical relativity

In numerical relativity, the projection method for outer boundary conditions has been used in references [102, 103, 421, 278, 225], the penalty FD one for multi-domain boundary conditions in [281Jump To The Next Citation Point, 141Jump To The Next Citation Point, 385Jump To The Next Citation Point, 145Jump To The Next Citation Point, 324Jump To The Next Citation Point, 325Jump To The Next Citation Point, 425, 256Jump To The Next Citation Point, 454, 426], and for spectral methods in – among many others – [130, 289, 90, 168Jump To The Next Citation Point, 306, 450, 402, 131, 149Jump To The Next Citation Point, 290, 97, 91, 31, 381, 150Jump To The Next Citation Point, 384Jump To The Next Citation Point].

[296Jump To The Next Citation Point] 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 7View Image 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 [296Jump To The Next Citation Point] 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.

View Image

Figure 7: Comparison of different numerical boundary approaches for an advection-diffusion equation where the initial data is perturbed, introducing an inconsistency with the boundary condition [296Jump To The Next Citation Point]. The SAT approach, besides guaranteeing time stability for general systems, “washes out” this inconsistency in time. “MODIFIED P” corresponds to a modified projection [226] for which this is solved at the expense of losing an energy estimate. Courtesy: Ken Mattsson. Reprinted with permission from [296]; copyright by Springer.

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, [286Jump To The Next Citation Point, 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, 413Jump To The Next Citation Point] 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, 302Jump To The Next Citation Point, 299Jump To The Next Citation Point, 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, [299] deals with systems of the form

autt = (bux)x, 0 ≤ x ≤ 1, (10.48 )
where a,b : [0,1 ] → ℝ are strictly positive, smooth functions, and difference and dissipative operators of the desired order satisfying SBP approximating (bu ) x x are constructed. See also [302] for generalizations, which include shift-type terms.

In [301] 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 [298].

The difficulty is not related to FDs: in [413Jump To The Next Citation Point] 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 [413] 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 [119].

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, [277Jump To The Next Citation Point] 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 [297]. 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 [39] would still need to be addressed or if they would be taken care of by the high-order accuracy.

10.3.1 Absorbing boundary conditions

Finally, we mention some results in numerical relativity concerning absorbing artificial boundaries. In [314], boundary conditions based on the work of [46], 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 [65]. 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 [267] was worked out in [33], where the accuracy of the boundary conditions was also tested.

In [366Jump To The Next Citation Point], 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 [366Jump To The Next Citation Point] was based on the harmonic formulation described in [286Jump To The Next Citation Point].

In Figure 8View Image, 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 Ψ0 at the boundary. The boundary surface is an approximate metric sphere of areal radius R = 41.9M, with M 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).

View Image

Figure 8: Comparison between boundary conditions in case (a) (solid) and case (b) (dotted); see the body of the text for more details. Four different resolutions are shown: (Nr, L) = (21,8), (31,10), (41,12 ) and (51,14), where N r and L refer to the number of collocation points in the radial and angular directions, where Chebyshev and spherical harmonics are used, respectively. Left panel: the difference Δ 𝒰 between the solution with outer boundary at R = 41.9M and the reference solution. Right panel: the constraint violation 𝒞 (see [366Jump To The Next Citation Point] for precise definitions of these quantities and further details). Courtesy: Oliver Rinne. Reprinted with permission from [366Jump To The Next Citation Point]; copyright by IOP.

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 [366Jump To The Next Citation Point], by computing the complex Weyl scalar Ψ4 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 R, which sets the Weyl scalar Ψ0 to zero, was estimated to be [88]

|γ(kR )| ≈ 3(kR )−4, (10.49 ) 2
for outgoing quadrupolar gravitational radiation with wave number k ≫ R− 1. Figure 9View Image shows the Weyl scalars Ψ0 and Ψ4 computed for the boundary conditions in case (b) and extracted at 1.9M inside the outer boundary. By computing their Fourier transform in time, the overall dependence of the ratio agrees very well with the predicted reflection coefficient.

For higher-order absorbing boundary conditions, which involve derivatives of the Weyl scalar Ψ0; see [369] and [365] for their numerical implementation.

View Image

Figure 9: Comparison of the time Fourier transform of Ψ0 and Ψ4 for two different radii of the outer boundary. The leveling off of Ψ0 for kM ≳ 3 is due to numerical roundoff effects (note the magnitude of Ψ 0 at those frequencies). Courtesy: Oliver Rinne. Reprinted with permission from [366Jump To The Next Citation Point]; copyright by IOP.

  Go to previous page Go up Go to next page