Which methods can ensure that physical quantities remain positive throughout a PDE simulation?

Question

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?

Answer 1

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.

Answer 2

Assuming we are solving hyperbolic equations without any source terms and assuming we provide physical initial conditions, making sure the numerical scheme we use is Total Variation Diminishing is a good way of ensuring the "physicality" of the computed solution. Since a TVD scheme preserves monotonicity, no new minima or maxima will be created and the solution will remain bounded by the initial values that we hopefully set correctly. Of course the issue is that TVD schemes are not the most obvious ones. Among linear schemes, only first order schemes are TVD (Godunov 1954). So since the 50's, a variety of non-linear TVD schemes have been developed to combine high-accuracy and monotonicity for the solution of hyperbolic equations.

For my applications, solving Navier-Stokes equations with large pressure/density gradients, we use a hybrid MUSCL-central scheme to capture the large gradients/discontinuities and retain good accuracy away from them. The first MUSCL scheme (MUSCL stands for Monotone Upstream-centered Schemes for Conservation Laws) was devised by Van Leer in 1979.

If you want to know more about this subject, please consult the works of Harten, Van Leer, Lax, Sod and Toro.

Answer 3

The above answers apply to time-dependent problems, but you could also demand positivity in a simple elliptic equation. In this case, you could formulate it as a variational inequality, giving bounds for the variables.

In PETSc, there are two VI solvers. One uses a reduced-space method, where variables in active constraints are removed from the system to be solved. The other uses a semi-smooth Newton method.

Answer 4

The answer to this question is that the system matrix $A$ of the linear equation system $$ Au = b $$ that we obtain from a discretization scheme must be a M-matrix. This means in particular that $A$ is inverse monotone and this means that inverse matrix $A^{-1}$ is a monotone matrix.

A monotone matrix $B\in\mathbb{R}^{n \times n}$ is a matrix, which has solely nonnegative entries and the symbolic expression for this is $B\geq0$. Furthermore, for a monotone matrix $B$ holds: $$ (B\geq 0) \qquad\Leftrightarrow\qquad (u\leq v ~\Rightarrow~ Bu\leq Bv, ~~\forall u,v \in \mathbb{R}^n ) $$

This condition applied to the inverse monotone system matrix $A$ means that it holds for the above linear equation system $$ 0 \leq b ~\Rightarrow~ 0=A^{-1}0 \leq A^{-1}b = u $$ Hence, the solution remains nonnegative, in case the right hand side $b$ of our system is nonnegative. Usually it is $b\geq 0$ in case of nonnegative initial values and nonnegative boundary data.

Commonly, discretization schemes that lead to a M-matrix are called monotone schemes and this are those schemes, which preserve non negativity.