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, humanreadable, 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 righthandsides in Cartesian coordinates agrees to roundofferror with the BSSN code in the Einstein Toolkit (Kranc/McLachlan), using a set of nontrivial initial data: a highlyspinning 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, 3connection 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
, andk
, which for 3 spatial dimensions ranges from 1 to 3. Notice also the standardized notationD
=down=covariant,d
=derivative.) 
Arbitraryorder 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 userdefined static, uniform coordinate system. Though the current codegeneration infrastructure relies on Mathematica, there is no dependency on Mathematicaonly functions, and our goal is to migrate away from closedsource, 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.

KreissOliger 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 logradial 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 RungeKutta method, which allows for coordinate singularities, such as those at r = 0 or sin(θ) = 0 in spherical coordinates. It turns out that fourth order RungeKutta is also suited for this purpose. The advantage of PIRK2 is that it requires fewer evaluations of the BSSN righthandsides, 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.