5.1 Computational boundaries

Boundary conditions are both the most important and the most difficult part of a theoretical treatment of most physical systems. Usually, that’s where all the physics is. And, in computational approaches, that’s usually where all the agony is. Computational boundaries for hyperbolic systems pose special difficulties. Even with an analytic form of the correct physical boundary condition in hand, there are seemingly infinitely more unstable numerical implementations than stable ones. In general, a stable problem places more boundary requirements on the finite difference equations than on the corresponding partial differential equations. Furthermore, the methods of linear stability analysis are often more unwieldy to apply to the boundary than to the interior evolution algorithm.

The von Neumann stability analysis of the interior algorithm linearizes the equations, while assuming a uniform grid with periodic boundary conditions, and checks that the discrete Fourier modes do not grow exponentially. There is an additional stability condition that a boundary introduces into this analysis. Consider the one-dimensional case. The mode ekx, with k real, is not included in the von Neumann analysis for periodic boundary conditions. However, for the half plane problem in the domain x ≤ 0, one can legitimately prescribe such a mode as initial data as long as k > 0 so that it has finite energy. Thus the stability of such boundary modes must be checked. In the case of an additional boundary, e.g. for a problem in the domain − 1 ≤ x ≤ 1, the Godunov–Ryaben’kii theory gives as a necessary condition for stability the combined von Neumann stability of the interior and the stability of the allowed boundary modes [224]. The Kreiss condition [132] strengthens this result by providing a sufficient condition for stability.

The correct physical formulation of any Cauchy problem for an isolated system also involves asymptotic conditions at infinity. These conditions must ensure not only that the total energy and energy loss by radiation are both finite, but they must also ensure the proper 1∕r asymptotic falloff of the radiation fields. However, when treating radiative systems computationally, an outer boundary is often established artificially at some large but finite distance in the wave zone, i.e. many wavelengths from the source. Imposing an appropriate radiation boundary condition at a finite distance is a difficult task even in the case of a simple radiative system evolving on a fixed geometric background. Gustaffson and Kreiss have shown in general that the construction of a nonreflecting boundary condition for an isolated system requires knowledge of the solution in a neighborhood of infinity[131].

When the system is nonlinear and not amenable to an exact solution, a finite outer boundary condition must necessarily introduce spurious physical effects into a Cauchy evolution. The domain of dependence of the initial Cauchy data in the region spanned by the computational grid would shrink in time along ingoing characteristics unless data on a worldtube traced out by the outer grid boundary is included as part of the problem. In order to maintain a causally sensible evolution, this worldtube data must correctly substitute for the missing Cauchy data which would have been supplied if the Cauchy hypersurface had extended to infinity. In a scattering problem, this missing exterior Cauchy data might, for instance, correspond to an incoming pulse initially outside the outer boundary. In a scalar wave problem with field Φ where the initial radiation is confined to a compact region inside the boundary, the missing Cauchy data outside the boundary would be Φ = Φ,t = 0 at the initial time t0. However, the determination of Cauchy data for general relativity is a global elliptic constraint problem so that there is no well defined scheme to confine it to a compact region. Furthermore, even in the scalar field case where Φ = Φ,t = 0 is appropriate Cauchy data outside the boundary at t0, it would still be a non-trivial evolution problem to correctly assign the associated boundary data for t ≥ t0.

It is common practice in computational physics to impose an artificial boundary condition (ABC), such as an outgoing radiation condition, in an attempt to approximate the proper data for the exterior region. This ABC may cause partial reflection of an outgoing wave back into the system [168Jump To The Next Citation Point154Jump To The Next Citation Point138Jump To The Next Citation Point202Jump To The Next Citation Point], which contaminates the accuracy of the interior evolution and the calculation of the radiated waveform. Furthermore, nonlinear waves intrinsically backscatter, which makes it incorrect to try to entirely eliminate incoming radiation from the outer region. The resulting error is of an analytic origin, essentially independent of computational discretization. In general, a systematic reduction of this error can only be achieved by moving the computational boundary to larger and larger radii. This is computationally very expensive, especially for 3-dimensional simulations.

A traditional ABC for the wave equation is the Sommerfeld condition. For a scalar field Φ satisfying the Minkowski space wave equation

ηαβ∂α∂βΦ = S, (35 )
with a smooth source S of compact support emitting outgoing radiation, the exterior retarded field has the form
f(t − r,𝜃,ϕ ) g(t − r,𝜃, ϕ) h(t,r,𝜃,ϕ) Φ = ------------+ ------2-----+ -----3----, (36 ) r r r
where f, g and h and their derivatives are smooth bounded functions. The simplest case is the monopole radiation
Φ = f-(t −-r) (37 ) r
which satisfies (∂ + ∂ )(rΦ ) = 0 t r. This motivates the use of the Sommerfeld condition
1- r (∂t + ∂r)(rΦ)|R = q(t,R, 𝜃,ϕ) (38 )
on a finite boundary r = R.

A homogeneous Sommerfeld condition, i.e.q = 0, is exact only in the spherically symmetric case. The Sommerfeld boundary data q in the general case (36View Equation) falls off as 1∕R3, so that a homogeneous Sommerfeld condition introduces an error, which is small only for large R. As an example, for the dipole solution

( ) f (t − r) f ′(t − r) f(t − r) ΦDipole = ∂z-------- = − ---------+ ----2--- cos𝜃 (39 ) r r r
we have
f (t − r) cos𝜃 q = -------------. (40 ) R3
A homogeneous Sommerfeld condition at r = R would lead to a solution &tidle;ΦDipole containing a reflected ingoing wave. For large R,
&tidle;Φ ∼ Φ + κF-(t +-r-−-2R-)cos-𝜃, (41 ) Dipole Dipole r
where ∂tf(t) = F (t) and the reflection coefficient has asymptotic behavior κ = O (1∕R2 ). More precisely, the Fourier mode
( ) &tidle; eiω(t−r) eiω(t+r−2R) ΦDipole(ω) = ∂z r + κ ω r , (42 )
satisfies the homogeneous boundary condition &tidle; (∂t + ∂r)(rΦDipole(ω)|R = 0 with reflection coefficient
1 1 κω = --2--2------------∼ ---2-2. (43 ) 2ω R + 2iωR − 1 2ω R

Much work has been done on formulating boundary conditions, both exact and approximate, for linear problems in situations that are not spherically symmetric. These boundary conditions are given various names in the literature, e.g., absorbing or non-reflecting. A variety of ABCs have been reported for linear problems. See the articles [107Jump To The Next Citation Point202245Jump To The Next Citation Point209Jump To The Next Citation Point43Jump To The Next Citation Point] for general discussions.

Local ABCs have been extensively applied to linear problems with varying success [168Jump To The Next Citation Point89Jump To The Next Citation Point35Jump To The Next Citation Point244Jump To The Next Citation Point138Jump To The Next Citation Point52155]. Some of these conditions are local approximations to exact integral representations of the solution in the exterior of the computational domain [89Jump To The Next Citation Point], while others are based on approximating the dispersion relation of the so-called one-way wave equations [168244]. Higdon [138] showed that this last approach is essentially equivalent to specifying a finite number of angles of incidence for which the ABCs yield perfect transmission. Local ABCs have also been derived for the linear wave equation by considering the asymptotic behavior of outgoing solutions [35], thus generalizing the Sommerfeld outgoing radiation condition. Although this type of ABC is relatively simple to implement and has a low computational cost, the final accuracy is often limited because the assumptions made about the behavior of the waves are rarely met in practice [107Jump To The Next Citation Point245Jump To The Next Citation Point].

The disadvantages of local ABCs have led some workers to consider exact nonlocal boundary conditions based on integral representations of the infinite domain problem [243107Jump To The Next Citation Point245Jump To The Next Citation Point]. Even for problems where the Green’s function is known and easily computed, such approaches were initially dismissed as impractical [89Jump To The Next Citation Point]; however, the rapid increase in computer power has made it possible to implement exact nonlocal ABCs for the linear wave equation and Maxwell’s equations in 3D [80Jump To The Next Citation Point126]. If properly implemented, this method can yield numerical solutions to a linear problem which converge to the exact infinite domain problem in the continuum limit, while keeping the artificial boundary at a fixed distance. However, due to nonlocality, the computational cost per time step usually grows at a higher power with grid size (𝒪 (N 4) per time step in three dimensions) than in a local approach [107Jump To The Next Citation Point80245Jump To The Next Citation Point].

The extension of ABCs to nonlinear problems is much more difficult. The problem is normally treated by linearizing the region between the outer boundary and infinity, using either local or nonlocal linear ABCs [245Jump To The Next Citation Point209]. The neglect of the nonlinear terms in this region introduces an unavoidable error at the analytic level. But even larger errors are typically introduced in prescribing the outer boundary data. This is a subtle global problem because the correct boundary data must correspond to the continuity of fields and their normal derivatives when extended across the boundary into the linearized exterior. This is a clear requirement for any consistent boundary algorithm, since discontinuities in the field or its derivatives would otherwise act as a spurious sheet source on the boundary, which contaminates both the interior and the exterior evolutions. But the fields and their normal derivatives constitute an overdetermined set of data for the boundary problem. So it is necessary to solve a global linearized problem, not just an exterior one, in order to find the proper data. The designation “exact ABC” is given to an ABC for a nonlinear system whose only error is due to linearization of the exterior. An exact ABC requires the use of global techniques, such as the difference potentials method, to eliminate back reflection at the boundary [245].

There have been only a few applications of ABCs to strongly nonlinear problems [107Jump To The Next Citation Point]. Thompson [240Jump To The Next Citation Point] generalized a previous nonlinear ABC of Hedstrom [137] to treat 1D and 2D problems in gas dynamics. These boundary conditions performed poorly in some situations because of their difficulty in adequately modeling the field outside the computational domain [240107]. Hagstrom and Hariharan [133] have overcome these difficulties in 1D gas dynamics by a clever use of Riemann invariants. They proposed a heuristic generalization of their local ABC to 3D, but this approach has not yet been validated.

In order to reduce the level of approximation at the analytic level, an artificial boundary for a nonlinear problem must be placed sufficiently far from the strong-field region. This sharply increases the computational cost in multi-dimensional simulations [89]. There is no numerical method which converges (as the discretization is refined) to the infinite domain exact solution of a strongly nonlinear wave problem in multi-dimensions, while keeping the artificial boundary fixed. Attempts to use compactified Cauchy hypersurfaces which extend the domain to spatial infinity have failed because the phase of short wavelength radiation varies rapidly in spatial directions [154]. Characteristic evolution avoids this problem by approaching infinity along the phase fronts.

CCM is a strategy that eliminates this nonlinear source of error. In the simplest version of CCM, Cauchy and characteristic evolution algorithms are pasted together in the neighborhood of a worldtube to form a global evolution algorithm. The characteristic algorithm provides an outer boundary condition for the interior Cauchy evolution, while the Cauchy algorithm supplies an inner boundary condition for the characteristic evolution. The matching worldtube provides the geometric framework necessary to relate the two evolutions. The Cauchy foliation slices the worldtube into spherical cross-sections. The characteristic evolution is based upon the outgoing null hypersurfaces emanating from these slices, with the evolution proceeding from one hypersurface to the next by the outward radial march described in Section 3.1). There is no need to truncate spacetime at a finite distance from the source, since compactification of the radial null coordinate used in the characteristic evolution makes it possible to cover the infinite space with a finite computational grid. In this way, the true waveform may be computed up to discretization error by the finite difference algorithm.

  Go to previous page Go up Go to next page