Go to previous page Go up Go to next page

8.7 Flow algorithms

Flow algorithms define a “flow” on 2-surfaces, i.e. they define an evolution of 2-surfaces in some pseudo-time λ, such that the apparent horizon is the λ → ∞ limit of a (any) suitable starting surface. Flow algorithms are different from other apparent horizon finding algorithms (except for zero-finding in spherical symmetry) in that their convergence does not depend on having a good initial guess. In other words, flow algorithms are global algorithms (Section 7.7).

To find the (an) apparent horizon, i.e. an outermost MOTS, the starting surface should be outside the largest possible MOTS in the slice. In practice, it generally suffices to start with a 2-sphere of areal radius substantially greater than 2 mADM.

The global convergence property requires that a flow algorithm always flow from a large starting surface into the apparent horizon. This means that the algorithm gains no particular benefit from already knowing the approximate position of the apparent horizon. In particular, flow algorithms are no faster when “tracking” the apparent horizon (repeatedly finding it at frequent intervals) in a numerical time evolution. (In contrast, in this situation a local apparent horizon finding algorithm can use the most recent previously-found apparent horizon as an initial guess, greatly speeding the algorithm’s convergence55).

Flow algorithms were first proposed for apparent horizon finding by Tod [157Jump To The Next Citation Point]. He initially considered the case of a time-symmetric slice (one where Kij = 0). In this case, a MOTS (and thus an apparent horizon) is a surface of minimal area and may be found by a “mean curvature flow”

i i ∂λx = − κs , (38 )
where xi are the spatial coordinates of a horizon-surface point, si is as before the outward-pointing unit 3-vector normal to the surface, and κ ≡ ∇ksk is the mean curvature of the surface as embedded in the slice. This is a gradient flow for the surface area, and Grayson [79] has proven that if the slice contains a minimum-area surface, this will in fact be the stable λ → ∞ limit of this flow. Unfortunately, this proof is valid only for the time-symmetric case.

For non-time-symmetric slices, Tod [157Jump To The Next Citation Point] proposed generalizing the mean curvature flow to the “expansion flow”

∂λxi = − Θsi. (39 )
There is no theoretical proof that this flow will converge to the (an) apparent horizon, but several lines of argument make this convergence plausible:

Numerical experiments by Bernstein [28Jump To The Next Citation Point], Shoemaker et al. [148Jump To The Next Citation Point149Jump To The Next Citation Point], and Pasch [118Jump To The Next Citation Point] show that in practice the expansion flow (39View Equation) does in fact converge robustly to the apparent horizon.

In the following I discuss a number of important implementation details for, and refinements of, this basic algorithm.

8.7.1 Implicit pseudo-time stepping

Assuming the Strahlkörper surface parameterization (4View Equation), the expansion flow (39View Equation) is a parabolic equation for the horizon shape function h.57 This means that any fully explicit scheme to integrate it (in the pseudo-time λ) must severely restrict its pseudo-time step Δλ for stability, and this restriction grows (quadratically) worse at higher spatial resolutions58. This makes the horizon finding process very slow.

To avoid this restriction, practical implementations of flow algorithms use implicit pseudo-time integration schemes; these can have large pseudo-time steps and still be stable. Because we only care about the λ → ∞ limit, a highly accurate pseudo-time integration is not important; only the accuracy of approximating the spatial derivatives matters. Bernstein [28] used a modified Du Fort–Frankel scheme [64]59 but found some problems with the surface shape gradually developing high-spatial-frequency noise. Pasch [118Jump To The Next Citation Point] reports that an “exponential” integrator (Hochbruck et al. [85]) works well, provided the flow’s Jacobian matrix is computed accurately60. The most common choice is probably that of Shoemaker et al. [148Jump To The Next Citation Point149Jump To The Next Citation Point], who use the iterated Crank–Nicholson (“ICN”) scheme61. They report that this works very well; in particular, they do not report any noise problems.

By refining his finite-element grid (Section 2.3) in a hierarchical manner, Metzger [109Jump To The Next Citation Point] is able to use standard conjugate-gradient elliptic solvers in a multigrid-like fashion62, using each refinement level’s solution as an initial guess for the next higher refinement level’s iterative solution. This greatly speeds the flow integration: Metzger reports that the performance of the overall surface-finding algorithm is “of the same order of magnitude” as that of Thornburg’s AHFinderDirect [156] elliptic-PDE apparent horizon finder (described in Section 8.5.7).

In a more general context than numerical relativity, Osher and Sethian [116] have discussed a general class of numerical algorithms for integrating “fronts propagating with curvature-dependent speed”. These flow a level-set function (Section 2.1) which implicitly locates the actual “front”.

8.7.2 Varying the flow speed

Another important performance optimization of the standard expansion flow (39View Equation) is to replace Θ in the right-hand side by a suitable nonlinear function of Θ, chosen so the surface shrinks faster when it is far from the apparent horizon. For example, Shoemaker et al. [148Jump To The Next Citation Point149Jump To The Next Citation Point] use the flow

[ ( ) ] ∂λxi = − (Θ − c)arctan2 Θ--−-c si (40 ) Θ0
for this purpose, where Θ 0 is the value of Θ on the initial-guess surface, and c (which is gradually decreased towards 0 as the iteration proceeds) is a “goal” value for Θ.

8.7.3 Surface representation and the handling of bifurcations

Since a flow algorithm starts with (topologically) a single large 2-sphere, if there are multiple apparent horizons present the surface must change topology (bifurcate) at some point in the flow. Depending on how the surface is represented, this may be easy or difficult.

Pasch [118] and Shoemaker et al. [148149] use a level-set function approach (Section 2.1). This automatically handles any topology or topology change. However, it has the drawback of requiring the flow to be integrated throughout the entire volume of the slice (or at least in some neighborhood of each surface). This is likely to be much more expensive than only integrating the flow on the surface itself. Shoemaker et al. also generate an explicit Strahlkörper surface representation (Section 2.2), monitoring the surface shape to detect an imminent bifurcation and reparameterizing the shape into 2 separate surfaces if a bifurcation happens.

Metzger [109Jump To The Next Citation Point] uses a finite-element surface representation (Section 2.3), which can represent any topology. However, if the flow bifurcates, then to explicitly represent each apparent horizon the code must detect that the surface self-intersects, which may be expensive.

8.7.4 Gundlach’s “fast flow”

Gundlach [80Jump To The Next Citation Point] introduced the important concept of a “fast flow”. He observed that the subtraction and inversion of the flat-space Laplacian in the Nakamura–Kojima–Oohara spectral integral-iteration algorithm (Section 8.4) is an example of “a standard way of solving nonlinear elliptic problems numerically, namely subtracting a simple linear elliptic operator from the nonlinear one, inverting it by pseudo-spectral algorithms and iterating”. Gundlach then interpreted the Nakamura–Kojima–Oohara algorithm as a type of flow algorithm where each pseudo-time step of the flow corresponds to a single functional-iteration step of the Nakamura–Kojima–Oohara algorithm.

In this framework, Gundlach defines a 2-parameter family of flows interpolating between the Nakamura–Kojima–Oohara algorithm and Tod’s [157] expansion flow (39View Equation),

∂λh = − A (1 − B Δ )− 1ρΘ, (41 )
where A ≥ 0 and B ≥ 0 are parameters, ρ > 0 is a weight functional which depends on h through at most 1st derivatives, Δ is the flat-space Laplacian operator, and −1 (1 − B Δ ) denotes inverting the operator (1 − B Δ ). Gundlach then argues that intermediate “fast flow” members of this family should be useful compromises between the speed of the Nakamura–Kojima–Oohara algorithm and the robustness of Tod’s expansion flow. Based on numerical experiments, Gundlach suggests a particular choice for the weight functional ρ and the parameters A and B. The resulting algorithm updates high-spatial-frequency components of h essentially the same as the Nakamura–Kojima–Oohara algorithm but should reduce low-spatial-frequency error components faster. Gundlach’s algorithm also completely avoids the need for numerically solving Equation (23View Equation) for the a00 coefficient in the Nakamura–Kojima–Oohara algorithm.

Alcubierre’s AHFinder [4Jump To The Next Citation Point] horizon finder includes an implementation of Gundlach’s fast flow algorithm63. AHFinder is implemented as a freely available module (“thorn”) in the External LinkCactus computational toolkit (see Table 2) and has been used by many research groups.

8.7.5 Summary of flow algorithms/codes

Flow algorithms are the only truly global apparent horizon finding algorithms and, as such, can be much more robust than local algorithms. In particular, flow algorithms can guarantee convergence to the outermost MOTS in a slice. Unfortunately, these convergence guarantees hold only for time-symmetric slices.

In the forms which have strong convergence guarantees, flow algorithms tend to be very slow. (Metzger’s algorithm [109Jump To The Next Citation Point] is a notable exception: It is very fast.) There are modifications which can make flow algorithms much faster, but then their convergence is no longer guaranteed. In particular, practical experience has shown that in some binary black hole coalescence simulations (Alcubierre et al. [5], Diener et al. [62]), “fast flow” algorithms (Section 8.7.4) can miss common apparent horizons which are found by other (local) algorithms.

Alcubierre’s apparent horizon finder AHFinder [4] includes a “fast flow” algorithm based on the work of Gundlach [80]63. It is implemented as a freely available module (“thorn”) in the External LinkCactus computational toolkit (see Table 2) and has been used by a number of research groups.

  Go to previous page Go up Go to next page