Design of standardized tests

Tests for numerical relativity should ideally satisfy a number of important properties. They should be broadly applicable, in the sense that they can be executed within a broad range of formulations. They should produce data which are unambiguous. They should specify free quantities such as the gauge. They should output variables which are independent of the details of the model and can be used for comparison. Finally, the tests and their output should provide some insight into both the characteristics of the code on which they are run and the comparative performance of other codes run on the same problem.

The first of these issues, universality, is not difficult to satisfy within a given class of codes. For instance, we can consider the class of ``'' codes which evolve spacelike surfaces on a finite domain with periodic boundary conditions. For most formulations it should be a straightforward technical problem to express some initial data set in variables appropriate for a given formulation, and to construct the desired output variables at each time step. Even within this class of codes, however, details of implementation might restrict the set of tests which a particular code could run. For instance, codes which are fixed to run on a spherical domain might be difficult to test directly using data specified on a cube with periodic boundaries.

More serious difficulties arise in cross-comparing between classes of codes based upon entirely different formulations. Direct comparison of a finite domain code with a compactified characteristic code or with a Cauchy code based upon hyperboloidal time slices extending to null infinity would be impossible, given the difficulties in defining equivalent initial data. However, for a Cauchy code run on a large spatial domain, one could imagine comparing physically motivated quantities, such as the extracted radiation, with the results of a compactified code.

Of course, each test must specify an appropriate set of output quantities for analysis and comparison. It is common practice to judge the performance of a relativity code by measuring the degree of violation of the constraint equations. A properly working code should satisfy each of its constraints at each time-step. Constraint violation is a clear indication of a problem and usually the first indicator of growing modes which will eventually kill a simulation. On the other hand, a code which accurately preserves the constraints, perhaps by enforcing them during the evolution, might suffer other losses of accuracy which produce a numerically generated spacetime that is unphysical, even though it lies on the constraint hypersurface (for an example, see [Siebel and Hübner(2001)]). Solutions to the constraint equations are not necessarily connected via the evolution equations (e.g., a Schwarzschild spacetime should not evolve into Minkowski space, even though both will satisfy constraints perfectly!). Indeed, mixing constraint equations with the evolution equations can have great effect on the numerical results [Yo et al.(2002)Yo, Baumgarte, and Shapiro]. For these reasons, it is important not to examine the constraints in isolation. Similarly, the length of time a code runs before it ``crashes'' is not an appropriate criterion of quality unless accompanied by some indication of how accurately the code reproduces the intended physics, for example the constant amplitude and phase of a wave in the linear regime propagating inside a box with periodic boundaries.

Ideally, one would like to compare the convergence of variables against exact solutions where the answer is known. These are the most unambiguous tests and the most important for debugging a code. The exact solution provides an explicit choice of lapse and shift which can be adopted by any code through the use of gauge source functions. Furthermore, the exact solution provides explicit boundary data.

However, exact solutions for complex physical problems in relativity are scarce. We need additional criteria by which to judge whether numerically generated spacetimes are ``correct''. For this one could include physical criteria, such as the energy balance between mass and radiation loss (for example, see [Brandt and Seidel(1995),Alcubierre et al.(2001)Alcubierre, Benger, Brügmann, Lanfermann, Nerger, Seidel, and Takahashi]). One could check some expected physical behavior of an evolution, such as the dispersion of a small amplitude Brill wave to evolve toward Minkowski spacetime, or the gravitational collapse of a high amplitude wave [Alcubierre et al.(2000)Alcubierre, Allen, Brügmann, Lanfermann, Seidel, Suen, and Tobias]. Unfortunately, such comparisons between codes can be misleading because of coordinate singularities produced by choice of gauge or by improper choice of boundary condition. For instance, the two spacetimes generated with the same initial data by two codes based upon harmonic time slicing might exhibit qualitatively different regularity properties because of different choices of shift. In such a case, the resulting difference in time slicings even complicates comparisons based upon curvature invariants because of the difficulty in identifying the same spacetime point in two different simulations. The same is true of spacetimes generated with the same initial data and gauge conditions but with different boundary conditions. Brill wave initial data which disperses to Minkowski space with a Sommerfeld boundary condition might collapse to a singularity with a Dirichlet boundary condition.

In the absence of exact solutions, convergence can be measured in the sense of Cauchy convergence. One might consider that a set of tests could be based upon the agreement between a number of independent convergent codes. However, caution needs to be taken with this approach since a common error could lead to acceptance of an incorrect numerical solution which then might systematically bias the modification of later codes to reproduce the accepted solution. Although convergence criteria are an important measure of the quality of any simulation, systematic problems can lead a code to converge to a physically irrelevant solution. Furthermore, agreement between codes is only straightforward when they use the same choice of lapse and shift, which is not always possible. Convergence tests should be accompanied by some expectation of how physical variables should behave. Initial data which have some physical interpretation, such as propagating waves, provide useful tests since they provide some picture of how the fields should behave. Any numerical relativity code should be able to reproduce such expected physical behavior if it is to be trusted when unexpected phenomena are encountered.

Ideally, tests should be ``simple'' as well as physical. For a simple enough test, some properties of a correct solution can be developed even if an exact solution is not known. Particular advantages and disadvantages of a code can be more readily isolated using simple tests for which the behavior of individual variables is understood. Of special importance, a simple test can be more readily implemented by a broad set of codes.

A characteristic feature of asymptotically flat systems is that energy gets lost through radiation. Correspondingly, the numerical treatment of such systems typically leads to an initial boundary value problem. However, in order to create a test suite which allows separation of problems associated with the boundary treatment from other problems, we believe it is fruitful to also consider tests with periodic boundaries. The use of periodic boundary conditions in testing numerical codes is a standard technique in computational science. In the context of general relativity, however, such tests have to be set up and interpreted with great care, since periodic boundaries create the cosmological context of a compact spatial manifold with the topology of a 3-torus. The resulting physics is qualitatively very different from asymptotically flat solutions. An example is the nonlinear stability of Minkowski spacetime, on which - at least implicitly - many numerical stability tests rely: For weak asymptotically flat data, it has been proved that perturbations decay to Minkowski spacetime [Christodoulou and Klainerman(1993),Friedrich(1986)]. This result does not hold in the context of periodic boundaries! For two Killing vectors, it has been shown [Andreasson et al.(2002)Andreasson, Rendall, and Weaver] that all initial data which are not flat fall into two classes which are related to each other by time reversal. Making the standard cosmological choice of time direction, all non-flat data have a crushing singularity in the past and exist globally in a certain sense in the future. In particular, the spacetime can be covered by constant mean curvature hypersurfaces whose mean curvature goes to zero in the future [Andreasson et al.(2002)Andreasson, Rendall, and Weaver,Jurke(2003)]. In the simulation of such a spacetime, although the initial data analytically implies expansion, discretization error or constraint violation can drive the numerically perturbed spacetime toward a singularity in the future.

When (in the present conventions) everywhere on a time slice, the pitfalls are easy to recognize. In that case, the singularity theorems (see e.g. Refs. [Hawking and Ellis(1973)] or [Wald(1984)]) imply that a singularity must form in finite proper time - even if the initial data are arbitrarily close to locally Minkowski data. Apart from such physical singularities, focusing effects can also create gauge singularities which could also lead to a code crash. Even in the case where the spacetime undergoes infinite expansion, the code can crash due to choice of evolution variables that become infinite.

These observations introduce various complications for the setup and interpretation of numerical tests with periodic boundaries. For ``linearized waves'' with initial data close to Minkowski data, constraint violation or other source of error might lead to either expansion or collapse, with different codes exhibiting different behavior. For Minkowski data in nonstandard coordinates the situation is similar: Individual runs should be expected to show instabilities but with grid refinement the runs should show convergence to Minkowski spacetime. For situations with on a whole time slice all codes should exhibit eventual collapse but in other cases the qualitative behavior may be code dependent. In order to interpret the results of a simulation it is thus vital to probe the underlying dynamics in terms of expansion versus collapse - in particular since many gauge conditions for the lapse involve . For this purpose, it is important to monitor such variables as , the eigenvalues of , and the total volume or proper time.

Since any given test will satisfy only a few of the desired criteria, some balance between redundancy and economy is needed. We can envisage a hierarchy of tests, starting from evolving flat space under various gauge assumptions, to linearized and then non-linear waves, to perturbations of a stationary black hole and then eventually to highly non-linear, dynamic black holes. Each successive test should introduce a new feature for which code performance can be isolated. The major goal of numerical relativity is the simulation of binary black holes. This requires special techniques, such as singularity excision, which by themselves are an extreme test for any code and can obscure the precise source of an instability arising in such a strong field regime. A standardized test suite should lead up to the binary problem through models of static, moving and perturbed black holes.

It is not our intention in this paper to present the specifics of a complete test suite. We concentrate here on an initial round of simple tests which serve to highlight certain important characteristics of the codes represented at the Mexico workshop and which can be readily performed with most codes. Even in trial implementations of these simple tests we have found complications that warn us that continuous feedback between design and experiment is absolutely necessary in developing a full test suite. The four tests in this initial round are: (i) the robust stability test, (ii) the gauge wave test, (iii) the linearized wave test and (iv) the Gowdy wave test.

Robust stability [Szilagyi et al.(2000)Szilagyi, Gómez, Bishop, and Winicour] is a particular good example of a testbed satisfying the above criteria. Random constraint violating initial data in the linearized regime is used to simulate unavoidable machine error. It can be universally applied, since any code can run perturbations of Minkowski space, and it is very efficient at revealing unstable modes, since at early times they are not hidden beneath a larger signal. Even in a convergence test based upon a non-singular solution, smoother truncation error dominates, at least until very late times, unless the resolution is very high. In stage I of the testbed, which we include in our first round of tests, evolution is carried out on a 3-torus (equivalent to periodic boundary conditions) in order to isolate problems which are independent of the boundary condition. In Stage II, one dimension of the 3-torus is opened up to form a 2-torus with plane boundaries, where random boundary data is applied at all times. This tests for stability in the presence of a smooth boundary. Boundary conditions which depend upon the direction of the outer normal, such as Neumann or Sommerfeld conditions, are best tested first with smooth boundaries in order to isolate problems. In Stage III, random data is applied at all faces of a cubic boundary, the common choice of outer boundary in simulating an isolated system. This tests for stability in the presence of edges and corners on the boundary. The test can be extended to Stage IV in which a spherical boundary is cut out of a Cartesian grid. These tests provide an efficient way to cull out unstable algorithms as a precursor to more time consuming convergence tests. The performance can be monitored by outputting the value of the Hamiltonian constraint.

The gauge wave testbed also very aptly fits our criteria. An exact wave-like solution is constructed by carrying out a coordinate transformation on Minkowski space. The solution can be carried out on a 3-torus, by matching the wavelength to the size of the evolution domain, or it can be carried out in the presence of boundaries. Since the gauge choice and boundary data are explicitly known, it is easy to carry out code comparisons. The evolution can be performed in the weak field regime or in the extreme strong field regime which borders on creating a coordinate singularity. Knowledge of a non-singular exact solution allows any instability to be attributed to code performance. Long, high resolution evolutions can be performed. Convergence criteria for the numerical solution can be easily incorporated into the test. In our first round of tests, we simulate a gauge wave of moderately non-linear amplitude propagating either parallel or diagonal to the generators of a 3-torus.

The linearized wave testbed uses a solution to the linearized Einstein equations and is complimentary to the first two. Run in the linearized regime, it provides an effectively exact solution of physical importance which can be used to check the amplitude and phase of a gravitational wave as it propagates on the 3-torus. If this test were run with a higher amplitude, so that nonlinear effects are not lost in machine noise, the constraint violation in the initial linearized data would introduce a complication. The tendency toward gravitational collapse could also amplify numerical error or the effect of a bad choice of gauge, e.g. the focusing effect inherent in Gaussian coordinates. A maximal slicing gauge cannot be used to avoid such problems, because, for any nontrivial solution, it would be inconsistent with periodic boundary conditions. While in the asymptotically flat context it is typical to use maximal initial hypersurfaces, it is known that a maximal Cauchy surface of topology has to correspond to locally Minkowski data (see e.g. Ref. [Schoen and Yau(1979)] and pages 2-3 of Ref. [Anderson(2001)]). This result is connected to the fact that does not admit metrics of positive Yamabe number, with the Yamabe number of flat being zero. In order to avoid these complications, our nonlinear wave tests, which complete the present round of testbeds, are based upon polarized Gowdy spacetimes. These spacetimes provide a family of exact solutions describing an expanding universe containing plane polarized gravitational waves [Gowdy(1971)] and a clear physical picture to test against the results of a simulation. We will carry out this test both in the expanding and collapsing time directions, which yields physically very different situations with potentially different mechanisms to trigger instabilities.

Except for the robust stability test, the remaining testbeds are based upon solutions with translational symmetry along one or two coordinate axes. By using the minimal required number of grid points along the symmetry axes, this allows the tests to be run with several grid sizes without exorbitant computational expense. Such tests could also be run with either 1D or 2D codes. In order to check that a 3D code is not taking undue advantage of the symmetry, the initial data can be superimposed with random noise as in the robust stability test. For purposes of economy we do not officially include this as part of the testbeds, but it is a useful practice which can completely alter the performance of a code.

Bibliography

Siebel and Hübner(2001)

F. Siebel and P. Hübner, Phys. Rev. D 64, 024021 (2001), gr-qc/0104079.

Yo et al.(2002)Yo, Baumgarte, and Shapiro

H.-J. Yo, T. Baumgarte, and S. Shapiro, Phys. Rev. D 66, 084026 (2002).

Brandt and Seidel(1995)

S. Brandt and E. Seidel, Phys. Rev. D 52, 870 (1995).

Alcubierre et al.(2001)Alcubierre, Benger, Brügmann, Lanfermann, Nerger, Seidel, and Takahashi

M. Alcubierre, W. Benger, B. Brügmann, G. Lanfermann, L. Nerger, E. Seidel, and R. Takahashi, Phys. Rev. Lett. 87, 271103 (2001), gr-qc/0012079.

Alcubierre et al.(2000)Alcubierre, Allen, Brügmann, Lanfermann, Seidel, Suen, and Tobias

M. Alcubierre, G. Allen, B. Brügmann, G. Lanfermann, E. Seidel, W.-M. Suen, and M. Tobias, Phys. Rev. D 61, 041501 (R) (2000), gr-qc/9904013.

Christodoulou and Klainerman(1993)

D. Christodoulou and S. Klainerman, The Global Nonlinear Stability of the Minkowski Space (Princeton University Press, Princeton, 1993).

Friedrich(1986)

H. Friedrich, Comm. Math. Phys. 107, 587 (1986).

Andreasson et al.(2002)Andreasson, Rendall, and Weaver

H. Andreasson, A. Rendall, and M. Weaver (2002), submitted, gr-qc/0211063.

Jurke(2003)

T. Jurke, Class. Quantum Grav. 20, 173 (2003), gr-qc/0210022.

Hawking and Ellis(1973)

S. W. Hawking and G. F. R. Ellis, The Large Scale Structure of Spacetime (Cambridge University Press, Cambridge, England, 1973).

Wald(1984)

R. M. Wald, General Relativity (The University of Chicago Press, Chicago, 1984).

Szilagyi et al.(2000)Szilagyi, Gómez, Bishop, and Winicour

B. Szilagyi, R. Gómez, N. T. Bishop, and J. Winicour, Phys. Rev. D 62 (2000), gr-qc/9912030.

Schoen and Yau(1979)

R. Schoen and S. T. Yau, Manuscr. Math. 28, 159 (1979).

Anderson(2001)

M. T. Anderson, Commun. Math. Phys. 222 222, 533 (2001).

Gowdy(1971)

R. H. Gowdy, Phys. Rev. Lett. 27, 826 (1971).