# BSSN Evolution

In this section, we integrate the 3+1 decomposition of Einstein's Equation in the BSSN formalism using SENR. The code is designed to be compact, human-readable, and easily extended.

# Code Features

• Designed to be easily extensible to arbitrary (even dynamical) reference metrics. We have confirmed that the code's evaluation of the BSSN right-hand-sides in Cartesian coordinates agrees to roundoff-error with the BSSN code in the Einstein Toolkit (Kranc/McLachlan), using a set of non-trivial initial data: a highly-spinning black hole in UIUC coordinates (Liu, Etienne, & Shapiro, PRD 80:121503 (2009); arXiv:1001.4077), with nontrivial initial data for the gauge.

• Compact and human readable, so that tensorial expressions can be modified or added on a whim. The code is only about 1500 lines long, partly because expressions are written in a form very similar to LaTeX. For example, 3-connection coefficients

\bar{Gamma}_{ijk} = 1/2 (\bar{\gamma}_{ij,k} + \bar{\gamma}_{ik,j} - \bar{\gamma}_{jk,i})

are very simply written in the code as

F3(i,j,k) { GammabarDDD[i][j][k] = 0.5 * (gammabarDDdD[i][j][k] + gammabarDDdD[i][k][j] - gammabarDDdD[j][k][i]); }

(Note that F3 denotes a nested loop over i, j, and k, which for 3 spatial dimensions ranges from 1 to 3. Notice also the standardized notation D=down=covariant, d=derivative.)

• Arbitrary-order finite difference code. Plug in any positive even number for the finite difference order, and all finite difference stencils are automatically generated (using a tiny automatic code generation utility).

• Supports upwinded/downwinded derivatives on shift advection terms, and the standard moving puncture gauge conditions (with support for advecting shift terms on gauge quantities).

• Formalism for automatically generating efficient code for rescaling gridfunctions a la Eqs. 20—22 of arXiv:1211.6632, as well as computing reference metric connection coefficients (Christoffel symbols) and corresponding covariant ("hatted") derivatives has been developed. We have developed a Mathematica notebook (available on the Download page) that is capable of generating C code that appropriately splits up the BSSN equations for the PIRK formalism, in the context of any user-defined static, uniform coordinate system. Though the current code-generation infrastructure relies on Mathematica, there is no dependency on Mathematica-only functions, and our goal is to migrate away from closed-source, expensive CAS solutions such as Mathematica, to maximize community adoption.

• Agrees to roundoff with the Baumgarte et al spherical BSSN code (arXiv:1211.6632). The initial data represent a strong random perturbation about flat spacetime. This choice eliminates all symmetry and tests every equation. To measure the expected roundoff error, the Baumgarte code is first run with the random initial data, and then rerun with the same initial data, but with a random perturbation in the least significant digit. The resulting evolved differences characterize the expected roundoff error. This is compared with the random initial data evolved by SENR, and is found to agree.

• OMP parallelization allows SENR to divide its effort across multiple CPU threads.

• Kreiss-Oliger dissipation is implemented on the BSSN fields to diffuse unresolved, high frequency modes. The dissipation order is automatically selected to match the finite difference stencil order. The overall strength of the dissipation is allowed to vary over space. For example, in a BH simulation, the dissipation is weak near the puncture and strong outside the horizon.

• Multiple coordinate system options are available, including Cartesian, cylindrical, and spherical, as well as nonuniform log-radial rescaling. The BSSN equations and initial data are automatically transformed into the desired basis. Implementing new distributions is straight forward.

• Supports PIRK2 and RK4 time integration. SENR was initially designed to integrate the BSSN equations using a second order partially implicit Runge-Kutta method, which allows for coordinate singularities, such as those at r = 0 or sin(θ) = 0 in spherical coordinates. It turns out that fourth order Runge-Kutta is also suited for this purpose. The advantage of PIRK2 is that it requires fewer evaluations of the BSSN right-hand-sides, but it is somewhat complicated to implement and suffers from a strong CFL restriction. On the other hand, RK4 requires roughly twice the number of function evaluations, but it has a simple implementation and the CFL factor can be more than doubled relative to PIRK2.