  We introduce NRPyElliptic, an elliptic solver for numerical relativity (NR)
built within the NRPy+ framework. As its first application, NRPyElliptic sets
up conformally flat, binary black hole (BBH) puncture initial data (ID) on a
single numerical domain, similar to the widely used TwoPunctures code. Unlike
TwoPunctures, NRPyElliptic employs a hyperbolic relaxation scheme, whereby
arbitrary elliptic PDEs are trivially transformed into a hyperbolic system of
PDEs. As consumers of NR ID generally already possess expertise in solving
hyperbolic PDEs, they will generally find NRPyElliptic easier to tweak and
extend than other NR elliptic solvers. When evolved forward in (pseudo)time,
the hyperbolic system exponentially reaches a steady state that solves the
elliptic PDEs. Notably NRPyElliptic accelerates the relaxation waves, which
makes it many orders of magnitude faster than the usual constant-wavespeed
approach. While it is still ${\sim}12$x slower than TwoPunctures at setting up
full-3D BBH ID, NRPyElliptic requires only ${\approx}0.3\%$ of the runtime for
a full BBH simulation in the Einstein Toolkit. Future work will focus on
improving performance and generating other types of ID, such as binary neutron
star.
