  We report on a new open-source, user-friendly numerical relativity code
package called SENR/NRPy+. Our code extends previous implementations of the
BSSN reference-metric formulation to a much broader class of curvilinear
coordinate systems, making it ideally-suited to modeling physical
configurations with approximate or exact symmetries. In the context of modeling
black hole dynamics, it is orders of magnitude more efficient than other widely
used open-source numerical relativity codes. NRPy+ provides a Python-based
interface in which equations are written in natural tensorial form and output
at arbitrary finite difference order as highly efficient C code, putting
complex tensorial equations at the scientist's fingertips without the need for
an expensive software license. SENR provides the algorithmic framework that
combines the C codes generated by NRPy+ into a functioning numerical relativity
code. We validate against two other established, state-of-the-art codes, and
achieve excellent agreement. For the first time we demonstrate--in the context
of puncture, trumpet, and dual black hole evolutions--nearly exponential
convergence of constraint violation and gravitational waveform errors to zero
as the order of spatial finite difference derivatives is increased, while
holding the spherical-like coordinate grids fixed at moderate resolution. Such
behavior outside the horizons is remarkable, as numerical errors do not
converge to zero inside horizons, and all points along the polar axis are
coordinate singularities. The formulation addresses such coordinate
singularities via cell-centered grids and a simple change of basis that
analytically regularizes tensor components with respect to the coordinates.
Future plans include extending this formulation to allow dynamical coordinate
grids and bispherical-like distribution of points to efficiently capture
orbiting compact binary dynamics.
