One of the most important open issues in the theoretical understanding of binary neutron star collisions is the amplification of magnetic fields after the merger. This amplification occurs first mainly at small scales, in a turbulent regime driven by the Kelvin-Helmholtz instability. This stage is followed by a large-scale ordering driven by winding and the magneto-rotational instability. However, the highest numerical resolution achieved in full GRMHD simulations O(10m) is extremely expensive (tens of millions of CPU hours) and are still far from capturing all the scales at play, possibly of a meter or less. One way to circumvent such limitation is by using large-eddy-simulations, a technique used in in computational fluid dynamics and other fields to model the scales smaller than the size of our numerical grid. By extending these LES to General Relativistic MHD, and using the gradient sub-grid-scale model, we can faithfully reproduce the magnetic amplification up to local values of 10^17 G during the first few milliseconds after the merger, but at a much lower computational cost. We also find convergence of the rms magnetic field of the remnant, at the order of 10^15-10^16 G. These results anticipates new simulations of the outflows and jets associated to these magnetic fields with a much more accurate initial conditions.