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 crosscomparing 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 timestep. 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
3torus. 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 nonflat 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 nonlinear waves, to perturbations of a stationary black
hole and then eventually to highly nonlinear, 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
nonsingular 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
3torus (equivalent to periodic boundary conditions) in order to isolate
problems which are independent of the boundary condition. In Stage II, one
dimension of the 3torus is opened up to form a 2torus 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 wavelike
solution is constructed by carrying out a coordinate transformation on
Minkowski space. The solution can be carried out on a 3torus, 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 nonsingular 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 nonlinear amplitude propagating
either parallel or diagonal to the generators of a 3torus.
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 3torus. 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 23 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),
grqc/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), grqc/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), grqc/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, grqc/0211063.
 Jurke(2003)

T. Jurke,
Class. Quantum Grav. 20,
173 (2003), grqc/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), grqc/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).
