Integrated Linear Reconstruction for Finite Volume Scheme on Arbitrary Unstructured Grids

TL;DR

This paper introduces an improved integrated linear reconstruction method for finite volume schemes on arbitrary unstructured grids, removing previous geometric restrictions and ensuring positivity-preserving properties.

Contribution

The paper presents a new convex quadratic programming-based reconstruction approach that satisfies the local maximum principle on arbitrary unstructured meshes without geometric constraints.

Findings

Ensures local maximum principle on arbitrary unstructured grids.

Reconstruction is parameter-free and efficient.

Finite volume scheme remains positivity-preserving for Euler equations.

Abstract

In [L. Chen and R. Li, Journal of Scientific Computing, Vol. 68, pp. 1172--1197, (2016)], an integrated linear reconstruction was proposed for finite volume methods on unstructured grids. However, the geometric hypothesis of the mesh to enforce a local maximum principle is too restrictive to be satisfied by, for example, locally refined meshes or distorted meshes generated by arbitrary Lagrangian-Eulerian methods in practical applications. In this paper, we propose an improved integrated linear reconstruction approach to get rid of the geometric hypothesis. The resulting optimization problem is a convex quadratic programming problem, and hence can be solved efficiently by classical active-set methods. The features of the improved integrated linear reconstruction include that i). the local maximum principle is fulfilled on arbitrary unstructured grids, ii). the reconstruction is parameter-free, and iii). the finite volume scheme is positivity-preserving when the reconstruction is generalized to the Euler equations. A variety of numerical experiments are presented to demonstrate the performance of this method.

Figures

Tables

Table 1. Table 1: Convergence test for the advection of double sine wave: uniform meshes.

	Cells	$L^1$ error	Order	$L^{\infty }$ error	Order	Reconstruction cost (s)
ILR	16$\times$16$\times$2	4.06e-2	2.33	1.71e-1	1.45	1
	32$\times$32$\times$2	1.16e-2	1.81	6.35e-2	1.43	7
	64$\times$64$\times$2	3.28e-3	1.82	3.00e-2	1.08	64
	128$\times$128$\times$2	8.61e-4	1.93	1.23e-2	1.28	444
MLP	16$\times$16$\times$2	6.19e-2	1.96	2.54e-1	1.13	1
	32$\times$32$\times$2	1.32e-2	2.23	9.70e-2	1.39	6
	64$\times$64$\times$2	2.69e-3	2.29	3.45e-2	1.49	51
	128$\times$128$\times$2	5.71e-4	2.24	1.19e-2	1.54	392
Barth’s limiter	16$\times$16$\times$2	1.36e-1	1.22	4.25e-1	0.82	1
	32$\times$32$\times$2	7.26e-2	0.90	2.65e-1	0.68	3
	64$\times$64$\times$2	3.46e-2	1.07	1.40e-1	0.92	29
	128$\times$128$\times$2	1.82e-2	0.93	9.93e-2	0.50	231
Unlimited	16$\times$16$\times$2	3.40e-2	2.17	7.51e-2	2.01	1
	32$\times$32$\times$2	7.63e-3	2.16	1.69e-2	2.15	3
	64$\times$64$\times$2	1.85e-3	2.05	3.94e-3	2.10	27
	128$\times$128$\times$2	4.58e-4	2.01	9.48e-4	2.05	203

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\nabla \cdot \mathbf{F}(\mathbf{u})=0.$$

$$\frac{\mathrm{d}{\mathbf{u}}_0}{\mathrm{d}t}+\frac{1}{|T_0|}{\oint }_{\partial T_0}\mathbf{F}(\mathbf{u})\cdot \mathbf{n}\mathrm{d}s=0,$$

$${\oint }_{\partial T_0}\mathbf{F}(\mathbf{u})\cdot \mathbf{n}\mathrm{d}s\approx \sum _{j=1}^J|e_j|\mathbf{F}(\mathbf{u}({\mathbf{z}}_j,t))\cdot {\mathbf{n}}_j,$$

$$\frac{\mathrm{d}{\mathbf{u}}_0}{\mathrm{d}t}+\frac{1}{|T_0|}\sum _{j=1}^J\mathcal{F}({\mathbf{u}}_j^{-},{\mathbf{u}}_j^{+};{\mathbf{n}}_j)|e_j|=0.$$

$$\frac{\mathrm{d}{\mathbf{u}}_h}{\mathrm{d}t}=\mathcal{L}({\mathbf{u}}_h)=-\sum _{T_0\in \mathcal{T}}\sum _{j=1}^J\frac{|e_j|}{|T_0|}\mathcal{F}({\mathbf{u}}_j^{-},{\mathbf{u}}_j^{+};{\mathbf{n}}_j)1_{T_0}.$$

$$\{\begin{array}{l}{\mathbf{u}}_h^{\ast }={\mathbf{u}}_h^n+\Delta t_n\mathcal{L}({\mathbf{u}}_h^n),\\ {\mathbf{u}}_h^{n+1}=\frac{1}{2}{\mathbf{u}}_h^n+\frac{1}{2}\left({\mathbf{u}}_h^{\ast }+\Delta t_n\mathcal{L}({\mathbf{u}}_h^{\ast })\right),\end{array}$$

$${\hat{u}}_0(\mathbf{x})=u_0+{\mathbf{L}}^{\top }(\mathbf{x}-{\mathbf{x}}_0),$$

$$\delta (\mathbf{L})=\sum _{j=1}^J({\hat{u}}_0({\mathbf{x}}_j)-u_j{)}^2.$$

$$m_j\le {\hat{u}}_0({\mathbf{z}}_j)\le M_j,j=1,2,\cdots ,J,$$

$$m_j=\min \{u_0,u_j\},M_j=\max \{u_0,u_j\},j=1,2,\cdots ,J.$$

$$\begin{array}{rl}\delta (\mathbf{L})& =\sum _{j=1}^J({\hat{u}}_0({\mathbf{x}}_j)-u_j{)}^2=\sum _{j=1}^J(u_0-u_j+{\mathbf{r}}_j^{\top }\mathbf{L}{)}^2\\ & =\sum _{j=1}^J\left((u_0-u_j{)}^2+2(u_0-u_j){\mathbf{r}}_j^{\top }\mathbf{L}+{\mathbf{L}}^{\top }{\mathbf{r}}_j{\mathbf{r}}_j^{\top }\mathbf{L}\right)\\ & ={\mathbf{L}}^{\top }\mathbf{G}\mathbf{L}+2{\mathbf{c}}^{\top }\mathbf{L}+\text{const,}\end{array}$$

$$m_j\le {\hat{u}}_0({\mathbf{z}}_j)=u_0+({\mathbf{z}}_j-{\mathbf{x}}_0{)}^{\top }\mathbf{L}=u_0+{\mathbf{a}}_j^{\top }\mathbf{L}\le M_j,$$

$$\mathbf{G}=\sum _{j=1}^J{\mathbf{r}}_j{\mathbf{r}}_j^{\top },\mathbf{c}=\sum _{j=1}^J(u_0-u_j){\mathbf{r}}_j,{\mathbf{r}}_j={\mathbf{x}}_j-{\mathbf{x}}_0,{\mathbf{a}}_j={\mathbf{z}}_j-{\mathbf{x}}_0,j=1,2,\cdots ,J.$$

$$\min$$

$$s.t.$$

$$\mathbf{A}={\left[{\mathbf{a}}_1,{\mathbf{a}}_2,\cdots ,{\mathbf{a}}_J\right]}^{\top },\overline{\mathbf{m}}={\left[{\overline{m}}_1,{\overline{m}}_2,\cdots ,{\overline{m}}_J\right]}^{\top },\overline{\mathbf{M}}={\left[{\overline{M}}_1,{\overline{M}}_2,\cdots ,{\overline{M}}_J\right]}^{\top },$$

$${\overline{m}}_j=m_j-u_0,{\overline{M}}_j=M_j-u_0,j=1,2,\cdots ,J.$$

$$\mathrm{Span}\{{\mathbf{r}}_1,{\mathbf{r}}_2,\cdots ,{\mathbf{r}}_J\}={\mathbb{R}}^2,$$

$$u_j^{-}={\hat{u}}_0({\mathbf{z}}_j)=u_0+{\mathbf{L}}^{\top }({\mathbf{z}}_j-{\mathbf{x}}_0),j=1,2,\cdots ,J.$$

$${\delta }_j=\{\begin{array}{ll}+1,& \text{if left-hand side of }j\text{-th constraint is active},\\ -1,& \text{if right-hand side of }j\text{-th constraint is active},\\ 0,& \text{if neither side of }j\text{-th constraint is active}.\end{array}$$

$$\begin{array}{rl}\min & \frac{1}{2}{\mathbf{p}}_k^{\top }\mathbf{G}{\mathbf{p}}_k+(\mathbf{G}{\mathbf{L}}_k+\mathbf{c}{)}^{\top }{\mathbf{p}}_k\\ s.t.& \mathbf{M}{\mathbf{p}}_k=0,\end{array}$$

$${\mathbf{M}}^{\top }\mathbf{\lambda }-\mathbf{G}{\mathbf{p}}_k=\mathbf{G}{\mathbf{L}}_{\mathbf{k}}+\mathbf{c},$$

$$\mathbf{M}{\mathbf{p}}_k=0.$$

$$(\mathbf{M}{\mathbf{G}}^{-1}{\mathbf{M}}^{\top })\mathbf{\lambda }=\mathbf{M}({\mathbf{L}}_k+{\mathbf{G}}^{-1}\mathbf{c}).$$

$$\mathbf{G}{\mathbf{p}}_k={\mathbf{M}}^{\top }\mathbf{\lambda }-\mathbf{G}{\mathbf{L}}_k-\mathbf{c}.$$

$$m_j\le u_0+{\mathbf{a}}_j^{\top }({\mathbf{L}}_k+{\alpha }_k{\mathbf{p}}_k)\le M_j,$$

$$m_j-u_0-{\mathbf{a}}_j^{\top }{\mathbf{L}}_k\le ({\mathbf{a}}_j^{\top }{\mathbf{p}}_k){\alpha }_k\le M_j-u_0-{\mathbf{a}}_j^{\top }{\mathbf{L}}_k.$$

$${\alpha }_k:=\min \left\{1,\underset{1\le j\le J}{\min }{\beta }_j\right\},{\beta }_j=\{\begin{array}{ll}\frac{M_j-u_0-{\mathbf{a}}_j^{\top }{\mathbf{L}}_k}{{\mathbf{a}}_j^{\top }{\mathbf{p}}_k},& {\mathbf{a}}_j^{\top }{\mathbf{p}}_k>0,\\ \frac{m_j-u_0-{\mathbf{a}}_j^{\top }{\mathbf{L}}_k}{{\mathbf{a}}_j^{\top }{\mathbf{p}}_k},& {\mathbf{a}}_j^{\top }{\mathbf{p}}_k<0,\\ +\infty ,& {\mathbf{a}}_j^{\top }{\mathbf{p}}_k=0.\end{array}$$

$${\sigma }_1=\frac{u_1-u_0}{x_1-x_0}=\frac{u_0-u_1}{h},{\sigma }_2=\frac{u_2-u_0}{x_2-x_0}=\frac{u_2-u_0}{h}.$$

$$\min$$

$$s.t.$$

$$2{\sigma }_2^{-}\le L\le 2{\sigma }_2^{+}.$$

$$L^{\ast }=\mathrm{minmod}\left(2{\sigma }_1,2{\sigma }_2,\frac{{\sigma }_1+{\sigma }_2}{2}\right),$$

Full text

\ams

65M08, 65M50, 76M12, 90C20

1 Introduction

The high-order finite volume schemes can be summarized by a reconstruct-evolve-average (REA) process, i.e. a piecewise polynomial is reconstructed in each cell with given cell averages, then the governing equation is evolved according to those polynomials, and finally the cell averages are recalculated. Among these three stages, reconstruction plays an important role in giving high-order solutions without numerical oscillations. Nowadays, second-order methods have been the workhorse for computational fluid dynamics [1]. To achieve second-order accuracy, a prediction of the gradient is obtained first. Due to the possible non-physical flow caused by the underestimation or overestimation of the gradient, it is necessary to limit the gradient in a proper way, and thus the prediction-limiting algorithm arises. In one-dimensional case, total variation diminishing (TVD) limiters are commonly used in designing high-resolution schemes for conservation laws. Unfortunately, it is very difficult to implement the TVD limiter for multi-dimensional problems, especially on unstructured grids. To get around this negative result, a new class of positive schemes has been proposed [2] which ensures a local maximum principle. Since the introduction of this idea, a large number of limiters have then been developed. These limiters include, among others, the Barth’s limiter [3], the Liu’s limiter [4], the maximum limited gradient (MLG) limiter [5] and the projected limited central difference (PLCD) limiter [6]. See [6] for a comprehensive comparison of these limiters. More recently, Park et al. successfully extended the multi-dimensional limiting process (MLP) introduced in [7] from structured grids to unstructured grids [8]. Li et al. proposed the weighted biased averaging procedure (WBAP) limiter [9] based upon the biased averaging procedure (BAP) limiter which was introduced in [10]. Towards the positivity-preserving property, which is crucial for the stability on solving the Euler equations, there have also been several pioneering works. For instance, Perthame and Shu [11] developed a general finite volume framework on preserving the positivity of density and pressure when solving the Euler equations. Motivated by this work, the framework was extended to the discontinuous Galerkin method on rectangular meshes [12] and on triangular meshes [13], and to the Runge-Kutta discontinuous Galerkin method [14]. More recently, a parametrized limiting technique was proposed in [15] to preserve the positivity property on solving the Euler equations on unstructured grids.

Besides the above prediction-limiting algorithm in the reconstruction, there are also methods which deliver the limited gradient in a single process. For example, Chen and Li [16] introduced the concept of integrated linear reconstruction (ILR), in which the limited gradient is computed by solving a linear programming on each cell using an efficient iterative method. A similar approach was given by May and Berger [17]. However, the fulfillment of local maximum principle of these methods requires certain geometric hypothesis on the grids [16]. Buffard and Clain [18] considered a monoslope MUSCL method, where a least-squares problem subjected to maximum principle constraints was imposed. They solved the optimization problem explicitly. However, the involved cases for triangular grids, discussed in the article, were rather complicated, let alone irregular grids with hanging nodes appeared in the computational practice such as mesh adaptation. This motivates us to discard the unsatisfactory geometric hypothesis on unstructured grids by imposing constraints on the quadrature points, and to solve the optimization problem iteratively. In our linear reconstruction, the gradient is indeed obtained by solving a quadratic programming problem using the active-set method. It can be verified that the resulting finite volume scheme for scalar conservation laws satisfies a local maximum principle on arbitrary unstructured grids. Besides, this linear reconstruction can be easily adapted to the Euler equations. And it can be shown that the numerical solutions preserve the positivity of density and pressure.

The remainder of this paper is organized as follows. In Section 2, we review the MUSCL-type finite volume scheme on unstructured grids. In Section 3, we describe the improved integrated linear reconstruction based on solving a series of quadratic programming problems. Section 4 is devoted to the discussion of local maximum principle for scalar conservation laws as well as the positivity-preserving property for the Euler equations. Numerical results are shown for the scalar problems and the Euler equations to validate the effectiveness and robustness of our method in Section 5. Finally, a short conclusion will be drawn in Section 6.

2 MUSCL-type finite volume scheme

In this section we briefly introduce the MUSCL-type finite volume methods for hyperbolic conservation laws on two-dimensional unstructured grids. The extension to three-dimensional cases is straightforward. Consider the following system of hyperbolic conservation laws:

[TABLE]

The computational domain is triangulated into an unstructured grid $\mathcal{T}$. To obtain the MUSCL-type finite volume scheme, we integrate (1) on some control volume $T_{0}\in\mathcal{T}$

[TABLE]

where the cell-averaged solution $\bm{u}_{0}(t)=\dfrac{1}{|T_{0}|}\displaystyle\int_{T_{0}}\bm{u}(\bm{x},t)\mathrm{d}\bm{x}$ and $\bm{n}$ is the unit outward normal over the boundary $\partial T_{0}$ consisting of $J$ edges $e_{1},e_{2},\cdots,e_{J}$.

Since we are concerned with second-order schemes, the boundary integral appeared in (2) is approximated by the midpoint quadrature rule, namely,

[TABLE]

where $\bm{z}_{j}$ and $\bm{n}_{j}$ represent the midpoint and normal of the edge $e_{j}~{}(j=1,2,\cdots,J)$. To maintain stability, the flux term $\bm{F}(\bm{u}(\bm{z}_{j},t))\cdot\bm{n}_{j}$ appeared in (3) is replaced by a numerical flux function $\mathcal{F}(\bm{u}_{j}^{-},\bm{u}_{j}^{+};\bm{n}_{j})$, where $\bm{u}_{j}^{\pm}$ denote values of numerical solutions at $\bm{z}_{j}$ outside/inside the cell $T_{0}$ respectively $(j=1,2,\cdots,J)$. The resulting semi-discrete scheme is

[TABLE]

To achieve second-order numerical accuracy, the values $\bm{u}_{j}^{\pm}$ are evaluated from the piecewise linear solution reconstructed from nearby cell averages. In Section 3 we will describe how to construct an appropriate piecewise linear solution which satisfies certain properties such as local maximum principle or positivity-preserving property. After that, we obtain a system of ordinary differential equations for $\bm{u}_{h}=\displaystyle\sum_{T_{0}\in\mathcal{T}}\bm{u}_{0}\mathbb{1}_{T_{0}}$:

[TABLE]

Finally, the system (4) is discretized by the SSP Runge-Kutta method to achieve second-order temporal accuracy [19]:

[TABLE]

where the time step length $\Delta t_{n}$ is determined through the CFL condition.

To effectively control the total amount of cells, we adopt the $h$-adaptive method following Li [20]. This method relaxes the regularity requirement of the mesh. Therefore, even if the background mesh is a regular mesh where any two adjacent cells share a whole edge, an irregular mesh may be produced during the refinement process. The irregular structure may lead to different resolution of the adjacent cells used for computing numerical fluxes. Fortunately, the strategy of mesh refinement proposed in [20] prevents the adjacent cells to differ more than one level of resolution, and thus there are essentially three possibilities which are depicted schematically in Fig. 1. The shaded triangle appeared in case b) is called twin triangle, since it can be divided into two smaller triangles by one of its medians. The twin triangle has one special edge shared by two neighbors. However, if we treat the twin triangle as a quadrilateral with $J=4$ edges, then all the computations can be carried out in the usual way.

3 Construction of the integrated linear reconstruction

In this section we describe the improved integrated linear reconstruction. Consider the solution $u(\bm{x},t)$ of a scalar conservation law. For simplicity we will omit the dependency in time throughout this discussion. The reconstructed linear function on a given cell $T_{0}\in\mathcal{T}$ with the centroid $\bm{x}_{0}$ can be formulated as

[TABLE]

where $\bm{L}$ is the gradient vector.

To compute the gradient $\bm{L}$ we require information from neighbors of $T_{0}$. In the interior of the computational domain the von Neumann neighbors are sufficient for reconstruction. Nevertheless, this is not the case on the boundary, where we need special treatment, as will be described in Section 3.3.

To achieve better accuracy without non-physical oscillation, we minimize the residuals over all gradient $\bm{L}$ subjected to some stability conditions. The solution of this optimization problem, if exists, is used to construct a piecewise linear solution. This is the rough idea of the integrated linear reconstruction.

3.1 Formulation of optimization problem

In the integrated linear reconstruction proposed in this paper, the objective function is the sum of squared residuals, rather than the sum of absolute residuals that was used in [16],

[TABLE]

Potentially, the smoothness of the reconstructed solutions can be improved. In [16], it has been shown that to theoretically guarantee the maximum principle, a certain hypothesis on the geometry of elements needs to be satisfied. The hypothesis can be smoothly satisfied with quality Delaunay triangular mesh, which can be generated by most mature tools such as EasyMesh [21]. However, it can be shown that the hypothesis can be easily violated on distorted meshes. To resolve this issue, we discard the constraints on the neighboring centroids and impose the following inequalities on the midpoint $\bm{z}_{j}$:

[TABLE]

where the lower and upper bounds are given by

[TABLE]

Now using the definition (6) of linear function $\hat{u}_{0}(\bm{x})$, the objective function (7) becomes

[TABLE]

and the constraints (8) become

[TABLE]

where

[TABLE]

Consequently, the linear reconstruction reduces to the following double-inequality constrained quadratic programming (QP) problem

[TABLE]

where

[TABLE]

It is reasonable to assume that

[TABLE]

which is true for most of the triangular grids in practical computations. Then the matrix $\mathbf{G}$ is positive definite and thus the existence and uniqueness of the minimizer $\bm{L}$ are guaranteed. With the gradient $\bm{L}$ in hand, we build the interior value $u_{j}^{-}$ as follows:

[TABLE]

In the following discussion we develop an efficient solver for the problem (9).

3.2 An active-set method for quadratic programming problem

The active-set method has been widely used since the 1970s and is effective for small- and medium-sized QP problems. This method updates the solution by solving a series of QP subproblems in which some of the inequality constraints are imposed as equalities. As for the double-inequality constrained problem (9), any constraint is assumed to be active on a single side on each step in the active-set method, even if both sides of this constraint hold with equality. Therefore, we introduce a balanced ternary variable $\delta_{j}$ to indicate the active state of $j$-th constraint, i.e.

[TABLE]

We generate iterates of (9) that remain feasible while steadily decreasing the objective function. To be more precise, given $\bm{L}_{k}$ obtained from the $k$-th iteration, we solve the following QP problem for the descending direction $\bm{p}_{k}$,

[TABLE]

where $\mathbf{M}=[\delta_{j}\bm{a}_{j}^{\top}]_{\delta_{j}\neq 0}$ is composed of normals of the active constraints.

The first-order necessary conditions for $\bm{p}_{k}$ to be a solution of (10) imply that there is a vector of Lagrange multipliers $\bm{\lambda}$ such that

[TABLE]

Multiply (11a) by $\mathbf{M}\mathbf{G}^{-1}$ and then add (11b) to obtain a linear system in the vector $\bm{\lambda}$ alone:

[TABLE]

We solve this symmetric positive definite system for $\bm{\lambda}$. Then $\bm{p}_{k}$ can be recovered from

[TABLE]

If the descending direction $\bm{p}_{k}$ is nonzero, we set $\bm{L}_{k+1}=\bm{L}_{k}+\alpha_{k}\bm{p}_{k}$ where the step-length parameter $\alpha_{k}$ is chosen to be the largest value in the range $[0,1]$ such that all constraints are satisfied. Indeed an explicit formula for $\alpha_{k}$ can be derived. To satisfy the $j$-th constraint, we have

[TABLE]

or

[TABLE]

Therefore, to maximize the decrement of objective function, the step-length parameter $\alpha_{k}$ can be chosen in the following way:

[TABLE]

If $\alpha_{k}<1$, then there exists an index $j$ such that $\beta_{j}=\alpha_{k}$. In this case the active indicator $\delta_{j}$ is adjusted to $-\operatorname{sgn}(\bm{a}_{j}^{\top}\bm{p}_{k})$ accordingly.

On the other hand, if $\bm{p}_{k}$ is zero, then we check the signs of the Lagrange multipliers. If all the multipliers are non-negative, then we have achieved optimality. Otherwise, we can find a feasible direction by dropping the constraint corresponding to the most negative multiplier.

The initial iterate $\bm{L}_{0}$ can be any feasible solution of (9) such as the null solution $\bm{L}_{0}=\bm{0}$. And the initial active indicators $\delta_{j}$ are all set to zero. The whole algorithm is sketched out in Appendix A. We refer the readers to [22] for more details about the active-set methods.

Remark 3.1** (Connection with MC limiter).**

Let us investigate a special case: the ILR on a one-dimensional grids with uniform spacing $h$. For each cell $T_{0}$, let the left and right neighbors be labeled with 1 and 2 respectively. Furthermore, the backward and forward slopes are denoted by $\sigma_{1}$ and $\sigma_{2}$ respectively, i.e.

[TABLE]

The problem (9) then reduces to the following univariate QP problem

[TABLE]

The optimal solution of (12)

[TABLE]

is the well-known MC (monotonized central-difference) limiter, where the minmod function

[TABLE]

Remark 3.2** (Connection with Barth’s limiter).**

If the gradient is further restricted to the direction of unlimited gradient, i.e. $\bm{L}=-\phi\mathbf{G}^{-1}\bm{c}$ for some $\phi\in[0,1]$, then we arrive at the Barth’s limiter [3]. A simple derivation yields the following explicit expression for $\phi$

[TABLE]

A comparison of ILR against Barth’s limiter will be provided in Section 5.

3.3 Boundary treatment

For the cells located around the domain boundary, the proposed reconstruction can not be implemented smoothly according to the procedure introduced in the previous subsection. On the one hand, for those cells with at least one edge on the boundary, the von Neumann neighbors are inadequate to give a reasonable reconstruction (Fig. 2 (a)). Thus, besides the von Neumann neighbors, we include those cells sharing at least one vertex with the current cell, which introduces the so-called Moore neighbors $T_{M_{1}},T_{M_{2}},\cdots,T_{M_{K}}$ (Fig. 2 (b)). On the other hand, it is hard to specify a suitable lower and upper bounds for edges on the boundary when the solution from the outside of the computational domain is unavailable. Here we adjust the bounds using solutions from the entire stencil, i.e.

[TABLE]

The resulting QP problem is still (9), with a new definition of $\mathbf{G}$ and $\bm{c}$

[TABLE]

Now the flux computation can be computed in the usual way. The boundary condition is imposed on the midpoint of the boundary edge. When periodic boundary condition is prescribed, no special treatment is needed.

4 Local maximum principle and positivity-preserving property

4.1 Local maximum principle

A more appealing feature of our reconstruction is that the resulting finite volume scheme satisfies a local maximum principle. We restrict ourselves to first-order Euler forward in time

[TABLE]

Second order SSP Runge-Kutta time discretization (5) will definitely keep the validity of the local maximum principle since it is a convex combination of two Euler forward steps. For clarity, the labels of edges are ordered in the way in Fig. 3. Here we have the following theorem on local maximum principle. The technique to prove the local maximum principle is similar to that of [13].

Theorem 4.1** (Local maximum principle).**

Suppose that $\mathcal{F}$ is a monotone Lipschitz continuous numerical flux function. Let $\mathcal{T}$ be a triangulation mixed with triangles and twin triangles. If for each cell $T_{0}\in\mathcal{T}$ the linear reconstruction satisfies

[TABLE]

then the finite volume scheme (13) fulfills the local maximum principle

[TABLE]

under the CFL-like condition

[TABLE]

where $h$ is the minimum size of cells measured by their diameters of inscribed circles.

Proof 4.2**.**

From Fig. 3 we have the following formula

[TABLE]

As a result, the finite volume scheme (13) can be split as

[TABLE]

with

[TABLE]

if $T_{0}$ is a triangle, or

[TABLE]

with

[TABLE]

if $T_{0}$ is a twin triangle.

It follows that the derivatives of $H_{1},\cdots,H_{J}$ with respect to their arguments $u_{1}^{\pm},\cdots,u_{J}^{\pm}$ are all non-negative provided that

[TABLE]

or, using the definition of $h$,

[TABLE]

Therefore, if we view the right-hand side of the scheme (15) or (16) as a function of $u_{1}^{\pm},\cdots,u_{J}^{\pm}$

[TABLE]

then the function $H$ is non-decreasing with respect to each argument. Moreover, we have $H(u,\cdots,u)=u$ due to the consistency of numerical flux function.

Denote

[TABLE]

then the trace values satisfy the inequalities

[TABLE]

Consequently the monotonicity of $H$ implies the local maximum principle

[TABLE]

Remark 4.3**.**

In Theorem 4.1, we focus on the case with triangular meshes generated by the adaptive mesh refinement. Indeed the local maximum principle is valid on any polygon meshes, as long as a quadrature formula like (14) can be found.

4.2 Positivity-preserving property

In this part we extend the previous framework of scalar conservation laws to the system of the Euler equations of ideal gases:

[TABLE]

[TABLE]

where $\rho$ is the density, $\bm{v}=(u,v)$ is the velocity, $E$ is the total energy, $p$ is the pressure and $\gamma>1$ is a constant.

The local maximum principle is no longer applicable for system of conservation laws. However, the positivity-preserving property is desired in the computation of the Euler equations. It means that the solution $\bm{u}$ belongs to the set of admissible states

[TABLE]

To simplify the discussion we consider the finite volume scheme

[TABLE]

with local Lax-Friedrichs flux

[TABLE]

where $c=\sqrt{{\gamma p}/{\rho}}$ is the speed of sound.

To guarantee positivity, we perform the reconstruction with respect to the primitive variables $\rho,u,v$ and $p$ respectively. The following theorem states that the finite volume scheme will preserve the positivity of numerical solutions.

Theorem 4.4** (Positivity-preserving).**

Let $\mathcal{T}$ be a triangulation mixed with triangles and twin triangles. Labels of edges are illustrated in Fig. 3. Suppose that for each cell $T_{0}\in\mathcal{T}$ the linear reconstruction with respect to the primitive variables satisfies $\bm{u}_{j}^{-}\in\mathcal{G}$ for $j=1,2,\cdots,J$. Then the finite volume scheme (17) with local Lax-Friedrichs flux (18) satisfies $\bm{u}_{0}^{n+1}\in\mathcal{G}$ under the CFL condition

[TABLE]

where $h$ has the same meaning as Theorem (4.1) and $0<\beta\leq 1$ is sufficiently small such that

[TABLE]

for any triangle $T_{0}$ and

[TABLE]

for any twin triangle $T_{0}$.

Proof 4.5**.**

We only prove the case with triangle here. The case with twin triangle can be carried out in a similar manner. Substituting (18) into (17) and rearranging terms yield

[TABLE]

where $\lambda=a\Delta t_{n}/2|T_{0}|$.

It can be verified that the set $\mathcal{G}$ is a convex cone. Note that $\lambda(|e_{1}|+|e_{2}|+|e_{3}|)={2a\Delta t_{n}}/{h}\leq\beta$. Since $\bm{u}_{1}^{-},\bm{u}_{2}^{-},\bm{u}_{3}^{-}\in\mathcal{G}$ and $\bm{u}_{0}^{n}-\lambda(|e_{1}|+|e_{2}|+|e_{3}|)(\bm{u}_{1}^{-}+\bm{u}_{2}^{-}+\bm{u}_{3}^{-})\in\mathcal{G}$, to guarantee $\bm{u}_{0}^{n+1}\in\mathcal{G}$ it suffices to verify that

[TABLE]

for any $\bm{u}\in\mathcal{G}$ and normal vector $\bm{n}$, that is, to check the positivity of the density and pressure. Actually,

[TABLE]

where $\rho(\bm{u})=\rho$ and $p(\bm{u})=(\gamma-1)\left(E-\dfrac{m^{2}+n^{2}}{2\rho}\right)$ for $\bm{u}=[\rho,~{}m,~{}n,~{}E]^{\top}$.

Remark 4.6**.**

Since $\bm{u}_{1}^{-},\bm{u}_{2}^{-},\cdots,\bm{u}_{J}^{-}$ are close to the cell average $\bm{u}_{0}^{n}$, we can expect that the value of the constant $\beta$ in (19) is close to $1/3$.

Remark 4.7**.**

Assume that at $n$-th time level $\bm{u}_{h}^{n}\in\mathcal{G}$. Then for any cell $T_{0}\in\mathcal{T}$, the primitive variables $\rho_{0}^{n}$ and $p_{0}^{n}$ must be positive. The constraints of ILR guarantee the positivity of trace values $\rho_{j}^{-}$ and $p_{j}^{-}$ for all $j$. As a result, $\bm{u}_{j}^{-}\in\mathcal{G}$ for all $j$. From Theorem 4.4 we conclude that $\bm{u}_{0}^{n+1}\in\mathcal{G}$ and hence at $(n+1)$-th time level $\bm{u}_{h}^{n+1}\in\mathcal{G}$. Therefore, the ILR for the Euler equations leads to positive numerical solutions as long as the initial solution is positive.

5 Numerical results

In this section we present several numerical examples to demonstrate the numerical performances of our integrated linear reconstruction on various unstructured grids.

5.1 Linear advection problem

We consider the linear advection equation

[TABLE]

with initial profile given by the double sine wave function [6, 8, 17, 16]

[TABLE]

The computational domain is $[0,1]\times[0,1]$. The periodic boundary condition is applied. We perform the convergence test on both uniform and non-uniform meshes (see Fig. 4 for the coarse mesh). The uniform mesh is generated by dividing rectangular cells along the diagonal direction, while the non-uniform mesh is generated by the Delaunay triangulation of EasyMesh [21]. The upwind flux is adopted as the numerical flux and the CFL number is $0.3$. We present quantitative comparisons of solution errors as well as the reconstruction cost until $t=1$ in Tables 2 and 2, where we include the MLP limiter (version u1) [8] as a representative example of limiters involving the Moore neighbors. It is found that even though MLP is slightly better than our ILR, they behave quite similar in the accuracy and efficiency. As a comparison, Barth’s limiter has the advantage on the efficiency, but it only gives rise to first-order accuracy. Moreover, it was observed in this numerical test that the iteration of QP solver converges in roughly two steps on average, which confirms the excellent performance of ILR.

5.2 Solid body rotation on the stretched mesh

To assess the behavior of ILR for non-uniform scalar flow field against an anisotropic mesh, we consider the following solid body rotation problem proposed by LeVeque [23]

[TABLE]

on a highly-stretched mesh in the computational domain $[0,1]\times[0,1]$ (Fig. 8). The highly-stretched mesh can be generated in the following three steps [24]:

A rectangular mesh is stretched toward the horizontal line $y=0.5$ by a factor $\beta=1.2$; 2. 2.

Irregularities are introduced by random shifts of interior nodes in $x$-directions; 3. 3.

Each distorted quadrilateral is further divided along its diagonal direction.

The initial profile consists of smooth hump, cone and slotted cylinder. These shapes are located within the circle of radius $r_{0}=0.15$ centered in $(x_{0},y_{0})$. In the other region the initial value vanishes. Denote by $r(x,y)=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}/r_{0}$ the normalized distance. The slotted cylinder is centered at $(x_{0},y_{0})=(0.5,0.75)$ and

[TABLE]

the sharp cone is centered at $(x_{0},y_{0})=(0.5,0.25)$ and

[TABLE]

and the smooth hump is centered at $(x_{0},y_{0})=(0.25,0.5)$ and

[TABLE]

In the computation we use a fixed time step $\Delta t=0.6\pi h$. The homogeneous Dirichlet boundary conditions are prescribed. To measure the oscillation of the numerical solution we compute the discrete total variation (TV) of the piecewise linear solution $\hat{u}_{h}=\displaystyle\sum_{T_{0}\in\mathcal{T}}\hat{u}_{0}\mathbb{1}_{T_{0}}$ in the following way

[TABLE]

The evolutions of the discrete total variation are depicted in Fig. 8 for different resolutions and in Fig. 8 for different methods. From here we can see that the total variation is decreasing as a whole except for individual tiny fluctuations when the slot rotates towards the stretching line $y=0.5$. The total variation will flatten out as mesh refines. Moreover, our ILR has the best performance on preserving the total variation of the solution on the highly-stretched mesh among the three methods. The snapshots in Fig. 8 show the shapes of the solutions with different methods at $t=2\pi$, corresponding to a complete rotation. We can observe that the result of ILR captures the initial profile much better than that of MLP or Barth’s limiter.

5.3 Double Mach reflection

As one benchmark of the Euler equations, the double Mach reflection problem [25] is considered here. The whole computational domain is $[0,4]\times[0,1]$. Initially, a right-moving Mach 10 shock is located at the beginning of the wall $(x,y)=(1/6,0)$, making a $60^{\circ}$ angle with the $x$-axis. The boundary setup is the same to [8]. The HLLC flux is used as the numerical flux and computation is carried out until $t=0.2$.

Fig. 10 shows the comparison of density contours computed on two successively refined Delaunay meshes. Obviously ILR provides a more stable structure for the Mach stem and slip line than Barth’s limiter does. Fig. 10 shows the close-up view around the Mach stem. Again ILR provides better resolution below the Mach stem, though the Barth’s limiter captures the shear layer instability from the shock triple point more accurately.

5.4 A Mach 3 wind tunnel with a step

This is another popular test for high-resolution schemes [25]. The problem begins with a uniform Mach 3 flow in a wind tunnel with a step of 1 length unit high and 3 length units long. The step is $0.2$ length units high and is located $0.6$ length units from the left end of the tunnel. Reflective boundary condition is applied along the walls of the tunnel. And inflow and outflow boundary conditions are used at the left and right boundaries respectively. To resolve the singularity at the corner of the step, we refine the initial mesh near the corner, as is shown in Fig. 13.

To capture the shock structure efficiently we apply the mesh adaptation. The adaptive indicator we use here is the jump of pressure. And the adaptation is carried out on each time level. The HLLC flux is used for this problem. Figs. 13 and 13 show the adaptive mesh and contour picture for the density at $t=4$. It can be observed that the ILR captures the discontinuities well and provides a high resolution of the slip line from the shock triple point.

5.5 Sedov point-blast wave

The Sedov point-blast wave is a self-similar problem with both low density and low pressure. The analytical solution of this problem is available and can be found on [26, 27]. Without positivity-preserving treatment the computation may break down due to the presence of negative density or pressure. Here, we regard this spherically symmetric problem as three-dimensional cylindrically symmetric flow. The governing equation is

[TABLE]

where

[TABLE]

with

[TABLE]

Let $\bm{U}=r\bm{u}$, then formally we obtain the Euler equations with a source term

[TABLE]

The finite volume discretization of (21) we apply is

[TABLE]

where $(w_{1},w_{2},w_{3})=(1/3,1/3,1/3)$ for $J=3$, and $(w_{1},w_{2},w_{3},w_{4})=(1/3,1/3,1/6,1/6)$ for $J=4$. It can be proved that the scheme (22) also preserves positivity provided that the time step length satisfies $a\Delta t_{n}\leq\beta h/4$. See Appendix B for details of the proof.

In the blast wave setup, a finite amount of energy $E_{0}=0.851072$ is released at the origin at $t=0$. The initial density $\rho_{0}=1$. These values are chosen so that the shock arrives at $r=1$ at the chosen time $t=1$ [27]. A triangular grid on $[0,1.2]\times[0,1.2]$ with uniform spacing $l=1.2/80$ is used for computation. Symmetric boundary conditions are imposed at the bottom and left edges. The total energy $E$ is a constant $E_{0}/2\pi l^{3}$ in the two cells at the left-bottom corner and $10^{-12}$ elsewhere (emulating a $\delta$-function at the origin). Again we use the HLLC flux as numerical flux. The contour lines of density and the profile along the radial direction are displayed in Fig. 14, which demonstrate the robustness and stability of ILR.

6 Conclusion

We present an integrated linear reconstruction for finite volume scheme on arbitrary unstructured grids. We use the sum of squared residuals as an objective function, and the inequalities on midpoints as linear constraints. The resulting optimization problem is a convex quadratic programming problem which can be efficiently solved using the active-set methods. Numerical tests are provided to show the capacity of this new technique to deal with computations on a variety of unstructured meshes. Extension to higher-order reconstruction is possible, by increasing the number of variables and constraints in the quadratic programming problems.

Acknowledgements

The first and third authors gratefully acknowledge the support of the National Natural Science Foundation of China (Grant No. 91330205, 11421110001, 11421101 and 11325102). The second author would like to thank the support of FDCT 050/2014/A1 from Macao SAR, MYRG2014-00109-FST from University of Macau, and the National Natural Science Foundation of China (Grant No. 11401608).

Appendix A Procedures of active-set methods

The procedures of active-set method for the double-inequality constrained quadratic programming (9) are listed in the following Algorithm 1.

There are two aspects regarding the numerical issue in our actual implementation. Firstly, although the matrix $\mathbf{M}$ in the active-set method is theoretically full rank at each iteration, we cannot guarantee it numerically. Therefore, the condition $\bm{a}_{j}^{\top}\bm{p}\gtrless 0$ in the computation of the coefficients $\beta_{j}$ is relaxed to $\bm{a}_{j}^{\top}\bm{p}\gtrless\epsilon U$ accordingly, where $\epsilon=10^{-12}$ and $U$ denotes the reference magnitude of the solution. Similarly, the termination condition of the Lagrange multipliers is relaxed to $\bm{\lambda}\geq\epsilon U$. Secondly, the maximum number of iterations is chosen as $6$, which is adequate for most two-dimensional numerical tests.

Appendix B Positivity of axisymmetric Euler equations

The positivity of axisymmetric Euler equations can be verified using the technique in [28]. For simplicity, we only consider the purely triangular grid. The finite volume scheme (22) can be rewritten as

[TABLE]

where $\tilde{\bm{U}}\in\mathcal{G}$ is the sum of last nine terms in (20), with $\bm{U}_{1}^{-},\bm{U}_{2}^{-}$ and $\bm{U}_{3}^{-}$ in place of $\bm{u}_{1}^{-},\bm{u}_{2}^{-}$ and $\bm{u}_{3}^{-}$ respectively. Therefore, to guarantee the positivity of $\bm{U}_{0}^{n+1}$, the time step length $\Delta t_{n}$ satisfies

[TABLE]

In the following derivation we drop the superscripts and subscripts. Let $\tau=2\Delta t_{n}/3\beta$. Obviously $\rho(\bm{U}+\tau\bm{s}(\bm{U},\bm{z}))=r\rho>0$. A direct algebraic calculation shows that

[TABLE]

and as a result,

[TABLE]

Since

[TABLE]

we obtain the following CFL-like condition

[TABLE]

where $r_{\min}$ is the minimum radius of quadrature points under consideration.