Go to previous page Go up Go to next page

8.5 Elliptic-PDE algorithms

The basic concept of elliptic-PDE algorithms is simple: We view the apparent horizon equation (16View Equation) as a nonlinear elliptic PDE for the horizon shape function h on the angular-coordinate space and solve this equation by standard finite-differencing techniques38, generally using Newton’s method to solve the resulting set of nonlinear algebraic (finite-difference) equations. Algorithms of this type have been widely used both in axisymmetry and in fully generic slices.

8.5.1 Angular coordinates, grid, and boundary conditions

In more detail, elliptic-PDE algorithms assume that the horizon is a Strahlkörper about some local coordinate origin, and choose an angular coordinate system and a finite-difference grid of Nang points on 2 S in the manner discussed in Section 2.2.

The most common choices are the usual polar-spherical coordinates (θ,φ) and a uniform “latitude/longitude” grid in these coordinates. Since these coordinates are “unwrapped” relative to the actual S2 trial-horizon-surface topology, the horizon shape function h satisfies periodic boundary conditions across the artificial grid boundary at φ = 0 and φ = 2π. The north and south poles θ = 0 and θ = π are trickier, but Huq et al. [89Jump To The Next Citation Point90Jump To The Next Citation Point], Shibata and UryĆ« [147Jump To The Next Citation Point], and Schnetter [132Jump To The Next Citation Point133Jump To The Next Citation Point]39 “reaching across the pole” boundary conditions for these artificial grid boundaries.

Alternatively, Thornburg [156Jump To The Next Citation Point] avoids the z axis coordinate singularity of polar-spherical coordinates by using an “inflated-cube” system of six angular patches to cover 2 S. Here each patch’s nominal grid is surrounded by a “ghost zone” of additional grid points where h is determined by interpolation from the neighboring patches. The interpatch interpolation thus serves to tie the patches together, enforcing the continuity and differentiability of h across patch boundaries. Thornburg reports that this scheme works well but was quite complicated to program.

Overall, the latitude/longitude grid seems to be the superior choice: it works well, is simple to program, and eases interoperation with other software.

8.5.2 Evaluating the expansion Θ

The next step in the algorithm is to evaluate the expansion Θ given by Equation (16View Equation) on the angular grid given a trial horizon surface shape function h on this same grid (6View Equation).

Most researchers compute Θ via 2-dimensional angular finite differencing of Equation (16View Equation) on the trial horizon surface. 2nd order angular finite differencing is most common, but Thornburg [156Jump To The Next Citation Point] uses 4th order angular finite differencing for increased accuracy.

With a (θ,φ ) latitude/longitude grid the Θ (h,∂uh, ∂uvh) function in Equation (16View Equation) is singular on the z axis (at the north and south poles θ = 0 and θ = π) but can be regularized by applying L’Hôpital’s rule. Schnetter [132Jump To The Next Citation Point133Jump To The Next Citation Point] observes that using a Cartesian basis for all tensors greatly aids in this regularization.

Huq et al. [89Jump To The Next Citation Point90Jump To The Next Citation Point] choose, instead, to use a completely different computation technique for Θ, based on 3-dimensional Cartesian finite differencing:

  1. They observe that the scalar field F defined by Equation (7View Equation) can be evaluated at any (3-dimensional) position in the slice by computing the corresponding (r,θ,φ) using the usual flat-space formulas, then interpolating h in the 2-dimensional (θ,φ) surface grid.
  2. Rewrite the apparent horizon condition (15View Equation) in terms of F and its (3-dimensional) Cartesian derivatives,
    Θ ≡ Θ(F, ∂iF,∂ijF;gij,∂kgij,Kij ) = 0. (27)
    Huq et al. [89Jump To The Next Citation Point90Jump To The Next Citation Point] give the Θ (F, ∂iF, ∂ijF ) function explicitly.
  3. For each (latitude/longitude) grid point on the trial horizon surface, define a 3×3 ×3-point local Cartesian grid centered at that point. The spacing of this grid should be such as to allow accurate finite differencing, i.e. in practice it should probably be roughly comparable to that of the underlying numerical-relativity simulation’s grid.
  4. Evaluate F on the local Cartesian grid as described in Step 1 above.
  5. Evaluate the Cartesian derivatives in Equation (27View Equation) by centered 2nd order Cartesian finite differencing of the F values on the local Cartesian grid.

Comparing the different ways of evaluating Θ, 2-dimensional angular finite differencing of Equation (16View Equation) seems to me to be both simpler (easier to program) and likely more efficient than 3-dimensional Cartesian finite differencing of Equation (27View Equation).

8.5.3 Solving the nonlinear elliptic PDE

A variety of algorithms are possible for actually solving the nonlinear elliptic PDE (16View Equation) (or (27View Equation) for the Huq et al. [89Jump To The Next Citation Point90Jump To The Next Citation Point] horizon finder).

The most common choice is to use some variant of Newton’s method. That is, finite differencing Equation (16View Equation) or (27View Equation) (as appropriate) gives a system of Nang nonlinear algebraic equations for the horizon shape function h at the N ang angular grid points; these can be solved by Newton’s method in Nang dimensions. (As explained by Thornburg [153Jump To The Next Citation Point, Section VIII.C], this is usually equivalent to applying the Newton–Kantorovich algorithm [37, Appendix C] to the original nonlinear elliptic PDE (16View Equation) or (27View Equation).)

Newton’s method converges very quickly once the trial horizon surface is sufficiently close to a solution (a MOTS). However, for a less accurate initial guess, Newton’s method may converge very slowly or even fail to converge at all. There is no usable way of determining a priori just how large the radius of convergence of the iteration will be, but in practice 14 to 13 of the horizon radius is often a reasonable estimate40.

Thornburg [153Jump To The Next Citation Point] described the use of various “line search” modifications to Newton’s method to improve its radius and robustness of convergence, and reported that even fairly simple modifications of this sort roughly doubled the radius of convergence.

Schnetter [132Jump To The Next Citation Point133Jump To The Next Citation Point] used the PETSc general-purpose elliptic-solver library [222324] to solve the equations. This offers a wide variety of Newton-like algorithms already implemented in a highly optimized form.

Rather than Newton’s method or one of its variants, Shibata et al. [146Jump To The Next Citation Point147Jump To The Next Citation Point] use a functional-iteration algorithm directly on the nonlinear elliptic PDE (16View Equation). This seems likely to be less efficient than Newton’s method but avoids having to compute and manipulate the Jacobian matrix.

8.5.4 The Jacobian matrix

Newton’s method, and all its variants, require an explicit computation of the Jacobian matrix

∂Θi- Jij = ∂hj , (28 )
where the indices i and j label angular grid points on the horizon surface (or equivalently on S2).

Notice that J includes contributions both from the direct dependence of Θ on h, ∂uh, and ∂uvh, and also from the indirect dependence of Θ on h through the position-dependence of the geometry variables g ij, ∂ g k ij, and K ij (since Θ depends on the geometry variables at the horizon surface position, and this position is determined by h). Thornburg [153Jump To The Next Citation Point] discusses this indirect dependence in detail.

There are two basic ways to compute the Jacobian matrix.

Numerical Perturbation:
 
The simplest way to determine the Jacobian matrix is by “numerical perturbation”, where for each horizon-surface grid point j, h is perturbed by some (small) amount ɛ at the j th grid point (that is, hi → hi + ɛδij), and the expansion Θ is recomputed41. The j th column of the Jacobian matrix (28View Equation) is then estimated as
Θi (h + ɛδij) − Θi(h) Jij ≈ ----------ɛ---------. (29)
Curtis and Reid [53], and Stoer and Bulirsch [150, Section 5.4.3] discuss the optimum choice of ɛ in this algorithm42.

This algorithm is easy to program but somewhat inefficient. It is used by a number of researchers including Schnetter [132Jump To The Next Citation Point133Jump To The Next Citation Point], and Huq et al. [89Jump To The Next Citation Point90Jump To The Next Citation Point].

Symbolic Differentiation:
 
A more efficient, although somewhat more complicated, way to determine the Jacobian matrix is the “symbolic differentiation” algorithm described by Thornburg [153Jump To The Next Citation Point] and also used by Pasch [118Jump To The Next Citation Point], Shibata et al. [146Jump To The Next Citation Point147Jump To The Next Citation Point], and Thornburg [156Jump To The Next Citation Point]. Here the internal structure of the finite differenced Θ (h ) function is used to directly determine the Jacobian matrix elements.

This algorithm is best illustrated by an example which is simpler than the full apparent horizon equation: Consider the flat-space Laplacian in standard (θ,φ) polar-spherical coordinates,

∂θh ∂φφh Δh ≡ ∂θθh + ----- + ---2--. (30) tan θ sin θ
Suppose we discretize this with centered 2nd order finite differences in θ and φ. Then neglecting finite-differencing truncation errors, and temporarily adopting the usual notation for 2-dimensional grid functions, hi,j ≡ h(θ= θi,φ= φj), our discrete approximation to Δh is given by
(Δh ) = hi−1,j-−-2hi,j-+-hi+1,j-+ --1--hi+1,j −-hi−1,j+ --1--hi,j−-1 −-2hi,j +-hi,j+1-. (31) i,j (Δ θ)2 tan θ 2Δ θ sin2 θ (Δ φ )2
The Jacobian of Δh is thus given by
( | --1--- ----1----- |||| (Δθ )2 ± 2 tan θ Δθ if (k,ℓ) = (i±1, j), ||| 1 ∂(Δh )(i,j) { --2-------2- if (k,ℓ) = (i,j±1 ), --∂h------= | sin θ(Δ φ) (32) (k,ℓ) ||| − ---2-- − -----2------ if (k,ℓ) = (i,j), ||| (Δ θ)2 sin2θ (Δ φ)2 |( 0 otherwise.
Thornburg [153Jump To The Next Citation Point] describes how to generalize this to nonlinear differential operators without having to explicitly manipulate the nonlinear finite difference equations.

8.5.5 Solving the linear equations

All the algorithms described in Section 8.5.3 for treating nonlinear elliptic PDEs require solving a sequence of linear systems of Nang equations in Nang unknowns. Nang is typically on the order of a few thousand, and the Jacobian matrices in question are sparse due to the locality of the angular finite differencing (see Section 8.5.4). Thus, for reasonable efficiency, it is essential to use linear solvers that exploit this sparsity. Unfortunately, many such algorithms/codes only handle symmetric positive-definite matrices while, due to the angular boundary conditions43 (see Section 8.5.1), the Jacobian matrices that arise in apparent horizon finding are generally neither of these.

The numerical solution of large sparse linear systems is a whole subfield of numerical analysis. See, for example, Duff, Erisman, and Reid [65], and Saad [130] for extensive discussions44. In practice, a numerical relativist is unlikely to write her own linear solver but, rather, will use an existing subroutine (library).

Kershaw’s [94Jump To The Next Citation Point] ILUCG iterative solver is often used; this is only moderately efficient, but is quite easy to program45. Schnetter [132Jump To The Next Citation Point133Jump To The Next Citation Point] reports good results with an ILU-preconditioned GMRES solver from the PETSc library. Thornburg [156Jump To The Next Citation Point] experimented with both an ILUCG solver and a direct sparse LU decomposition solver (Davis’ UMFPACK library [5758565554]), and found each to be more efficient in some situations; overall, he found the UMFPACK solver to be the best choice.

8.5.6 Sample results

As an example of the results obtained with this type of apparent horizon finder, Figure 10View Image shows the numerically-computed apparent horizons (actually, MOTSs) at two times in a head-on binary black hole collision. (The physical system being simulated here is very similar to that simulated by Matzner et al. [108], a view of whose event horizon is shown in Figure 5View Image.)

View Image

Figure 10: This figure shows the numerically-computed apparent horizons (actually MOTSs) at two times in a head-on binary black hole collision. The black holes are colliding along the z axis. Figure reprinted with permission from [156Jump To The Next Citation Point]. © 2004 by IOP Publishing Ltd.

As another example, Figure 11View Image shows the time dependence of the irreducible masses of apparent horizons found in a (spiraling) binary black hole collision, simulated at several different grid resolutions, as found by both AHFinderDirect and another External LinkCactus apparent horizon finder, AHFinder46. For this evolution, the two apparent horizon finders give irreducible masses which agree to within about 2% for the individual horizons and 0.5% for the common horizon.

View Image

Figure 11: This figure shows the irreducible masses (√area-∕(16π)) of individual and common apparent horizons in a binary black hole collision, as calculated by two different apparent horizon finders in the External LinkCactus toolkit, AHFinder and AHFinderDirect. (AHFinderDirect was also run in simulations at several different resolutions.) Notice that when both apparent horizon finders are run in the same simulation (resolution dx = 0.080), there are only small differences between their results. Figure reprinted with permission from [5Jump To The Next Citation Point]. © 2005 by the American Physical Society.

As a final example, Figure 8View Image shows the numerically-computed event and apparent horizons in the collapse of a rapidly rotating neutron star to a Kerr black hole. (The event horizons were computed using the EHFinder code described in Section 5.3.3.)

8.5.7 Summary of elliptic-PDE algorithms/codes

Elliptic-PDE apparent horizon finders have been developed by many researchers, including Eardley [67]47, Cook [505251], and Thornburg [153Jump To The Next Citation Point] in axisymmetry, and Shibata et al. [146147], Huq et al. [8990], Schnetter [132Jump To The Next Citation Point133Jump To The Next Citation Point], and Thornburg [156Jump To The Next Citation Point] in fully generic slices.

Elliptic-PDE algorithms are (or can be implemented to be) among the fastest horizon finding algorithms. For example, running on a 1.7 GHz processor, Thornburg’s AHFinderDirect [156Jump To The Next Citation Point] averaged 1.7 s per horizon finding, as compared with 61 s for an alternate “fast-flow” apparent horizon finder AHFinder (discussed in more detail in Section 8.7)48. However, achieving maximum performance comes at some cost in implementation effort (e.g. the “symbolic differentiation” Jacobian computation discussed in Section 8.5.4).

Elliptic-PDE algorithms are probably somewhat more robust in their convergence (i.e. they have a slightly larger radius of convergence) than other types of local algorithms, particularly if the “line search” modifications of Newton’s method described by Thornburg [153Jump To The Next Citation Point] are implemented49. Their typical radius of convergence is on the order of 1 3 of the horizon radius, but cases are known where it is much smaller. For example, Schnetter, Herrmann, and Pollney [135Jump To The Next Citation Point] report that (with no “line search” modifications) it is only about 10% for some slices in a binary black hole coalescence simulation.

Schnetter’s TGRapparentHorizon2D [132133Jump To The Next Citation Point] and Thornburg’s AHFinderDirect [156Jump To The Next Citation Point] are both elliptic-PDE apparent horizon finders implemented as freely available modules (“thorns”) in the External LinkCactus computational toolkit (see Table 2). Both work with either the PUGH unigrid driver or the External LinkCarpet mesh-refinement driver for External LinkCactus. TGRapparentHorizon2D is no longer maintained, but AHFinderDirect is actively supported and is now used by many different research groups50.


  Go to previous page Go up Go to next page