- Email: [email protected]
Limiting strategies for polynomial reconstructions in the finite volume approximation of the linear advection equation
Applied Numerical Mathematics 49 (2004) 277–289 www.elsevier.com/locate/apnum
Limiting strategies for polynomial reconstructions in the finite volume approximation of the linear advection equation Enrico Bertolazzi a,∗ , Gianmarco Manzini b a Dipartimento di Ingegneria Meccanica e Strutturale, Università di Trento, via Mesiano 77, I-38050 Trento, Italy b Istituto di Analisi Numerica IAN-CNR, via Ferrata 1, I-27100 Pavia, Italy
Abstract A cell-centered semi-discrete finite volume method is proposed to accurately solve the time-dependent scalar advection equation. The spatial accuracy is ensured by a piecewise linear reconstruction which requires a suitable limiting strategy to control spurious numerical oscillations. Three different approaches are conceived to limit the approximate solution. 2004 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Finite volume; Non-oscillatory reconstructions; Limiters
1. Introduction The present work details with the numerical approximation of the solution of the linear advection of a scalar quantity C in a constant velocity field V. The construction of an accurate non-oscillatory scheme is presented for the two-dimensional equation; schemes for one and three space dimensions can be similarly derived. First-order accurate schemes are developed in the finite volume framework by using a piecewise constant representation of the approximate solution within mesh cells. Higher-order approximations are usually achieved by implementing polynomial reconstructions. Nonetheless, it is a well-known fact that even the simplest linear schemes are oscillatory near discontinuities. An effective way of controlling the process of formation of spurious numerical oscillations consists in introducing a solution dependent limiter in the reconstruction algorithm. That is, to avoid the limitation of the Godunov theorem [1], appropriate non-linear schemes must also be devised for linear problems. There exists a large literature about non-linear schemes for hyperbolic conservation * Corresponding author.
E-mail address: [email protected] (E. Bertolazzi). 0168-9274/$30.00 2004 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2003.12.007
278
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
laws; we mention, for example, the total variation diminishing (TVD) methods [2–7] and the essentially non-oscillatory (ENO/WENO) methods [8–11]. The aim of this paper is to discuss some limiting strategies for the linear reconstruction in the framework of second-order accurate finite volume approximation methods. To this purpose, we propose and analyze three different limiting constraints for which some important semi-discrete properties are maintained. An important issue that we want also to focus on, is that with proper initial and boundary conditions an initially non-negative scalar solution C must remain non-negative. This property is crucial when C represents a non-negative quantity (i.e., mass density) because negative solutions are physically unacceptable. The model problem is described in Section 2, while the finite volume formulation and the piecewiselinear reconstruction algorithm are detailed in Sections 3 and 4. The first limiting condition is discussed in Section 5; then the method is reformulated in a more compact matrix-like form in Section 6. Two stronger limiting constraints are introduced and their consequences analyzed in Section 7. The effectiveness of our approach is shown by a set of numerical experiments in Section 8, followed by the conclusions in Section 9.
2. The scalar advection equation Our model problem is the two-dimensional scalar advection equation, which reads as Ct + ∇ · (CV) = 0
on Ω,
(1)
where C(t, x) is a scalar quantity advected by an assigned steady velocity field V throughout a polygonal domain Ω. The boundary of Ω, indicated by ∂Ω, can be split into an inflow and an outflow part defined by inflow boundary:
Γ − = {x ∈ ∂Ω | V · n < 0},
outflow boundary:
Γ + = {x ∈ ∂Ω | V · n 0},
and such that ∂Ω = Γ − ∪ Γ + . The model problem described by Eq. (1) is completed by a suitable initial condition C(0, x) = C0 (x),
x ∈ Ω,
and a Dirichlet boundary condition on the inflow boundary, that is C(t, x) = g(x),
(t, x) ∈ R+ × Γ − .
Proposition 1. Under the assumption g(x) 0,
x ∈ ∂Γ − ,
C0 (x) 0,
x ∈ Ω,
the following properties hold:
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
(a) (b) (c)
0 C(t, x) M, d CL1 (Ω) + γ C, vn = 0, dt d C2L2 (Ω) + (γ C)2 , vn = 0, dt
279
(2)
where
M = max C0 L∞ (Ω) , gL∞ (Γ − ) .
Remark 2. The symbol u, v denotes the familiar scalar product of two scalar functions u and v, i.e., u, v = uv ds, ∂Ω
γ is the trace operator and vn = V · n.
3. Discrete spaces and norms Let Ω be a polygonal domain and {Th } a family of triangulations assumed, as usual in finite element approximation, regular and conformal in the sense specified by Ciarlet for h → 0. The discretization parameter is the maximum diameter of the cells in Th . Any triangulation Th is also supposed weakly acute, i.e., the maximum angle allowed for a triangle is π/2. We shall use the notations and conventions listed above. • Triangles are conventionally labelled by an integer index. The generic triangle Ti has area |Ti | and centroid xi . • The positive 1-dimensional Lebesgue measure intersection of two triangles or of a triangle and the border ∂Ω is an edge. The former case defines an internal edge, the latter case a boundary edge. The internal edge shared by triangles Ti and Tj is denoted by eij . For the sake of notation consistency, a boundary edge is also denoted by the same symbol eij . In this latter case i is the unique triangle the edge belongs to and j a specific boundary identifier like a fictitious external triangle. • eij has 1-dimensional Lebesgue measure (length) |eij |, mid-point xij and unit orthogonal vector nij . This latter is arbitrarily oriented from Ti to Tj if eij is an internal edge, and outward directed if eij is on the boundary. • Eh and Th are, respectively, the set of the mesh edges and triangles. Eh denotes the subset of mesh edges belonging to the boundary ∂Ω. • Th (i) is the index set of the triangles adjacent to Ti and Th (i) the index set of the “fictitious triangles”. Let us now introduce the finite dimensional vector spaces of piecewise polynomial functions Pk Eh = f : f |eij ∈ Pk (R) , Pk (Th ) = f : f |Ti ∈ Pk R2 , whose elements are indicated by bold letters. In particular, fi indicates the value f (Ti ) when f ∈ P0 (Th ) and fij the value f(eij ) when f ∈ P0 (Eh ).
280
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
In P0 (Th ) and P0 (Eh ) we introduce, respectively, the mesh-dependent norm and the mesh-dependent scalar product
1/p p |Ti ||fi | , f, g = |eij |fij gij , fp = eij ∈Eh
T i ∈ Th
which are discrete versions of the Lp norm and the L2 scalar product. Let us also define the two projection maps 1 p 0 π(f )|Ti = πi (f ) = f (x) dx, π : L (Ω) → P (Th ), |T i | Ti 1 p − 0
f (x) dx. π(f )|eij = πij (f ) = π : L (Γ ) → P Eh , |eij | eij
The linear reconstruction procedure on unstructured meshes can be viewed as the mapping R : P0(Th ) → Pk (Th ),
R(f)|Ti = Ri (f),
Ti ∈ T h ,
(3)
where Ri (f)(x) = fi + Gi (f) · (x − xi ),
x ∈ Ti .
(4)
The term Gi (f) denotes the cell-centered reconstructed gradient. It is worth noting that this definition is consistent with the initial cell-averaged approximate values in the sense that cell-averages are preserved by the reconstruction process. Formally, for any f ∈ Lp (Ω) we have 1 R π(f ) (x) dx, Ti ∈ Th , (5) πi (f ) = |T i | Ti
which is equivalent to set π ◦ R = id.
4. The finite-volume formulation The finite-volume approximation of Eq. (1) logically proceeds in two steps. In the first step, Eq. (1) is re-formulated on the triangulated domain by integrating within each control volume Ti and then by applying the Gauss divergence theorem [12]. We have d C dx + CV · nij ds = 0, Ti ∈ Th . (6) dt
Ti
j ∈Th (i)∪Th (i)eij
The second step consists in approximating (6) by the following discrete relation dci + Ghij = 0, Ti ∈ Th , |T i | dt
(7)
j ∈Th (i)∪Th (i)
which holds for every t > 0. In Eq. (7), ci ≈ πi (C) approximates the cell-averaged value of the advected quantity C, while Ghij is the discrete flux integral which approximates both the internal and the boundary
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
281
Fig. 1. Definition of qij .
edge advective flux h Gij ≈ C(·, x)V · nij ds,
j ∈ Th (i) ∪ Th (i).
(8)
eij
Line integrals (8) are computed by applying a single-point quadrature rule to the standard upstream numerical flux [7]. The numerical flux integral takes the form Ghij (c) = |eij |Ri (c)(qij )vij+ + |eij |Rj (c)(qij )vij− ,
j ∈ Th (i)
(9)
at internal edges and the form Ghij (c) = |eij |Ri (c)(qij )vij+ + |eij |g(qij )vij− ,
j ∈ Th (i)
(10)
at boundary edges. In Eqs. (9) and (10) we have set vij = V · nij and vij± = (vij ± |vij |)/2. In the case of an internal edge eij , the quadrature point qij is defined as the intersection of eij and the segment connecting xi to xj , see Fig. 1. When a boundary edge is considered, qij is the edge middle point. The discrete initial solution is given by projecting the analytical initial solution onto the piecewise constant space P0 (Th ), i.e., ci (0) = πi (C0 ),
Ti ∈ T h .
Two discrete trace operators can be defined in P0 (Th ) and P1 (Th ). These are γ 0 : P0 (Th ) → P0 Eh , γij0 (c) = cij , and
γ : P (Th ) → P Eh , 1
1
0
γij1 (c)
=
g(qij ) Ri (c)(qij )
if eij ⊂ Γ − , otherwise.
5. A first condition on the reconstruction In order to control the numerical oscillations, we assume that a limiter is introduced into the reconstruction procedure [13]. The limited reconstructed solution within Ti is again denoted by Ri . A first condition of the reconstructed solution is that the sign of the slope in the direction qij − xi is the same of the jump cj − ci (or γij1 (c) − ci for the boundary); thus we assume that sign Ri (c)(qij ) − ci = sign(cj − ci ),
282
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
for internal edges, while sign Ri (c)(qij ) − ci = sign γij1 (c) − ci , for boundary edges. This condition implies that there exists a non-negative scalar coefficient σij (c) such that: cj − ci , if j ∈ Th (i), (COND.1) Ri (c)(qij ) = ci + σij (c) × 1 γij (c) − ci if j ∈ Th (i). 6. The semi-discrete formulation Assuming condition (COND.1) and introducing the scalar quantity v˜ ij (c) = vij+ σij (c) + vj+i σj i (c), in (9) the internal edge numerical flux integral becomes Ghij (c) = |eij |vij+ ci − |eij |vj+i cj + |eij |v˜ij (c)(cj − ci ),
j ∈ Th (i).
The internal flux balance of the cell Ti can be written in a matrix–vector product form Ghij (c) = Gc − G(c)c |i , Ti ∈ T h , j ∈Th (i)
by introducing + if j ∈ Th (i), −|eij |vj i + G|ij = |eik |vik if i = j, k∈Th (i) 0 otherwise, if j ∈ Th (i), −|eij |v˜j i (c) G|ij (c) = |e |v˜ (c) if i = j, k∈Th (i) ik ik 0 otherwise. Similarly, condition (COND.1) in (10) yields the boundary edge numerical flux integral in the form Ghij (c) = |eij |vij γij1 (c) + |eij |vj+i 1 − σij (c) ci − γij1 (c) , j ∈ Th (i), and its the contribution to the numerical flux balance of Ti takes the more compact vector form Ghij (c) = bi (c), Ti ∈ Th . j ∈Th (i)
Consequently, the FV method can be re-formulated as dc = T−1 G(c)c − Gc − b(c) , dt where Tij = δij |Ti | is the mass matrix and the vector b(c) collects all of the contributions from the boundaries. features the following properties. The matrices G and G(c)
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
283
• G is a singular M-matrix; is a singular Stieltjes matrix. • G(c) are sparse matrices with the same sparse symmetric pattern. The fact that G and Moreover, G and G(c) G(c) have such a structure is important in time discretization because it allows for the construction of economic and robust semi-implicit schemes [14,15].
7. Stronger constraints in the limiting process A discrete version of the analytical properties (2(a)–(b)) stated in Proposition 1 can be recovered when a stronger constraint is assumed in the limiting process. Let us take 1 if vij > 0, (COND.2) 0 σij (c) < ∞ otherwise, in place of (COND.1). Then, it follows that if c 0 and ci = 0 ⇒ Fi (c) 0, where the notation x 0 means that every entry of the vector x is non-negative. This implies a nonnegativity result because c(t) 0 for all t > 0 if c(0) 0. Let us introduce the following mesh-dependent semi-norms |eij ||vij |(ci − cj )2 , |c|22,U = eij ∈EhI
|c|22,R =
2|eij | vij+ σij (c) + vj+i σj i (c) (ci − cj )2 .
eij ∈EhI
The seminorm |c|2,U is related to the upwind (U) numerical diffusion, while the seminorm |c|2,R measures the numerical anti-diffusion related to the reconstruction algorithm (R). A stronger result holds for the FV approximate solution, that can be formally stated as follows. Proposition 3. Under the same assumptions of Proposition 1 and taking (COND.2) satisfied by σij (c) the FV approximate solution c satisfies (a ) 0 ci (t) M, (b ) dtd c1 + γ 1 (c), vn = 0, (c ) dtd c22 + γ 1 (c)2 , vn = B(c) − A(c), where A(c) = |c|22,U − |c|22,R ,
B(c) =
2 γ 0 − γ 1 (c) , vn .
We remark that A(c) can be negative, i.e., |c|22,U < |c|22,R , and that B(c) = O(h2 ).
284
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Clearly, it is a highly desirable feature that A(c) be non-negative and as small as possible in order to reduce the difference A(c) − B(c), that is to control the numerical diffusion by the numerical antidiffusion. Sufficient conditions can be formally stated that are capable of ensuring this feature. Let us substitute the constraint (COND.2) by the stronger one: 1/2 if vij > 0, (COND.3) 0 σij (c) < ∞ otherwise. Then, the following proposition holds. Proposition 4. Under the assumptions of Proposition 1 and (COND.3) on σij (c) the FV approximate solution c satisfies (i) |c|2,U |c|2,R , that is A(c) 0, (ii) A(c) = O(h2 ), on smooth solutions.
8. Numerical tests Test case 1 In the first test case, the domain is Ω = (0, 1) × (0, 1), the velocity field is V = (1, 1)T , the initial condition is C0 (x) = 0 and the inlet boundary condition is given by taking g(x) = 1 on Γ − . The same solution is computed on a regular and an irregular mesh. The two meshes have almost the same numbers of triangles. The regular mesh is generated by splitting Ω in 25 × 25 equally sized squares and then splitting each square in two triangles. The irregular mesh is generated by the mesh generator code TRIANGLE (see Acknowledgements) forcing irregularity while the unstructured solution code was developed by using P2MESH software system [16]. Test case 2 In the second test case, the domain is Ω = (−1, 10) × (− 12 , 12 ), the velocity field is V = (1, 0)T and the inlet boundary condition is given by taking g(x) = 0 on Γ − . The initial condition is given by 1 if x ∈ [a, b] × [c, d], C0 (x) = χ[− 1 , 1 ]×[− 1 , 1 ] (x), χ[a,b]×[c,d] (x) = 2 2 2 2 0 elsewhere, and the exact solution is C(t, x) = C0 (x − tv). Test case 3 In the third test case, we take the same domain, velocity field and inlet boundary condition of the second test case, but we consider a different initial solution. This latter one is 2 C0 (x, y) = 1 − x 2 1 − 4y 2 χ[−1,1]×[− 1 , 1 ] (x). 2 2
In both second and third test cases calculations are performed on a regular mesh of 2200 triangles.
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
285
Fig. 2. Solution of test case 1 at t = 0.5 s (left: regular mesh; right: irregular mesh).
Fig. 3. Solution of test case 2 at t = 0 s and t = 8 s (top: piecewise constant reconstruction; bottom: piecewise linear reconstruction).
Fig. 4. Solution of test case 3 at t = 0 s and t = 8 s (top: piecewise constant reconstruction; bottom: piecewise linear reconstruction).
Fig. 2 shows the piecewise linear reconstructed solution of test case 1 computed at t = 0.5 s. In both the regular and irregular mesh calculations we forced (COND.2) on the reconstruction algorithm and no overshoot or undershoot is produced. Fig. 3 shows the superimposed contours of the computed solutions
286
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Fig. 5. L∞ norm of the solution for test case 2 (left) and 3 (right).
of test case 2 at t = 0 (on the left) and t = 8 s (on the right). As expected the less diffusive solution is obtained with linear reconstruction. As in the first test case, no overshoot or undershoot are produced. Fig. 4 shows the superimposed contours of the computed solutions of test case 3 at t = 0 (on the left) and t = 8 s (on the right). As in the first test case, no overshoot or undershoot are produced. Fig. 5 shows the L∞ norm of the numerical solution of the test cases 2 and 3. The exact solution shows a constant value equal to 1. As expected, the piecewise constant reconstruction produces a more dissipative solution which results in a faster decay of the L∞ norm. Fig. 6 shows the L2 norm of the numerical solution of test cases 2 and 3. The exact solution shows a constant value equal to 1 for test case 2 and about 0.6584 for test case 3. As expected, the piecewise constant reconstruction produces a more dissipative solution which results in a faster decay of the L2 norm. Fig. 7 shows the logarithm of A(c) of the computed solution for the test cases 2 and 3. Ideally for a non-dissipative scheme A(c) ≡ 0 moreover experimentally for this test case is enough condition (COND.2) to have A(c) 0. It is well revealed the effect of linear reconstruction which lower A(c) such a values. Fig. 8 shows the logarithm of A(c) of the numerical solution of test case 3 by using four different meshes. The meshes are composed by 550, 2200, 8800 and 35200 triangles, corresponds to a mesh size resolution h = 0.2, 0.1, 0.05, 0.025. Since the solution of test case 3 is smooth, A(c) = O(h2 ) as can be seen in Fig. 8. Finally, the L2 norm of the difference between the exact and numerical solution is given in Table 1. The estimated order of accuracy of the methods is also displayed in Table 1. Table shows that the measured order is lower than the expected one. This effect is due to the non-linearity introduced by the limiter.
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Fig. 6. L2 norm of the solution for test case 2 (left) and 3 (right).
Fig. 7. A(c) values for test case 2 (left) and 3 (right). Greater values are for the piecewise constant reconstruction.
287
288
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Fig. 8. A(c) values for test case 3. Lower values are for the finer meshes.
Table 1 L2 norm of the error and numerical estimation of convergence order for the test case 3 h
0.2 0.1 0.05 0.025
1st order
2nd order
Eh = c − π(C)2
E log2 ( Eh/2 ) h
Eh = c − π(C)2
3.85 × 10−1 2.93 × 10−1 2.00 × 10−1 1.23 × 10−1
– 0.40 0.55 0.69
2.37 × 10−1 7.07 × 10−2 2.17 × 10−2 6.92 × 10−3
E
log2 ( Eh/2 ) h – 1.74 1.65 1.71
9. Conclusions The numerical dissipation of the first-order upwind scheme and the second-order linearly reconstructed scheme are expressed in terms of suitable mesh-dependent seminorms of the discrete solution field c. The anti-diffusive effect of the reconstruction process is outlined in a dimensionally independent way and is related to the properties (a), (b) and (c) stated in Proposition 1 for the analytical solution. In the case of the linear reconstruction algorithm, the discrete version of these properties does not hold. However, the limiting strategy (COND.2) restores the discrete version of points (a), (b) and (c).
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
289
Acknowledgements The unstructured Delaunay grids were generated by the mesh generator TRIANGLE, a code implemented by Shewchuck, see the URL: http://almond.srv.cs.cmu.edu/afs/cs/project/quake/public/ www/triangle.html.
References [1] S.K. Godunov, Finite difference method for numerical computation of discontinous solution of the equations of fluid dynamics, Mat. Sb. 47. [2] P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21 (1984) 995–1011. [3] A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (3) (1983) 357–393. [4] R.J. LeVeque, Numerical Methods for Conservation Laws, second ed., Birkhäuser, Basel, 1992. [5] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, second ed., Springer, Berlin, 1999, a practical introduction. [6] E. Godlewski, P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer, New York, 1996. [7] C. Hirsch, Numerical Computation of Internal and External Flows, Wiley, Chichester, 1990. [8] A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high-order accurate essentially nonoscillatory schemes. III, J. Comput. Phys. 71 (2) (1987) 231–303. [9] T. Sonar, On families of pointwise optimal finite volume ENO approximations, SIAM J. Numer. Anal. 35 (6) (1998) 2350–2369 (electronic). [10] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1) (1996) 202–228. [11] C.-W. Shu, High order ENO and WENO schemes for computational fluid dynamics, in: T.J. Barth, H. Deconinck (Eds.), High-Order Methods for Computational Physics, Springer, Berlin, 1999, pp. 439–582. [12] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, vol. VII, North-Holland, Amsterdam, 2000, pp. 713–1020. [13] M.E. Hubbard, Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids, J. Comput. Phys. 155 (1) (1999) 54–74. [14] E. Bertolazzi, A finite volume scheme for two-dimensional chemically reactive hypersonic flow, Internat. J. Numer. Methods Heat Fluid Flow 8 (8) (1998) 888–933. [15] E. Bertolazzi, G. Manzini, A triangle-based unstructured finite-volume method for chemically reactive hypersonic flows, J. Comput. Phys. 166 (1) (2001) 84–115. [16] E. Bertolazzi, G. Manzini, Algorithm 817 p2mesh: Generic object-oriented interface between 2-d unstructured meshes and fem/fvm-based pde solvers, ACM Trans. Math. Software (TOMS) 28 (1) (2002) 101–132.
Limiting strategies for polynomial reconstructions in the finite volume approximation of the linear advection equation Enrico Bertolazzi a,∗ , Gianmarco Manzini b a Dipartimento di Ingegneria Meccanica e Strutturale, Università di Trento, via Mesiano 77, I-38050 Trento, Italy b Istituto di Analisi Numerica IAN-CNR, via Ferrata 1, I-27100 Pavia, Italy
Abstract A cell-centered semi-discrete finite volume method is proposed to accurately solve the time-dependent scalar advection equation. The spatial accuracy is ensured by a piecewise linear reconstruction which requires a suitable limiting strategy to control spurious numerical oscillations. Three different approaches are conceived to limit the approximate solution. 2004 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Finite volume; Non-oscillatory reconstructions; Limiters
1. Introduction The present work details with the numerical approximation of the solution of the linear advection of a scalar quantity C in a constant velocity field V. The construction of an accurate non-oscillatory scheme is presented for the two-dimensional equation; schemes for one and three space dimensions can be similarly derived. First-order accurate schemes are developed in the finite volume framework by using a piecewise constant representation of the approximate solution within mesh cells. Higher-order approximations are usually achieved by implementing polynomial reconstructions. Nonetheless, it is a well-known fact that even the simplest linear schemes are oscillatory near discontinuities. An effective way of controlling the process of formation of spurious numerical oscillations consists in introducing a solution dependent limiter in the reconstruction algorithm. That is, to avoid the limitation of the Godunov theorem [1], appropriate non-linear schemes must also be devised for linear problems. There exists a large literature about non-linear schemes for hyperbolic conservation * Corresponding author.
E-mail address: [email protected] (E. Bertolazzi). 0168-9274/$30.00 2004 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2003.12.007
278
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
laws; we mention, for example, the total variation diminishing (TVD) methods [2–7] and the essentially non-oscillatory (ENO/WENO) methods [8–11]. The aim of this paper is to discuss some limiting strategies for the linear reconstruction in the framework of second-order accurate finite volume approximation methods. To this purpose, we propose and analyze three different limiting constraints for which some important semi-discrete properties are maintained. An important issue that we want also to focus on, is that with proper initial and boundary conditions an initially non-negative scalar solution C must remain non-negative. This property is crucial when C represents a non-negative quantity (i.e., mass density) because negative solutions are physically unacceptable. The model problem is described in Section 2, while the finite volume formulation and the piecewiselinear reconstruction algorithm are detailed in Sections 3 and 4. The first limiting condition is discussed in Section 5; then the method is reformulated in a more compact matrix-like form in Section 6. Two stronger limiting constraints are introduced and their consequences analyzed in Section 7. The effectiveness of our approach is shown by a set of numerical experiments in Section 8, followed by the conclusions in Section 9.
2. The scalar advection equation Our model problem is the two-dimensional scalar advection equation, which reads as Ct + ∇ · (CV) = 0
on Ω,
(1)
where C(t, x) is a scalar quantity advected by an assigned steady velocity field V throughout a polygonal domain Ω. The boundary of Ω, indicated by ∂Ω, can be split into an inflow and an outflow part defined by inflow boundary:
Γ − = {x ∈ ∂Ω | V · n < 0},
outflow boundary:
Γ + = {x ∈ ∂Ω | V · n 0},
and such that ∂Ω = Γ − ∪ Γ + . The model problem described by Eq. (1) is completed by a suitable initial condition C(0, x) = C0 (x),
x ∈ Ω,
and a Dirichlet boundary condition on the inflow boundary, that is C(t, x) = g(x),
(t, x) ∈ R+ × Γ − .
Proposition 1. Under the assumption g(x) 0,
x ∈ ∂Γ − ,
C0 (x) 0,
x ∈ Ω,
the following properties hold:
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
(a) (b) (c)
0 C(t, x) M, d CL1 (Ω) + γ C, vn = 0, dt d C2L2 (Ω) + (γ C)2 , vn = 0, dt
279
(2)
where
M = max C0 L∞ (Ω) , gL∞ (Γ − ) .
Remark 2. The symbol u, v denotes the familiar scalar product of two scalar functions u and v, i.e., u, v = uv ds, ∂Ω
γ is the trace operator and vn = V · n.
3. Discrete spaces and norms Let Ω be a polygonal domain and {Th } a family of triangulations assumed, as usual in finite element approximation, regular and conformal in the sense specified by Ciarlet for h → 0. The discretization parameter is the maximum diameter of the cells in Th . Any triangulation Th is also supposed weakly acute, i.e., the maximum angle allowed for a triangle is π/2. We shall use the notations and conventions listed above. • Triangles are conventionally labelled by an integer index. The generic triangle Ti has area |Ti | and centroid xi . • The positive 1-dimensional Lebesgue measure intersection of two triangles or of a triangle and the border ∂Ω is an edge. The former case defines an internal edge, the latter case a boundary edge. The internal edge shared by triangles Ti and Tj is denoted by eij . For the sake of notation consistency, a boundary edge is also denoted by the same symbol eij . In this latter case i is the unique triangle the edge belongs to and j a specific boundary identifier like a fictitious external triangle. • eij has 1-dimensional Lebesgue measure (length) |eij |, mid-point xij and unit orthogonal vector nij . This latter is arbitrarily oriented from Ti to Tj if eij is an internal edge, and outward directed if eij is on the boundary. • Eh and Th are, respectively, the set of the mesh edges and triangles. Eh denotes the subset of mesh edges belonging to the boundary ∂Ω. • Th (i) is the index set of the triangles adjacent to Ti and Th (i) the index set of the “fictitious triangles”. Let us now introduce the finite dimensional vector spaces of piecewise polynomial functions Pk Eh = f : f |eij ∈ Pk (R) , Pk (Th ) = f : f |Ti ∈ Pk R2 , whose elements are indicated by bold letters. In particular, fi indicates the value f (Ti ) when f ∈ P0 (Th ) and fij the value f(eij ) when f ∈ P0 (Eh ).
280
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
In P0 (Th ) and P0 (Eh ) we introduce, respectively, the mesh-dependent norm and the mesh-dependent scalar product
1/p p |Ti ||fi | , f, g = |eij |fij gij , fp = eij ∈Eh
T i ∈ Th
which are discrete versions of the Lp norm and the L2 scalar product. Let us also define the two projection maps 1 p 0 π(f )|Ti = πi (f ) = f (x) dx, π : L (Ω) → P (Th ), |T i | Ti 1 p − 0
f (x) dx. π(f )|eij = πij (f ) = π : L (Γ ) → P Eh , |eij | eij
The linear reconstruction procedure on unstructured meshes can be viewed as the mapping R : P0(Th ) → Pk (Th ),
R(f)|Ti = Ri (f),
Ti ∈ T h ,
(3)
where Ri (f)(x) = fi + Gi (f) · (x − xi ),
x ∈ Ti .
(4)
The term Gi (f) denotes the cell-centered reconstructed gradient. It is worth noting that this definition is consistent with the initial cell-averaged approximate values in the sense that cell-averages are preserved by the reconstruction process. Formally, for any f ∈ Lp (Ω) we have 1 R π(f ) (x) dx, Ti ∈ Th , (5) πi (f ) = |T i | Ti
which is equivalent to set π ◦ R = id.
4. The finite-volume formulation The finite-volume approximation of Eq. (1) logically proceeds in two steps. In the first step, Eq. (1) is re-formulated on the triangulated domain by integrating within each control volume Ti and then by applying the Gauss divergence theorem [12]. We have d C dx + CV · nij ds = 0, Ti ∈ Th . (6) dt
Ti
j ∈Th (i)∪Th (i)eij
The second step consists in approximating (6) by the following discrete relation dci + Ghij = 0, Ti ∈ Th , |T i | dt
(7)
j ∈Th (i)∪Th (i)
which holds for every t > 0. In Eq. (7), ci ≈ πi (C) approximates the cell-averaged value of the advected quantity C, while Ghij is the discrete flux integral which approximates both the internal and the boundary
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
281
Fig. 1. Definition of qij .
edge advective flux h Gij ≈ C(·, x)V · nij ds,
j ∈ Th (i) ∪ Th (i).
(8)
eij
Line integrals (8) are computed by applying a single-point quadrature rule to the standard upstream numerical flux [7]. The numerical flux integral takes the form Ghij (c) = |eij |Ri (c)(qij )vij+ + |eij |Rj (c)(qij )vij− ,
j ∈ Th (i)
(9)
at internal edges and the form Ghij (c) = |eij |Ri (c)(qij )vij+ + |eij |g(qij )vij− ,
j ∈ Th (i)
(10)
at boundary edges. In Eqs. (9) and (10) we have set vij = V · nij and vij± = (vij ± |vij |)/2. In the case of an internal edge eij , the quadrature point qij is defined as the intersection of eij and the segment connecting xi to xj , see Fig. 1. When a boundary edge is considered, qij is the edge middle point. The discrete initial solution is given by projecting the analytical initial solution onto the piecewise constant space P0 (Th ), i.e., ci (0) = πi (C0 ),
Ti ∈ T h .
Two discrete trace operators can be defined in P0 (Th ) and P1 (Th ). These are γ 0 : P0 (Th ) → P0 Eh , γij0 (c) = cij , and
γ : P (Th ) → P Eh , 1
1
0
γij1 (c)
=
g(qij ) Ri (c)(qij )
if eij ⊂ Γ − , otherwise.
5. A first condition on the reconstruction In order to control the numerical oscillations, we assume that a limiter is introduced into the reconstruction procedure [13]. The limited reconstructed solution within Ti is again denoted by Ri . A first condition of the reconstructed solution is that the sign of the slope in the direction qij − xi is the same of the jump cj − ci (or γij1 (c) − ci for the boundary); thus we assume that sign Ri (c)(qij ) − ci = sign(cj − ci ),
282
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
for internal edges, while sign Ri (c)(qij ) − ci = sign γij1 (c) − ci , for boundary edges. This condition implies that there exists a non-negative scalar coefficient σij (c) such that: cj − ci , if j ∈ Th (i), (COND.1) Ri (c)(qij ) = ci + σij (c) × 1 γij (c) − ci if j ∈ Th (i). 6. The semi-discrete formulation Assuming condition (COND.1) and introducing the scalar quantity v˜ ij (c) = vij+ σij (c) + vj+i σj i (c), in (9) the internal edge numerical flux integral becomes Ghij (c) = |eij |vij+ ci − |eij |vj+i cj + |eij |v˜ij (c)(cj − ci ),
j ∈ Th (i).
The internal flux balance of the cell Ti can be written in a matrix–vector product form Ghij (c) = Gc − G(c)c |i , Ti ∈ T h , j ∈Th (i)
by introducing + if j ∈ Th (i), −|eij |vj i + G|ij = |eik |vik if i = j, k∈Th (i) 0 otherwise, if j ∈ Th (i), −|eij |v˜j i (c) G|ij (c) = |e |v˜ (c) if i = j, k∈Th (i) ik ik 0 otherwise. Similarly, condition (COND.1) in (10) yields the boundary edge numerical flux integral in the form Ghij (c) = |eij |vij γij1 (c) + |eij |vj+i 1 − σij (c) ci − γij1 (c) , j ∈ Th (i), and its the contribution to the numerical flux balance of Ti takes the more compact vector form Ghij (c) = bi (c), Ti ∈ Th . j ∈Th (i)
Consequently, the FV method can be re-formulated as dc = T−1 G(c)c − Gc − b(c) , dt where Tij = δij |Ti | is the mass matrix and the vector b(c) collects all of the contributions from the boundaries. features the following properties. The matrices G and G(c)
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
283
• G is a singular M-matrix; is a singular Stieltjes matrix. • G(c) are sparse matrices with the same sparse symmetric pattern. The fact that G and Moreover, G and G(c) G(c) have such a structure is important in time discretization because it allows for the construction of economic and robust semi-implicit schemes [14,15].
7. Stronger constraints in the limiting process A discrete version of the analytical properties (2(a)–(b)) stated in Proposition 1 can be recovered when a stronger constraint is assumed in the limiting process. Let us take 1 if vij > 0, (COND.2) 0 σij (c) < ∞ otherwise, in place of (COND.1). Then, it follows that if c 0 and ci = 0 ⇒ Fi (c) 0, where the notation x 0 means that every entry of the vector x is non-negative. This implies a nonnegativity result because c(t) 0 for all t > 0 if c(0) 0. Let us introduce the following mesh-dependent semi-norms |eij ||vij |(ci − cj )2 , |c|22,U = eij ∈EhI
|c|22,R =
2|eij | vij+ σij (c) + vj+i σj i (c) (ci − cj )2 .
eij ∈EhI
The seminorm |c|2,U is related to the upwind (U) numerical diffusion, while the seminorm |c|2,R measures the numerical anti-diffusion related to the reconstruction algorithm (R). A stronger result holds for the FV approximate solution, that can be formally stated as follows. Proposition 3. Under the same assumptions of Proposition 1 and taking (COND.2) satisfied by σij (c) the FV approximate solution c satisfies (a ) 0 ci (t) M, (b ) dtd c1 + γ 1 (c), vn = 0, (c ) dtd c22 + γ 1 (c)2 , vn = B(c) − A(c), where A(c) = |c|22,U − |c|22,R ,
B(c) =
2 γ 0 − γ 1 (c) , vn .
We remark that A(c) can be negative, i.e., |c|22,U < |c|22,R , and that B(c) = O(h2 ).
284
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Clearly, it is a highly desirable feature that A(c) be non-negative and as small as possible in order to reduce the difference A(c) − B(c), that is to control the numerical diffusion by the numerical antidiffusion. Sufficient conditions can be formally stated that are capable of ensuring this feature. Let us substitute the constraint (COND.2) by the stronger one: 1/2 if vij > 0, (COND.3) 0 σij (c) < ∞ otherwise. Then, the following proposition holds. Proposition 4. Under the assumptions of Proposition 1 and (COND.3) on σij (c) the FV approximate solution c satisfies (i) |c|2,U |c|2,R , that is A(c) 0, (ii) A(c) = O(h2 ), on smooth solutions.
8. Numerical tests Test case 1 In the first test case, the domain is Ω = (0, 1) × (0, 1), the velocity field is V = (1, 1)T , the initial condition is C0 (x) = 0 and the inlet boundary condition is given by taking g(x) = 1 on Γ − . The same solution is computed on a regular and an irregular mesh. The two meshes have almost the same numbers of triangles. The regular mesh is generated by splitting Ω in 25 × 25 equally sized squares and then splitting each square in two triangles. The irregular mesh is generated by the mesh generator code TRIANGLE (see Acknowledgements) forcing irregularity while the unstructured solution code was developed by using P2MESH software system [16]. Test case 2 In the second test case, the domain is Ω = (−1, 10) × (− 12 , 12 ), the velocity field is V = (1, 0)T and the inlet boundary condition is given by taking g(x) = 0 on Γ − . The initial condition is given by 1 if x ∈ [a, b] × [c, d], C0 (x) = χ[− 1 , 1 ]×[− 1 , 1 ] (x), χ[a,b]×[c,d] (x) = 2 2 2 2 0 elsewhere, and the exact solution is C(t, x) = C0 (x − tv). Test case 3 In the third test case, we take the same domain, velocity field and inlet boundary condition of the second test case, but we consider a different initial solution. This latter one is 2 C0 (x, y) = 1 − x 2 1 − 4y 2 χ[−1,1]×[− 1 , 1 ] (x). 2 2
In both second and third test cases calculations are performed on a regular mesh of 2200 triangles.
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
285
Fig. 2. Solution of test case 1 at t = 0.5 s (left: regular mesh; right: irregular mesh).
Fig. 3. Solution of test case 2 at t = 0 s and t = 8 s (top: piecewise constant reconstruction; bottom: piecewise linear reconstruction).
Fig. 4. Solution of test case 3 at t = 0 s and t = 8 s (top: piecewise constant reconstruction; bottom: piecewise linear reconstruction).
Fig. 2 shows the piecewise linear reconstructed solution of test case 1 computed at t = 0.5 s. In both the regular and irregular mesh calculations we forced (COND.2) on the reconstruction algorithm and no overshoot or undershoot is produced. Fig. 3 shows the superimposed contours of the computed solutions
286
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Fig. 5. L∞ norm of the solution for test case 2 (left) and 3 (right).
of test case 2 at t = 0 (on the left) and t = 8 s (on the right). As expected the less diffusive solution is obtained with linear reconstruction. As in the first test case, no overshoot or undershoot are produced. Fig. 4 shows the superimposed contours of the computed solutions of test case 3 at t = 0 (on the left) and t = 8 s (on the right). As in the first test case, no overshoot or undershoot are produced. Fig. 5 shows the L∞ norm of the numerical solution of the test cases 2 and 3. The exact solution shows a constant value equal to 1. As expected, the piecewise constant reconstruction produces a more dissipative solution which results in a faster decay of the L∞ norm. Fig. 6 shows the L2 norm of the numerical solution of test cases 2 and 3. The exact solution shows a constant value equal to 1 for test case 2 and about 0.6584 for test case 3. As expected, the piecewise constant reconstruction produces a more dissipative solution which results in a faster decay of the L2 norm. Fig. 7 shows the logarithm of A(c) of the computed solution for the test cases 2 and 3. Ideally for a non-dissipative scheme A(c) ≡ 0 moreover experimentally for this test case is enough condition (COND.2) to have A(c) 0. It is well revealed the effect of linear reconstruction which lower A(c) such a values. Fig. 8 shows the logarithm of A(c) of the numerical solution of test case 3 by using four different meshes. The meshes are composed by 550, 2200, 8800 and 35200 triangles, corresponds to a mesh size resolution h = 0.2, 0.1, 0.05, 0.025. Since the solution of test case 3 is smooth, A(c) = O(h2 ) as can be seen in Fig. 8. Finally, the L2 norm of the difference between the exact and numerical solution is given in Table 1. The estimated order of accuracy of the methods is also displayed in Table 1. Table shows that the measured order is lower than the expected one. This effect is due to the non-linearity introduced by the limiter.
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Fig. 6. L2 norm of the solution for test case 2 (left) and 3 (right).
Fig. 7. A(c) values for test case 2 (left) and 3 (right). Greater values are for the piecewise constant reconstruction.
287
288
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
Fig. 8. A(c) values for test case 3. Lower values are for the finer meshes.
Table 1 L2 norm of the error and numerical estimation of convergence order for the test case 3 h
0.2 0.1 0.05 0.025
1st order
2nd order
Eh = c − π(C)2
E log2 ( Eh/2 ) h
Eh = c − π(C)2
3.85 × 10−1 2.93 × 10−1 2.00 × 10−1 1.23 × 10−1
– 0.40 0.55 0.69
2.37 × 10−1 7.07 × 10−2 2.17 × 10−2 6.92 × 10−3
E
log2 ( Eh/2 ) h – 1.74 1.65 1.71
9. Conclusions The numerical dissipation of the first-order upwind scheme and the second-order linearly reconstructed scheme are expressed in terms of suitable mesh-dependent seminorms of the discrete solution field c. The anti-diffusive effect of the reconstruction process is outlined in a dimensionally independent way and is related to the properties (a), (b) and (c) stated in Proposition 1 for the analytical solution. In the case of the linear reconstruction algorithm, the discrete version of these properties does not hold. However, the limiting strategy (COND.2) restores the discrete version of points (a), (b) and (c).
E. Bertolazzi, G. Manzini / Applied Numerical Mathematics 49 (2004) 277–289
289
Acknowledgements The unstructured Delaunay grids were generated by the mesh generator TRIANGLE, a code implemented by Shewchuck, see the URL: http://almond.srv.cs.cmu.edu/afs/cs/project/quake/public/ www/triangle.html.
References [1] S.K. Godunov, Finite difference method for numerical computation of discontinous solution of the equations of fluid dynamics, Mat. Sb. 47. [2] P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21 (1984) 995–1011. [3] A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (3) (1983) 357–393. [4] R.J. LeVeque, Numerical Methods for Conservation Laws, second ed., Birkhäuser, Basel, 1992. [5] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, second ed., Springer, Berlin, 1999, a practical introduction. [6] E. Godlewski, P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer, New York, 1996. [7] C. Hirsch, Numerical Computation of Internal and External Flows, Wiley, Chichester, 1990. [8] A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high-order accurate essentially nonoscillatory schemes. III, J. Comput. Phys. 71 (2) (1987) 231–303. [9] T. Sonar, On families of pointwise optimal finite volume ENO approximations, SIAM J. Numer. Anal. 35 (6) (1998) 2350–2369 (electronic). [10] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1) (1996) 202–228. [11] C.-W. Shu, High order ENO and WENO schemes for computational fluid dynamics, in: T.J. Barth, H. Deconinck (Eds.), High-Order Methods for Computational Physics, Springer, Berlin, 1999, pp. 439–582. [12] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, vol. VII, North-Holland, Amsterdam, 2000, pp. 713–1020. [13] M.E. Hubbard, Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids, J. Comput. Phys. 155 (1) (1999) 54–74. [14] E. Bertolazzi, A finite volume scheme for two-dimensional chemically reactive hypersonic flow, Internat. J. Numer. Methods Heat Fluid Flow 8 (8) (1998) 888–933. [15] E. Bertolazzi, G. Manzini, A triangle-based unstructured finite-volume method for chemically reactive hypersonic flows, J. Comput. Phys. 166 (1) (2001) 84–115. [16] E. Bertolazzi, G. Manzini, Algorithm 817 p2mesh: Generic object-oriented interface between 2-d unstructured meshes and fem/fvm-based pde solvers, ACM Trans. Math. Software (TOMS) 28 (1) (2002) 101–132.