  Numerical relativity is central to the investigation of astrophysical sources
in the dynamical and strong-field gravity regime, such as binary black hole and
neutron star coalescences. Current challenges set by gravitational-wave and
multi-messenger astronomy call for highly performant and scalable codes on
modern massively-parallel architectures. We present GR-Athena++, a
general-relativistic, high-order, vertex-centered solver that extends the
oct-tree, adaptive mesh refinement capabilities of the astrophysical
(radiation) magnetohydrodynamics code Athena++. To simulate dynamical
space-times GR-Athena++ uses the Z4c evolution scheme of numerical relativity
coupled to the moving puncture gauge. We demonstrate stable and accurate binary
black hole merger evolutions via extensive convergence testing, cross-code
validation, and verification against state-of-the-art effective-one-body
waveforms. GR-Athena++ leverages the task-based parallelism paradigm of
Athena++ to achieve excellent scalability. We measure strong scaling
efficiencies above $95\%$ for up to $\sim 1.2\times10^4$ CPUs and excellent
weak scaling is shown up to $\sim 10^5$ CPUs in a production binary black hole
setup with adaptive mesh refinement. GR-Athena++ thus allows for the robust
simulation of compact binary coalescences and offers a viable path towards
numerical relativity at exascale.
