  General relativistic magnetohydrodynamic (GRMHD) simulations represent a
fundamental tool to probe various underlying mechanisms at play during binary
neutron star (BNS) and neutron star (NS) - black hole (BH) mergers.
Contemporary flux-conservative GRMHD codes numerically evolve a set of
conservative equations based on `conserved' variables which then need to be
converted back into the fundamental (`primitive') variables. The corresponding
conservative-to-primitive variable recovery procedure, based on root-finding
algorithms, constitutes one of the core elements of such GRMHD codes. Recently,
a new robust, accurate and efficient recovery scheme called RePrimAnd was
introduced, which has demonstrated the ability to always converge to a unique
solution. The scheme provides fine-grained error policies to handle invalid
states caused by evolution errors, and also provides analytical bounds for the
error of all primitive variables. In this work, we describe the technical
aspects of implementing the RePrimAnd scheme into the GRMHD code Spritz. To
check our implementation as well as to assess the various features of the
scheme, we perform a number of GRMHD tests in three dimensions. Our tests,
which include critical cases such as a NS collapse to a BH as well as the early
evolution (~50 ms) of a Fishbone-Moncrief BH-accrection disk system, show that
RePrimAnd is able to support magnetized, low density environments with
magnetic-to-fluid pressure ratios as high as 10^4, in situations where the
previously used recovery scheme fails.
