The tremendous challenge of comparing our theoretical models with the gravitational-wave observations in the new era of multimessenger astronomy requires accurate and fast numerical simulations of complicated physical systems described by the Einstein and the matter equations. These requirements can only be satisfied if the simulations can be parallelized efficiently on a large number of processors and advanced computational strategies are adopted. To achieve this goal we have developed Simflowny, an open platform for scientific dynamical models which automatically generates parallel code for different simulation frameworks, allowing the use of HPC infrastructures to non-specialist scientists. One of these frameworks is SAMRAI, a mature patch-based structured adaptive mesh refinement infrastructure, capable of reaching exascale in some specific problems.Here we present the numerical techniques that we have implemented on this framework by using Simflowny in order to perform fast, efficient, accurate and highly-scalable simulations. These techniques involve high-order schemes for smooth and non-smooth solutions, Adaptive Mesh Refinement with arbitrary resolution ratios and an optimal strategy for the sub-cycling in time. We validate the automatically generated codes for the SAMRAI infrastructure with some simple test examples (i.e., wave equation and Newtonian MHD) and finally with the Einstein equations.

In several relativistic astrophysics scenarios, the understanding of the rich magnetohydrodynamics is hampered by the limitations set by the achievable numerical resolution. In these cases, it is a tremendous challenge to accurately simulate numerically all the relevant scales. 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 us 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. One of the main applications is the binary neutron star mergers, where the collision excites a nontrivial amplification at small spatial scales. Motivated by this scenario, 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 that our extension of the gradient subgrid-scale model to general relativistic magnetohydrodynamics is a promising approach for studying binary neutron star mergers and other relevant astrophysical scenarios.

The detection of binary neutron star mergers represents one of the most important and complex astrophysical discoveries of the recent years. One of the unclear aspects of the problem is the turbulent magnetic field amplification, initially triggered by the Kelvin-Helmholtz instability at much smaller scales than any reachable numerical resolution nowadays. Here we present numerical simulations of the first 10 milliseconds of a binary neutron star merger. First, we confirm in detail how the simulated amplification depends on the numerical resolution and is distributed on a broad range of scales, as expected from turbulent magnetohydrodynamics theory. We find that an initial large-scale magnetic field of 10 11 G inside each star is amplified in the remnant to root-mean-square values above 10 16 G within the first 5 milliseconds for our highest-resolution run. Then, we run large eddy simulations, exploring the performance of the subgrid-scale gradient model, already tested successfully in previous turbulent box simulations. We show that the addition of this model is especially important in the induction equation, since it leads to an amplification of the magnetic field comparable to a higher-resolution run, but with a greatly reduced computational cost. In the first 10 milliseconds, there is no clear hint for an ordered, large-scale magnetic field, which should indeed occur in longer timescales through magnetic winding and the magnetorotational instability.

In the interior of neutron stars, the induction equation regulates the long-term evolution of the magnetic fields by means of resistivity, Hall dynamics and ambipolar diffusion. Despite the apparent simplicity and compactness of the equation, the dynamics it describes is not trivial and its understanding relies on accurate numerical simulations. While a few works in 2D have reached a mature stage and a consensus on the general dynamics at least for some simple initial data, only few attempts have been performed in 3D, due to the computational costs and the need for a proper numerical treatment of the intrinsic non-linearity of the equation.Here, we carefully analyze the general induction equation, studying its characteristic structure, and we present a new Cartesian 3D code, generated by the user-friendly, publicly available Simflowny platform. The code uses high-order numerical schemes for the time and spatial discretization, and relies on the highly-scalable SAMRAI architecture for the adaptive mesh refinement. We present the application of the code to several benchmark tests, showing the high order of convergence and accuracy achieved and the capabilities in terms of magnetic shock resolution and three-dimensionality. This paper paves the way for the applications to a realistic, 3D long-term evolution of neutron stars interior and, possibly, of other astrophysical sources.