Testing the accuracy and stability of spectral methods in numerical relativity

The accuracy and stability of the Caltech-Cornell pseudospectral code is evaluated using the Kidder, Scheel, and Teukolsky (KST) representation of the Einstein evolution equations. The basic “Mexico City tests” widely adopted by the numerical relativity community are adapted here for codes based on spectral methods. Exponential convergence of the spectral code is established, apparently limited only by numerical roundoff error or by truncation error in the time integration. A general expression for the growth of errors due to finite machine precision is derived, and it is shown that this limit is achieved here for the linear plane-wave test.

Figure 1

Constraints for Minkowski with random noise.—Higher resolutions are expected to have larger constraints because more closely spaced points result in larger derivatives. The constraints do not grow in time.

Figure 2

Error energy for Minkowski with random noise.—The linear increase in time is due to a nonzero average in the random noise added to $K_{ij}$. This average approaches zero as resolution is increased, since there are more points over which to average. The flat line shows the evolution when the average value of $K_{ij}$ is set to zero in the initial data.

Figure 3

Phase error for 1D sinusoidal linear wave.—Phase error at $t=25$ crossing times for various time steps, and several time-stepping algorithms. These tests were all run with three points in the $x$-direction. The dashed line indicates the expected accuracy limit due to roundoff error, cf. Eq. (26).

Figure 4

Constraints for 1D Gaussian linear wave.—Here we see the exponential convergence of the constraints with higher spatial resolution. At late times, the higher resolutions grow sublinearly in time, probably because of accumulated machine roundoff error.

Figure 5

Error energy for 1D Gaussian linear wave.—The solid lines show the error energy at $1/2$ crossing times, with clearly visible exponential convergence at early times. The dashed lines show the error energy at integer crossing times for the same resolutions. The smallness of the error energy at early times demonstrates the low dispersion of the numerical method, as explained in the main text. At later times, the error is dominated by the quadratic growth explained in the text.

Figure 6

Constraints for 2D linear waves.—The solid lines represent the Gaussian wave, while the dashed lines represent the sinusoidal wave. The sinusoid is fully resolved spatially with 3 points. Going to higher resolutions merely introduces spatial errors in the unnecessary basis functions, which leads to an increase in the constraints with resolution.

Figure 7

Error energy for 2D linear waves.—The solid lines represent the Gaussian wave, while dashed lines represent the sinusoidal wave, both observed at $1/2$ crossing times. As in the 1D case (Fig. 5), both sets of evolutions converge to quadratic growth of the error caused by the Hamiltonian constraint, explained in the text.

Figure 8

Constraints for high-amplitude 2D gauge wave.—Dashed lines indicate the unfiltered behavior; solid lines indicate the filtered behavior. Note that, despite an effective loss of resolution, filtering greatly improves the stability of the evolution.

Figure 9

Error energy for high-amplitude 2D gauge wave.—As in Fig. 8, dashed lines are unfiltered, and solid lines are filtered. The growth of the filtered error energy is exponential in time. For the highest resolution the time step was cut in half ($dt=dx/80$) to reduce time-discretization error to the same level as spatial discretization error. The dotted line shows the same evolution with time step $dt=dx/40$, which is dominated by time-discretization error.

Figure 10

Constraints for shifted gauge wave.—Solid lines indicate $A=0.5$, and dashed lines indicate $A=0.1$. For both amplitudes we filter out the top $1/3$ spectral coefficients.

Figure 11

Error energy for shifted gauge wave.—Solid lines indicate $A=0.5$, while dashed lines indicate $A=0.1$. The growth in the $A=0.1$ runs is roughly linear in time, accelerating to quadratic at later times. The dotted line indicates the standard time step ($dt=dx/40$) with 33 points, which is dominated by temporal discretization error, while the blue dashed curve uses $dt=dx/80$.

Figure 12

Constraints for expanding Gowdy spacetime.—At early times, the exponential convergence of spectral methods is clearly visible. Soon, however, the evolutions are dominated by constraints growing roughly as $e^{t/5}$.

Figure 13

Error energy for expanding Gowdy spacetime.—The error energy converges with increased spatial resolution, but $\Vert \delta \mathcal{U}{\Vert }_2/\Vert \mathcal{U}{\Vert }_2$ grows like $e^{t/5}$.

Figure 14

Constraints for collapsing Gowdy spacetime.—Note that the simulation starts at $\tau \sim 9.875$, and proceeds backwards.

Figure 15

Error energy for collapsing Gowdy spacetime.—The simulation starts at $\tau \sim 9.875$, and proceeds backwards.

Figure 16

Comparing constraints for 1D Gaussian linear waves.—Hamiltonian constraint norms (ragged curves) are much smaller than $\Vert \mathcal{C}\Vert$ for this test, so by themselves they are not good diagnostics of constraint violations.