Numerical relativity (NR) simulations are crucial for studying the coalescence of compact binaries. Based on NR data, we produce a model for the mass and spin of the remnant black hole (BH) for the coalescence of black hole-neutron star systems, discussing its crucial role in gravitational wave (GW) modeling and in the parameter estimation of the two signals GW200105 and GW200115. In the context of binary neutron star merger simulations, we perform the first systematic study comparing results obtained with various neutrino treatments, the presence of turbulent viscosity and different grid resolutions. We find that the time of BH formation after merger is heavily affected by grid resolution and turbulent viscosity. An early BH formation limits matter ejection from the accretion disc, as the BH swallows a significant portion of it. Our results indicate that more reliable kilonova light curves are obtained only if the various ejecta components are present. Moreover, robust r-process nucleosynthesis yields require inclusion of both neutrino emission and reabsorption in simulations. Advanced neutrino schemes and turbulent viscosity in simulations resolved beyond current standards appear necessary for reliable astrophysical predictions. To carry out computationally demanding simulations of growing complexity, next-generation NR codes that can efficiently leverage the latest pre-exascale many-core and heterogeneous infrastructures are required. To this end we develop GR-Athena++, a new dynamical spacetime solver built on top of Athena++, that shows high-order convergence properties and excellent parallel scalability up to O(105) cores in full 3D binary black hole (BBH) merger simulations. Finally we present GR-AthenaK, the first performance-portable spacetime solver, obtained by refactoring GR-Athena++ with the Kokkos programming model. We demonstrate the correctness and convergence properties of GR-AthenaK with BBH runs on GPUs. GR-AthenaK shows a speedup ∼50 on one GPU compared to GR-Athena++ on a single CPU core.