Physical quantities like pressure, density, energy, temperature, and concentration should always be positive, but numerical methods sometimes compute negative values during the solution process. This is not okay because the equations will compute complex or infinite values (typically crashing the code). Which numerical methods can be used to guarantee that these quantities remain positive? Which of these methods is most efficient?
Which methods can ensure that physical quantities remain positive throughout a PDE simulation?
pde;cfd;hyperbolic pde
The most common method is to reset negative values to some small, positive number.  Of course, this is not a mathematically sound solution.  A better general approach that may work and is easy, is to reduce the size of your time step.Negative values often arise in the solution of hyperbolic PDEs, because the appearance of shocks can lead to oscillations, which will tend to create negative values if there are near-vacuum states near the shock.  Using a total variation diminishing (TVD) or other non-oscillatory (ENO, WENO) method can reduce this tendency.  Those methods are based on using nonlinear limiters to compute derivatives of the solution.  However, you may still get negative values for several reasons:If you use the method of lines and apply a high-order time integrator.  Most TVD schemes are provably TVD only in the semi-discrete sense or with Euler's method.  For higher order time integration, you should use a strong stability preserving (SSP) time discretization; these schemes are also known as contractive or monotonicity preserving.  There is a recent book on the subject by Sigal Gottlieb, Chi-Wang Shu, and myself.If you don't use local characteristic decomposition for systems of equations, your solution will not be TVD (TVD schemes only possess that property for scalar problems).  So it's best to reconstruct/interpolate in characteristic variables.If you have a nonlinear system, negative values can arise even if you use local characteristic decomposition.  For instance, any linearized Riemann solver (such as a Roe solver) for the shallow water equations or the Euler equations can be shown to generate negative values in sufficiently challenging conditions.  A solution is to use an HLL solver (or a variant of HLL); some of those are provably positive.TVD schemes are only second order; higher order non-oscillatory schemes like WENO do not strictly satisfy TVD or maximum principles.  But a new modification of those high-order schemes does; it is developed in several recent papers by Xiangxiong Zhang (a student of Chi-Wang Shu).There are, of course, many other specialized approaches for particular equations, such as in David George's GeoClaw code, which uses a Riemann solver with extra non-physical waves to enforce positivity.