  One of the main challenges in the numerical modeling of binary neutron-star
(BNS) mergers is long-term simulations of the post-merger remnant over
timescales of the order of seconds. When this modeling includes all the aspects
of complex physics, the computational costs can easily become enormous. To
address this challenge in part, we have developed a novel hybrid approach in
which the solution from a general-relativistic magnetohydrodynamics (GRMHD)
code solving the full set of the Einstein equations in Cartesian coordinates is
coupled with another GRMHD code in which the Einstein equations are solved
under the Conformally Flat Condition (CFC). The latter approximation has a long
history and has been shown to provide an accurate description of compact
objects in non-vacuum spacetimes. An important aspect of the CFC is that the
elliptic equations need to be solved only for a fraction of the steps needed
for the underlying HD/MHD evolution, thus allowing for a gain in computational
efficiency that can be up to a factor of $\sim 6~(230)$ in three-dimensional
(two-dimensional) simulations. We present the basic features of the new code,
the strategies necessary to interface it when importing both two- and
three-dimensional data, and a novel and robust approach to the recovery of the
primitive variables. To validate our new framework, we have carried out code
tests with various coordinate systems and different numbers of spatial
dimensions, involving a variety of astrophysical scenarios, including the
evolution of the post-merger remnant of a BNS merger over a timescale of one
second. \texttt{BHAC+}, can accurately reproduce the evolution of compact
objects in non-vacuum spacetimes and that, when compared with the evolution in
full general relativity, the CFC reproduces accurately both the gravitational
fields and the matter variables at a fraction of the computational costs.
