Shock with Confidence: Formal Proofs of Correctness for Hyperbolic Partial Differential Equation Solvers

For the purposes of this paper, we consider homogeneous systems of (potentially non-linear) first-order hyperbolic PDEs, expressed in the generic conservation law form:

where ${\mathbf{U}}$ represents the vector of conserved variables and ${\mathbf{F}\left(\mathbf{U}\right)}$ represents the matrix of fluxes of those variables. We can discretize such a system of PDEs in both space and time, using a uniform grid spacing of ${\Delta x}$ and a time-step of ${\Delta t}$, indexed by subscripts ($i$) and superscripts ($n$):

respectively. Specifically, we employ a finite volume discretization, in which the conserved variable vector ${\mathbf{U}}$ at time ${t^{n}}$ is taken to be constant within each cell ${x_{i}}$ (and therefore piecewise constant across the entire domain), and we adopt the shorthand ${\mathbf{U}_{i}^{n}=\mathbf{U}\left(x_{i},t^{n}\right)}$ to denote this value. The value of the conserved variable vector ${\mathbf{U}_{i}^{n}}$ within each cell can then be evolved in time by means of the explicitly conservative update formula(Toro, 2009):

where ${\mathbf{F}_{i+\frac{1}{2}}}$ and ${\mathbf{F}_{i-\frac{1}{2}}}$ represent inter-cell fluxes, i.e. the interpolated fluxes of the conserved variables between cells ${x_{i}}$ and ${x_{i+1}}$, and between cells ${x_{i-1}}$ and ${x_{i}}$, respectively. Different choices of inter-cell flux function give rise to different finite volume schemes, exhibiting different numerical properties(LeVeque, 2011).

In what follows, we shall focus on four one-dimensional model equation systems in particular, covering all possible combinations of linear vs. non-linear and scalar vs. vector. The first two are purely 1D scalar equations of the form:

namely the linear advection equation with flux ${f\left(u\right)=au}$ (for an arbitrary constant advection velocity ${a\in\mathbb{R}}$), and the non-linear inviscid Burgers’ equation with flux ${f\left(u\right)=\frac{1}{2}u^{2}}$. In each case, the conserved variable $u$ simply represents an arbitrary advected quantity. The second two are 1D vector equation systems of the form:

namely the linear perfectly hyperbolic Maxwell’s equations(Munz et al., 2000), with conserved variable vector:

and flux vector:

and the non-linear isothermal Euler equations, with conserved variable vector:

and flux vector:

For the perfectly hyperbolic Maxwell’s equations, the conserved variables correspond to the $x$, $y$ and $z$ components of the electric and magnetic field vectors ${\mathbf{E}}$ and ${\mathbf{B}}$, as well as the correction potentials ${\phi}$ and ${\psi}$ for cleaning divergence errors in ${\mathbf{E}}$ and ${\mathbf{B}}$, respectively. ${c,\chi,\gamma\in\mathbb{R}}$ are arbitrary constants representing the speed of light, and the propagation speeds of electric and magnetic divergence errors, respectively. For the isothermal Euler equations, the conserved variables correspond to the fluid density ${\rho}$, and the $x$, $y$ and $z$ components of the fluid momentum ${\rho u}$, ${\rho v}$ and ${\rho w}$ (with $u$, $v$ and $w$ representing the $x$, $y$ and $z$ components of the fluid velocity, which are not conserved). ${v_{th}\in\mathbb{R}}$ is an arbitrary constant representing the thermal velocity of the fluid.

One of the simplest choices of inter-cell flux function is the Lax-Friedrichs (finite difference) flux(Lax, 1954)(LeVeque, 1992):

A finite volume solver based on a Lax-Friedrichs inter-cell flux will satisfy the following properties(Breuß, 2004):

The significance of hyperbolicity-preservation is that it is a necessary condition for the solver to remain deterministic, i.e. to prevent discrete solutions from becoming multi-valued(Hartman, 1960). Note that the flux Jacobian ${\mathbf{J}_{\mathbf{F}}}$ being symmetric/Hermitian is a sufficient but not necessary condition for the solver to preserve hyperbolicity. The significance of the CFL stability condition is that, at least for a Lax-Friedrichs solver, CFL stability subsumes all other standard notions of numerical stability (e.g. ${L^{1}}$, ${L^{2}}$ and ${L^{\infty}}$ stability) within a single condition(Courant et al., 1967). Finally, the significance of the local Lipschitz continuity condition, as well as the slightly stronger convexity condition, is that these constraints are sufficient to guarantee that the solver will converge to the physically correct (i.e. thermodynamically consistent) solution in the presence of weak, nonlinear waves, such as shock waves in hydrodynamics(Breuß, 2004).

A more sophisticated choice of inter-cell flux function is the Roe (linearized Riemann problem) flux(Roe, 1981):

In the above, we assume an inter-cell Roe matrix ${\mathbf{A}\left(\mathbf{U}_{i}^{n},\mathbf{U}_{i+1}^{n}\right)}$, which is a linearized approximation to the flux Jacobian ${\mathbf{J}_{\mathbf{F}}}$ that is taken to be constant between cells ${x_{i}}$ and ${x_{i+1}}$, satisfying consistency with the exact Jacobian in the appropriate limit:

such that ${\lambda_{p}}$ are the eigenvalues of ${\mathbf{A}\left(\mathbf{U}_{i}^{n},\mathbf{U}_{i+1}^{n}\right)}$, ${\mathbf{r}_{p}}$ are the (right) eigenvectors of ${\mathbf{A}\left(\mathbf{U}_{i}^{n},\mathbf{U}_{i+1}^{n}\right)}$, and ${\alpha_{p}}$ are the components of the inter-cell jump in the conserved variables ${\mathbf{U}_{i+1}^{n}-\mathbf{U}_{i}^{n}}$, as represented in the ${\mathbf{r}_{p}}$ eigenbasis:

A finite volume solver based on a Roe inter-cell flux will satisfy the following properties(Wesseling, 2001):

Note that Theorem 2.5 regarding Roe solvers is essentially identical in content to Theorem 2.1 regarding Lax-Friedrichs solvers, but with the exact flux Jacobian ${\mathbf{J}_{\mathbf{F}}}$ replaced by the linearized inter-cell approximation ${\mathbf{A}\left(\mathbf{U}_{i}^{n},\mathbf{U}_{i+1}^{n}\right)}$. Note, moreover, that a Lax-Friedrichs solver is guaranteed to conserve the components of ${\mathbf{U}}$ automatically, and does not require any additional conditions (such as those in Theorem 2.6) to hold. Finally, note that the CFL stability and local Lipschitz continuity criteria for Lax-Friedrichs solvers apply to Roe solvers too, i.e. a Roe solver for a given PDE system will not be CFL stable unless the corresponding Lax-Friedrichs solver is also CFL stable, etc.

The Lax-Friedrichs and Roe solvers described above are, naively, only first-order accurate in space. In order to achieve second-order spatial accuracy, it is necessary to replace the piecewise constant approximations of the conserved variable vector ${\mathbf{U}}$ with piecewise linear approximations instead(Van Leer, 1979). At the boundary between cells ${x_{i}}$ and ${x_{i+1}}$, the left- and right-sided extrapolated values of ${\mathbf{U}}$ are given by:

and:

respectively, while at the boundary between cells ${x_{i-1}}$ and ${x_{i}}$, the left- and right-sided extrapolated values of ${\mathbf{U}}$ are given by:

and:

respectively, where ${\mathbf{r}_{i}}$ is a ratio of successive gradients of ${\mathbf{U}}$ (computed componentwise):

The function ${\boldsymbol{\phi}\left(\mathbf{r}_{i}\right)}$ in the above is a flux limiter(LeVeque, 2011), intended to limit the spatial derivatives arising within the reconstruction so as to damp any spurious oscillations that might otherwise appear in the vicinity of steep gradients or true discontinuities. Since all quantities are computed componentwise, it is sufficient to think of the flux limiter as being a scalar function of a scalar ratio of gradients ${\phi\left(r\right)}$.

A second-order spatially accurate scheme can now be constructed by taking the left- and right-extrapolated values ${\mathbf{U}_{i+\frac{1}{2}}^{L}}$ and ${\mathbf{U}_{i+\frac{1}{2}}^{R}}$ and evolving them forward by a half time-step ${\frac{\Delta t}{2}}$, to obtain the evolved states ${\overline{\mathbf{U}_{i+\frac{1}{2}}^{L}}}$ and ${\overline{\mathbf{U}_{i+\frac{1}{2}}^{R}}}$:

and:

respectively. The inter-cell flux ${\mathbf{F}_{i+\frac{1}{2}}}$ for the second-order scheme is now evaluated in the usual way, but with the left and right cell states ${\mathbf{U}_{i}^{n}}$ and ${\mathbf{U}_{i+1}^{n}}$ replaced by the evolved boundary-extrapolated states ${\overline{\mathbf{U}_{i+\frac{1}{2}}^{L}}}$ and ${\overline{\mathbf{U}_{i+\frac{1}{2}}^{R}}}$, respectively, i.e. for the case of Lax-Friedrichs fluxes, one has:

and for the case of Roe fluxes, one has:

for Roe matrix ${\mathbf{A}\left(\overline{\mathbf{U}_{i+\frac{1}{2}}^{L}},\overline{\mathbf{U}_{i+\frac{1}{2}}^{R}}\right)}$.

Flux limiters, and the resulting second-order extrapolations that they enable, satisfy the following properties:

In what follows, we shall focus upon four standard flux limiters in particular, all of which had previously been implemented as part of the Gkeyll code, namely the minmod limiter(Roe, 1986):

the monotonized-centered limiter(Van Leer, 1977):

the superbee limiter(Roe, 1986):

and the van Leer limiter(Van Leer, 1974):

Before we proceed with formally verifying various properties of finite volume schemes in Racket, it is first necessary for us to be able to generate reliable C implementations which certifiably match the symbolic Racket expressions being reasoned about. To this end, we introduce a general data structure for representing hyperbolic PDE systems in Racket, consisting of four lists of symbolic Racket expressions: cons-exprs representing the conserved variable vector, flux-exprs representing the flux vector, max-speed-exprs representing the maximum wave-speed estimates (used for enforcing the CFL stability condition), and parameters representing any additional simulation parameters (such as equation of state variables). For example, the equations representing the density ${\rho}$ and $x$-momentum ${\rho u}$ components of the isothermal Euler equation system are represented in our Racket implementation as follows, assuming a thermal velocity ${v_{th}=1.0}$:

Each of the Racket expressions in these lists may then be converted recursively into a string representing a functionally equivalent C expression, using the function convert-expr:

The base case of convert-expr converts any symbol or number directly into its corresponding string. Any single- or multi-argument function, e.g. (abs arg) or (max arg1 arg2), is converted into its C equivalent, e.g. fabs(arg) or fmax(arg1, arg2), with convert-expr being called on all interior expressions. Finally, any elementary arithmetic operation, e.g. (+ arg1 arg2 ...), has convert-expr mapped over each argument, and the results are interspersed with the corresponding arithmetic symbol, e.g. +, in C. We have generally erred on the side of over-generating parentheses in the resulting C code, in order to guarantee that the order of operations remains identical between the C and Racket versions of mathematical expressions. Conditional expressions in Racket are either converted into the corresponding ternary operators in C:

or are converted directly into if statements. Finally, the strings generated by convert-expr are spliced into a template for a generic finite volume solver using format. Appropriate templates exist for both entirely standalone solvers, and for bespoke solver modules that can be integrated into the larger Gkeyll code.

At its core, our automated theorem-proving framework is based upon a symbolic simplification algorithm, which aims to reduce every symbolic Racket expression to some canonical algebraic form. Although such simplification algorithms are standard in computer algebra, our particular application makes this task non-trivial in two respects. First, we wish for our simplification algorithm to be based on a globally confluent(Robinson and Voronkov, 2001) and strongly normalizing(Baader and Nipkow, 1999) underlying rewriting system, since such rewriting systems exhibit highly desirable correctness and termination properties for the type of algebraic/equational theorem-proving with which we are concerned. This places certain restrictions on the kinds of rewriting rules that we are able to include, since many algebraically correct transformations, such as those corresponding to commutativity of addition or multiplication:

cannot be safely included without risking breaking the strong normalization property of the rewriting system (and hence termination of the theorem-prover). Second, since the variables over which we are reasoning typically represent floating-point numbers, and specifically doubles in C, many standard algebraic rules, such as associativity of addition or multiplication:

cannot safely be assumed to hold (except in certain restricted cases) due to the accumulation of truncation errors. For instance, in floating-point arithmetic(Goldberg, 1991):

To this end, we have tried wherever possible to admit only to those algebraic transformations that are permitted under the IEEE 754 standard for floating-point arithmetic(IEEE Computer Society et al., 2008).

The particular collection of algebraic rewriting rules used within our theorem-prover is somewhat complex and ad hoc, so we shall only summarize the salient elements here. The basic structure consists of a recursively-defined symbolic-simp-rule function of the form:

The first few rules shown above are examples of elementary algebraic properties(Bläsius and Bürckert, 1992), such as the existence of 0 as a (left) additive identity, the existence of 1 as a (left) multiplicative identity, the existence of 0 as a (left) annihilator for multiplication, etc. The next rule is an example of how symbolic-simp-rule is mapped over each term within an elementary arithmetic operation such as (+ arg1 arg2 ...) (similar mappings are performed for operations such as (abs arg) or (max arg1 arg2), with each subexpression being recursively simplified). Finally, if none of the rewriting rules match the expression, then the expression has reached a normal form and is returned verbatim. The symbolic simplifier itself then consists of a single function symbolic-simp which recursively calls symbolic-simp-rule until a fixed point is achieved (i.e. until the expression stops changing):

In order to facilitate the reduction of all symbolic expressions to a canonical form, rewriting rules exist to move all numerical constants or coefficients to the left of non-numerical expressions within sums or products:

to collect “like” terms together within sums and differences (via the distributive property):

and to evaluate any arithmetic expressions or operations involving purely numerical values directly:

etc. A more complicated set of rules and heuristics exist regarding whether and when to expand brackets, factorize subexpressions, and so on. Algebraic rules are also defined for standard mathematical functions such as (sqrt ...) or (abs ...), stating for instance that the square root of the square of a quantity, or the square of the square root of a quantity, is equal to the quantity itself:

or that the absolute value of the negation of a quantity is equal to the quantity itself:

or that the (max ...) and (min ...) functions satisfy (at least in the binary case):

and:

respectively, i.e:

etc.

In order to compute symbolic derivatives of arbitrary Racket expressions, we implement a minimalistic automatic differentiation algorithm, again restricting ourselves to assume only those algebraic properties which hold for arbitrary floating-point numbers. Any Racket expression expr may then be differentiated symbolically with respect to the variable var, by means of the function symbolic-diff:

The base case of symbolic-diff evaluates the derivative of any symbol or numerical constant to 0, unless that symbol matches the variable with respect to which one is differentiating, in which case it evaluates to 1. For any sum of the form (+ arg1 arg2 ...), symbolic-diff is mapped over each term, reflecting the linearity of differentiation (and likewise for differences). For any product of the form (* arg1 arg2 ...), symbolic-diff is applied to each term in the product separately (with all other terms being kept fixed), with the results then being summed together:

reflecting the product rule of differentiation, etc. Certain standard mathematical functions, such as (sqrt ...), also have their derivatives specifically encoded wherever they are well-defined:

With the ability to differentiate arbitrary scalar functions thus in place, the symbolic Jacobian of a list of symbolic Racket expressions exprs, evaluated with respect to a list of symbolic Racket variables vars, may now be computed via a straightforward symbolic-jacobian function:

Likewise for the symbolic gradient of a single symbolic Racket expression expr, with respect to a list of symbolic Racket variables vars, via the symbolic-gradient function:

The symbolic Hessian (symbolic-hessian) of a symbolic Racket expression may therefore be computed as a composition of the two:

In order to prove that a given Lax-Friedrichs solver meets the CFL stability criterion listed in Theorem 2.2, it is sufficient to compute the flux Jacobian ${\mathbf{J}_{\mathbf{F}}}$ by evaluating:

and then to compute its symbolic eigenvalues, which for instance can be done in the case of a 2x2 matrix using:

where we have simply encoded the explicit solutions ${\lambda_{\pm}}$ of the characteristic polynomial for the 2x2 matrix:

i.e:

We can now proceed to apply (abs ...) to each symbolic eigenvalue, and then compare the absolute eigenvalues against the expressions in max-speed-exprs (after mapping symbolic-simp over the two lists of expressions, reducing them to their respective normal forms) in order to confirm that they are indeed identical. Since the time-steps ${\Delta t}$ within the generated C code are computed directly from the elements of max-speed-exprs as:

where ${\left\lvert a\right\rvert}$ is the largest value in max-speed-exprs, this condition is sufficient to guarantee CFL stability under Theorem 2.2 provided that ${0<C_{CFL}\leq 1}$, which is also verified by the theorem-prover during the initialization step:

In order to prove that a given Lax-Friedrichs solver meets the hyperbolicity-preservation criterion listed in Theorem 2.1, one must, furthermore, determine whether each of the symbolic eigenvalues of the flux Jacobian ${\mathbf{J}_{\mathbf{F}}}$ corresponds to a real number. For this purpose, we introduce an is-real function:

The base case of is-real treats any real numerical constant as a real number, and any simulation parameter or conserved variable is also assumed to be real by default. This latter condition is then enforced during the simulation initialization step itself:

We enforce that the set of real numbers is closed under operations like addition, subtraction, and multiplication, e.g:

etc., and also that the set of reals remains closed under division so long as the denominator is non-zero, and under square roots so long as the argument is non-negative; these latter two conditions are enforced via the additional functions is-non-zero and is-non-negative, which will be described momentarily.

Moreover, strict hyperbolicity-preservation of a given Lax-Friedrichs solver can be proven (in the 2x2 case) by proving that the symbolic eigenvalues of the flux Jacobian ${\mathbf{J}_{\mathbf{F}}}$ are not only real but also distinct, using the are-distinct function:

The base case of are-distinct treats any pair of non-equal numerical constants as distinct. Pairs of expressions of the form ${\left\{x,-x\right\}}$ or ${\left\{x+y,x-y\right\}}$ are then also treated as distinct, provided that the expressions $x$ or $y$ are themselves non-zero, respectively:

where the function is-non-zero treats all non-zero numerical constants as non-zero, as its base case:

It also treats all simulation parameters as non-zero, and enforces this condition, using similar logic to what was previously described for is-real. Finally, it treats the product of any two non-zero expressions as non-zero:

Finally, in order to prove that a given Lax-Friedrichs solver meets the local Lipschitz continuity criterion listed in Theorems 2.3 and 2.4, it is sufficient to compute the symbolic Hessian ${\mathbf{H}_{f}}$ by evaluating:

for each flux-expr $f$ in the list flux-exprs, and then to compute its symbolic eigenvalues using symbolic-eigvals2 (at least in the 2x2 case), in order to confirm that each symbolic eigenvalue is non-negative. Non-negativity can be checked using the is-non-negative function:

whose base case treats all non-negative numerical constants as non-negative. Moreover, we enforce that the sum, product or quotient of two non-negative numbers is always non-negative, e.g:

etc. Although we have described these techniques in the context of proving the hyperbolicity-preserving and strict hyperbolicity-preserving properties of Lax-Friedrichs solvers via Theorem 2.1, we note that the same techniques can be used to prove the hyperbolicity-preserving and strict hyperbolicity-preserving properties of Roe solvers also, via Theorem 2.5, by exploiting the fact that a valid symbolic Roe matrix ${\mathbf{A}\left(\mathbf{U}_{i}^{n},\mathbf{U}_{i+1}^{n}\right)}$ can be calculated as an average of the two symbolic flux Jacobians:

and the reality and distinctness of its eigenvalues can therefore be verified in exactly the same way as for ${\mathbf{J}_{\mathbf{F}}}$. The flux conservation/jump continuity condition can be verified purely algebraically, by simply confirming that the two sides of the equation in Theorem 2.6 both reduce to the same canonical form using symbolic-simp.

In order to determine whether a given flux limiter ${\phi\left(r\right)}$ satisfies the symmetry condition presented in Theorem 2.7, it is necessary to perform a ${r\to\frac{1}{r}}$ variable transformation throughout the expression for ${\phi\left(r\right)}$, in order to determine whether this transformed expression is algebraically equivalent to ${\frac{\phi\left(r\right)}{r}}$. For this purpose, we use the recursively-defined variable-transform function:

which applies itself recursively to all subexpressions, and replaces any subexpression it finds which matches var with a new subexpression matching new-var. For the second-order total variation diminishing (TVD) condition presented in Theorem 2.8, we use a strictly stronger set of conditions which imply (but are not equivalent to) the Sweby criteria, namely the limiting conditions that:

assuming that ${\phi\left(r\right)=0}$ for all ${r<0}$, combined with the condition that ${\phi\left(r\right)}$ be a (non-strictly) concave function:

These symbolic limits of the form ${\lim\limits_{x\to x_{0}}\left[f\left(x\right)\right]}$ are evaluated via the evaluate-limit function, which first performs a variable transformation ${x\to x_{0}}$ (using the extended reals ${x\in\mathbb{R}\cup\left\{-\infty,+\infty\right\}}$), and then recursively simplifies until the limiting expression achieves a fixed point:

where evaluate-limit-rule encodes valid algebraic simplification rules over the extended reals (a specialized subset of the algebraic rules encoded by symbolic-simp-rule).

For the two scalar PDEs (i.e. the linear advection and inviscid Burgers’ equations), the hyperbolicity-preservation, CFL stability and local Lipschitz continuity properties of the Lax-Friedrichs solver, and the hyperbolicity-preservation and flux conservation (jump continuity) properties of the Roe solver, can all be proved directly and without further complication, as outlined in Tables 1 and 2. For the two vector PDE systems (i.e. the perfectly hyperbolic Maxwell’s equations and isothermal Euler equations), slightly more work is required. In order to ensure that the theorem-prover never needs to manipulate any higher-dimensional (i.e. beyond 2x2) matrices, we first “factorize” these equation systems into coupled pairs of PDEs, with each pair being dealt with independently. Maxwell’s equations effectively “factorize” into four coupled pairs of linear advection equations, for the ${E^{y}}$ and ${B^{z}}$ components, the ${E^{z}}$ and ${B^{y}}$ components, the ${E^{x}}$ and ${\phi}$ components, and the ${B^{x}}$ and ${\psi}$ components, respectively. The isothermal Euler equations effectively “factorize” into a non-linear pair of equations for the ${\rho}$ and ${\rho u}$ components, and a linear pair of equations for the ${\rho v}$ and ${\rho w}$ components (since the ${\rho v}$ and ${\rho w}$ momentum components can each be regarded as being advected linearly with a speed of $u$ at each time-step). As shown in Tables 3 and 4, for all four coupled systems constituting the perfectly hyperbolic Maxwell’s equations, the hyperbolicity-preservation, strict hyperbolicity-preservation, CFL stability and local Lipschitz continuity properties of the Lax-Friedrichs solver, and the hyperbolicity-preservation, strict hyperbolicity-preservation and flux conservation (jump continuity) properties of the Roe solver, can all be proved unproblematically. However, for the two coupled systems constituting the isothermal Euler equations, several properties cannot immediately be proven: for the Lax-Friedrichs solver, proofs cannot be found for the local Lipschitz continuity property for the ${\rho}$ and ${\rho u}$ components, or for the strict hyperbolicity-preservation property for the ${\rho v}$ and ${\rho w}$ components, while for the Roe solver, proofs cannot be found for any of the properties for the ${\rho}$ and ${\rho u}$ components, or for the strict hyperbolicity-preservation property for the ${\rho v}$ and ${\rho w}$ components. It is worth expanding upon each of these “failure” cases in greater detail.