An Integrated Linear Reconstruction for Finite Volume Scheme on Unstructured Grids

Abstract

Linear reconstruction based on local cell-averaged values is the most commonly adopted technique to achieve a second-order accuracy when one uses the finite volume scheme on unstructured grids. For solutions with discontinuities appearing in such as conservation laws, a certain limiter has to be applied to the predicted gradient to prevent numerical oscillations. We propose in this paper a new formulation for linear reconstruction on unstructured grids, which integrates the prediction of the gradient and the limiter together. By solving on each cell a tiny linear programming problem without any parameters, the gradient is directly obtained which satisfies the monotonicity condition. It can be shown that the resulting numerical scheme with our new method fulfils a discrete maximum principle with fair relaxed geometric constraints on grids. Numerical results demonstrate that our method achieves satisfactory numerical accuracy with theoretical guarantee of local discrete maximum principle.

Appendix: All-Inequality Simplex Method

Consider a general LP problem in the all-equality form (20). The Lagrangian of (20) is

$$\begin{aligned} L({{\varvec{x}}},\varvec{\lambda })={{\varvec{c}}}^\top {{\varvec{x}}}+ \varvec{\lambda }^\top ({{\varvec{b}}}-{{\varvec{A}}} {\mathbf {x}}) ={{\varvec{b}}}^\top \varvec{\lambda }+{{\varvec{x}}}^\top ({{\varvec{c}}}-{{\varvec{A}}}^\top \varvec{\lambda }), \end{aligned}$$

and the dual problem of (20) is

$$\begin{aligned} \text {min } {{\varvec{b}}}^\top \varvec{\lambda }\quad \text {s.t. } {{\varvec{A}}}^\top \varvec{\lambda }={{\varvec{c}}}\text { and } \varvec{\lambda }\ge 0. \end{aligned}$$

(39)

By the weak duality theorem, one knows that if ${{\varvec{x}}}$ is feasible for (20) and $\varvec{\lambda }$ is feasible for (39), then

$$\begin{aligned} {{\varvec{c}}}^\top {{\varvec{x}}}\le {{\varvec{b}}}^\top \varvec{\lambda }. \end{aligned}$$

Moreover, ${{\varvec{x}}}$ is optimal for (20) provided that the equality above holds. Now we start from a vertex, and go from vertex to vertex until this equality holds. The iteration takes the form

$$\begin{aligned} {{\varvec{x}}}:={{\varvec{x}}}+\alpha {{\varvec{p}}}, \end{aligned}$$

where the step length $\alpha$ and the searching direction ${{\varvec{p}}}$ need to be determined at each step.

Denote by ${{\varvec{a}}}_i^\top$ the i-th row of ${{\varvec{A}}}$. The active matrix ${{\varvec{M}}}$ consists of d rows of the matrix ${{\varvec{A}}}$, indicating that corresponding constraints at the current searching point ${{\varvec{x}}}$ are active. The active matrix ${{\varvec{M}}}$ needs to be non-singular.

Suppose that at a given step the active Lagrange multipliers $\tilde{\varvec{\lambda }}$ satisfy

$$\begin{aligned} {{\varvec{M}}}^\top \tilde{\varvec{\lambda }}={\varvec{c}} \text { and }\tilde{\varvec{\lambda }}\ge 0, \end{aligned}$$

then ${{\varvec{c}}}^\top {{\varvec{x}}}={{\varvec{b}}}^\top \varvec{\lambda }$ with all the inactive Lagrange multipliers vanishing, and we end up with an optimal point ${{\varvec{x}}}$. Otherwise, we remove the constraint j such that the component $\lambda _j<0$ and determine the searching direction ${\varvec{p}}$ through

$$\begin{aligned} {{\varvec{M}}}{{\mathbf {p}}}=-{{\varvec{e}}}_j, \end{aligned}$$

where ${{\varvec{e}}}_j$ is the j-th coordinate vector.

Assuming that the point in the next step ${{\varvec{x}}}+\alpha _i{{\varvec{p}}}$ satisfies the constraint i, then we have

$$\begin{aligned} {{\varvec{a}}}_i^\top ({{\varvec{x}}}+\alpha _i{{\varvec{p}}})=b_i, \end{aligned}$$

which yields the step length with regard to i to be

$$\begin{aligned} \alpha _i=\dfrac{b_i-{{\varvec{a}}}_i^\top {{\varvec{x}}}}{{{\varvec{a}}}_i^\top {{\varvec{p}}}}. \end{aligned}$$

Since $b_i-{{\varvec{a}}}_i^\top {{\varvec{x}}}\ge 0$, it suffices to consider those constraints satisfying ${{\varvec{a}}}_i^\top {{\varvec{p}}}>0$. The actual step length $\alpha$ is then taken as the minimum of all the estimated step lengths

$$\begin{aligned} \alpha =\min _{i\in D} \{\alpha _i\}. \end{aligned}$$

Finally we choose an index k satisfying $\alpha _k=\alpha$ and update the active matrix by replacing j-th row of ${{\varvec{M}}}$ with k-th row of ${{\varvec{A}}}$. So far we finish one loop. The whole algorithm is illustrated in Algorithm 1. Note that we use the negation of both the searching direction ${{\varvec{p}}}$ and step length $\alpha$.