Large eddy simulation is the most viable approach for the accurate simulation of turbulent reacting flows. With the proper combination of high-order, low-dissipation numerical methods, physics-based subgrid-scale models, and boundary conditions it is becoming possible to simulate many combustion flows at relevant conditions. However, non-premixed flows are a particular challenge because the thickness of the fuel/oxidizer interface scales inversely with Reynolds number. Most numerical methods diffuse the interface, resulting in artificial mixing and spurious reactions. Furthermore, when higher-order numerical methods are used, there are often aphysical undershoots and overshoots in the scalar variables (e.g. species mass fractions or progress variable). In this paper, we attempt to construct a numerical method that mitigates this issue.