Scalar Wave

On this page, we present a number of tests that validate our new, dynamical-background-coordinate techniques for SENR. Tests below are performed in the context of a scalar wave evolution on a flat spacetime background.

Uncompactified TwoPunctures Convergence

Figure 9 demonstrates that our scalar wave evolutions exhibit exponential convergence with increasing finite difference stencil order in our new, uncompactified TwoPunctures coordinate system. Note that we hold the grid resolution and time step fixed. Figure 10 shows the evolution of the scalar wave amplitude and relative error with 8th-order finite-difference stencils. The simulation runs for three light-crossing times with no evidence of a build-up in errors, and the outer boundary is set to the exact solution.

Fig. 9: Convergence of the scalar wave relative error with increasing finite difference stencil order on the uncompactified TwoPunctures grid. (LEFT) The mean of the logarithm of the relative error over the entire 3D grid versus stencil order. The data are fit by an exponential function E ∝ 10ax, where x is the finite difference stencil order and a is a fit parameter. (RIGHT) The mean log relative error along a radial line from the origin to the outer boundary, for different stencil orders. The grid resolution and time step are fixed between runs.
Fig. 10: Evolution of the solution amplitude and relative error of the scalar wave on the uncompactified TwoPunctures grid. (LEFT) The solution amplitude. (RIGHT) The relative error. (TOP) The full uncompactified TwoPunctures grid, extending out to ρ ~ 10. The coordinate distribution behaves like the original TwoPunctures distribution near the origin, and approaches constant radial grid spacing far from the origin. (BOTTOM) A zoom-in view showing one unit area around the right puncture, located at (x, ρ) = (1, 0).

Magnified TwoPunctures Coordinates

One application of SENR will be to model the evolution of compact binaries, composed of neutron stars and black holes. We seek coordinates that are approximately spherical in the vicinity of each compact object, and approach logarithmically distributed spherical coordinates in the far field region. A promising candidate is found in the TwoPunctures coordinates, which are a compactified prolate spheroidal distribution, shown in Fig. 7.

Fig. 7: The original TwoPunctures coordinate distribution.

The TwoPunctures coordinates have foci located at (x, ρ) = (±b, 0), which will correspond to the centers of the compact objects. About each focus, we apply a sort of "magnification transformation" to decrease the local coordinate density. This is critical for increasing our CFL-limited time step to reasonable values. Figure 8 shows the how TwoPunctures coordinates are modified with varying magnification parameters.

Fig. 8: Click the figure to see the animated version.
Our modified TwoPunctures coordinate distribution with varying coordinate density near each focus. The highest density case corresponds to the original TwoPunctures coordinate system (see Fig. 7 above).

Back to Top

Exponential Convergence with Increased Finite Difference Order

For nonuniform coordinate distributions, it can be tricky to define a fair convergence test because a global doubling of coordinate resolution over the entire domain does not necessarily equate to a local doubling. To circumvent this ambiguity, we can perform the analog of a spectral convergence test by keeping constant the coordinate resolution and varying the finite difference stencil order. For a finite difference stencil of order n with step size h, the truncation error scales as Ehn. Thus, we expect exponential convergence with increasing stencil order. Figure 6 shows such a convergence test resulting from integrating the relative error over the same dynamical coordinate grids in the animation below. figure for

Fig. 6: Relative error between exact and numerical solutions versus the finite difference stencil order. Each point represents the relative error integrated over the physical domain. The data are fit by an exponential relationship E ∝ 10ax, where x is the finite difference stencil order and a is a fit parameter.

Back to Top

Logarithmic Dynamical Coordinate Convergence

The physical domain is taken to be flat spacetime. The physical radial coordinate is logarithmically distributed according to the rule that the cell size at point n is related to its neighbor at n - 1 by Δrn = κΔrn-1, where κ > 1 is the physical step size growth factor and Δr0 = Δrmin is given. The angular coordinate φ rotates at a constant velocity. The scalar wave propagates at the speed of light, taking time tc to cross the domain. Figure 5 shows the evolution of a linear combination of eigensolutions with exact outer boundary conditions, as well as a surface plot of the relative error between the numerical and exact solutions. The plot at the bottom of Fig. 5 shows the relative error for two different resolutions, with the high resolution data scaled by a constant corresponding to the convergence factor, demonstrating second-order convergence.

Fig. 5: Click the figure to see the animated version.
An equatorial slice of the -plane from 0 < r < 10. The grid is spaced logarithmically in r-direction, and has a constant rotation frequency about the φ-axis. (LEFT) The numerical solution amplitude; (RIGHT) the relative error; (BOTTOM) the relative error versus radius for two different resolutions, scaled by the convergence order to show second-order convergence. The plus marks on the circular plots sample the locations of the dynamical physical grid points. The time is expressed in terms of tc, the time required for a luminal signal to cross the grid.

Back to Top

Dynamical Coordinates

This extension moves to a 3+1-dimensional evolution. The physical domain is taken to be flat spacetime with somewhat arbitrary coordinate distribution, which is mapped on to a uniform unit sphere grid. The physical coordinates are allowed to be nonuniform and dynamical in time. In this example, the physical radial coordinate is logarithmically spaced from 0.25 < r < 10, and the φ coordinate rotates at a constant angular velocity. The solution is a linear combination of the first three eigenmodes. Figure 4 shows the solution amplitude and absolute error (colors) during the evolution, as well as the dynamical coordinate points (pluses). All finite difference derivatives are evaluated on the static, uniform coordinates (R, Θ, Φ) with ordinary cell-centered finite difference stencils.
Fig. 4: Click the figure to see the animated version.
An equatorial slice of the -plane from 0 < r < 10. The grid is spaced logarithmically in r-direction, and has a constant rotation frequency about the φ-axis. (LEFT) The numerical solution amplitude; (RIGHT) the absolute error. The plus marks indicate the location of the dynamical physical grid points.

Back to Top

Uniform Static Coordinate Convergence

We implement second-order spatial finite difference derivatives and a second-order PIRK integrator. The outer boundary condition is given by the exact eigensolution solution, which in this case is the fundamental mode. The physical domain spans 0 < r < 10, which is remapped to the uniform domain 0 < R < 1 with step size ΔR by a constant rescaling Δr = 10ΔR. Shown in Fig. 3 is the solution amplitude at two resolutions, the corresponding absolute errors, and the measured convergence order N. This is a 1+1-dimensional evolution. The scalar wave propagates at the speed of light, taking time tc to cross the domain.
Fig. 3: Click the figure to see the animated version.
(TOP) Solution amplitude u at two resolutions, where ΔR is the cell spacing on the uniform; (MIDDLE) the absolute error E between the numerical and exact solutions; (BOTTOM) the measured convergence order N.

Back to Top

Benefits of PIRK

The fully explicit Runge-Kutta numerical integration algorithm is widely utilized to approximate the time evolution of systems of differential equations. We implement a second-order scheme to evolve the scalar wave equation in three dimensions on a spherical domain. Coordinate singularities in the spherical Laplace operator at the origin and along the symmetry axis spoil the evolution stability, as shown in Fig. 1. For this example, a known exact solution to the wave equation is used as the initial condition. During the subsequent evolution, the error is evaluated as the difference between the exact and numerical solutions.

Fig. 1: Absolute error of the wave equation solution evolved using the second-order, fully explicit Runge-Kutta algorithm. The gray scale color represents the base-10 logarithm of the absolute difference between the exact and numerical solutions. The plot shows an arbitrary one-dimensional slice of the spherical domain evolving in time. Note the off-scale error contaminating the evolution, originating from the coordinate singularity located at the origin.

The instabilities arising from singularities in the differential operators are remedied by employing a partially implicit Runge-Kutta (PIRK) evolution method. The PIRK algorithm evolves the regular fields explicitly and the singular fields implicitly using a two-step iterative scheme. Figure 2 shows the result of evolving the same initial condition at the same resolution as shown in Fig. 1, but using PIRK to perform the numerical integration.

Fig. 2: Absolute error of the wave equation solution evolved using the second-order, PIRK algorithm. The gray scale color represents the base-10 logarithm of the absolute difference between the exact and numerical solutions. The plot shows an arbitrary one-dimensional slice of the spherical domain evolving in time. Note that the errors remain small for the entire evolution, and in particular that points near the origin are well behaved.

Back to Top