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. Sharp interfaces can also be present in the initial or boundary conditions. When higher-order numerical methods are used, there are often aphysical undershoots and overshoots in the scalar variables (e.g. passive scalars, species mass fractions or progress variable). These numerical issues are especially prominent when low-dissipation methods are used, since sharp jumps in flow variables are not always coincident with regions of strong variation in the scalar fields: consequently, special detection mechanisms and dissipative fluxes are needed. Most numerical methods diffuse the interface, resulting in artificial mixing and spurious reactions. In this paper, we propose a numerical method that mitigates this issue. We present methods for passive and active scalars, and demonstrate their effectiveness with several examples.