The detection of binary neutron star mergers represents one of the most important astrophysical discoveries of the recent years. Due to the extreme matter and gravity conditions and the rich dynamics developed, it becomes a tremendous challenge to accurately simulate numerically all the scales present during the collision. Here we present how to study such systems by using large eddy simulations with a self-consistent subgrid-scale gradient model, that we generalized to the special relativistic case in a previous work and now extend to the general relativistic case. Adapted from nonrelativistic scenarios, the so-called gradient model allows to capture part of the effects of the hidden dynamics on the resolved scales, by means of a physically-agnostic, mathematically-based Taylor expansion of the nonlinear terms in the conservative evolution equations' fluxes. We assess the validity of this approach in bounding-box simulations of the magnetic Kelvin-Helmholtz instability. Several resolutions and a broad range of scenarios are considered in order to carefully test the performance of the model under three crucial aspects: (i) highly curved backgrounds, (ii) jumps on the fluid density profiles and (iii) strong shocks. The results suggest our extension of the gradient subgrid-scale model to general relativistic magnetohydrodynamics is a promising approach for studying binary neutron stars mergers, and potentially to other relevant astrophysical scenarios.

Performing accurate large eddy simulations in compressible, turbulent magnetohydrodynamics is more challenging than in non-magnetized fluids due to the complex interplay between kinetic, magnetic and internal energy at different scales. Here we extend the sub-grid-scale gradient model, so far used in the momentum and induction equations, to account also for the unresolved scales in the energy evolution equation of a compressible ideal MHD fluid with a generic equation of state. We assess the model by considering box simulations of the turbulence triggered across a shear layer by the Kelvin-Helmholtz instability, testing cases where the small-scale dynamics cannot be fully captured by the resolution considered, such that the efficiency of the simulated dynamo effect depends on the resolution employed. This lack of numerical convergence is actually a currently common issue in several astrophysical problems, where the integral and fastest-growing-instability scales are too far apart to be fully covered numerically. We perform a-priori and a-posteriori tests of the extended gradient model. In the former, we find that, for many different initial conditions and resolutions, the gradient model outperforms other commonly used models in terms of correlation with the residuals coming from the filtering of a high-resolution run. In the second test, we show how a low-resolution run with the gradient model is able to quantitatively reproduce the evolution of the magnetic energy (the integrated value and the spectral distribution) coming from higher-resolution runs. This extension is the first step towards the implementation in relativistic magnetohydrodynamics.