High-Order Central ENO Finite-Volume Scheme with Adaptive Mesh

High-Order Central ENO Finite-Volume Scheme with Adaptive Mesh Refinement for the Advection-Diffusion Equation Lucian Ivan and Clinton P.T. Groth University of Toronto Institute for Aerospace Studies, 4925 Dufferin Street, Toronto, Ontario, M3H 5T6, Canada [email protected], [email protected]

1 Scope High-order methods are being actively pursued in an effort to reduce the cost of large-scale scientific computing applications. Moreover, for complex flows both high-order discretizations and adaptive mesh refinement (AMR) may be required. For hyperbolic conservation laws, the challenge has been to achieve accurate discretizations while coping in a reliable and robust fashion with discontinuities and shocks. For elliptic equations, it is desirable that the discretization procedure remain accurate while satisfying a maximum principle, even on stretched/distorted meshes [Coi94, DABLP99]. In spite of the many successful high-order methods proposed for both structured and unstructured meshes [HEOC87, Abg94, HS99, CS89, Bar93, OV02] and their application to complex engineering problems even in combination with AMR procedures [MJ06, WA06, BC05], there is still no consensus for robust, efficient, and high-order-accurate schemes that fully deal with all of the aforementioned issues and are applicable to more arbitrary meshes. This study considers the development and application of a high-order central essentially non-oscillatory (CENO) finite-volume procedure with AMR to the advection-diffusion equation. A consistent order of discretization error is sought for both hyperbolic and elliptic terms.

2 High-Order CENO Scheme The two-dimensional advection-diffusion equation considered is given by ∂u + ∇ · (V(x, y, u) u) = ∇ · (κ(x, y, u) ∇u) + φ(x, y, u) , !" # !" # !" # ∂t advective(hyperbolic)

dif f usive(elliptic)

(1)

source

where u is the solution (a scalar quantity), V is the advection velocity vector, κ is the diffusion coefficient, φ is a non-linear source term and x and y are the two spatial coordinates. Based on the relative magnitudes of the advective

444

L. Ivan and C.P.T. Groth

and diffusive fluxes, the solutions of this equation can range from those having a more hyperbolic nature and governed by wave propagation phenomena to those having a more elliptic nature and governed by diffusive processes. High-order solutions of (1) are sought here by applying a finite-volume spatial discretization procedure in conjunction with high-order polynomial solution reconstruction, upwind discretization of the hyperbolic flux, and centrally weighting discretization of the elliptic flux. The semi-discrete form of the finite-volume formulation applied to (1) for a cell (i, j) of a twodimensional multi-block mesh composed of quadrilateral computational cells is given by Nf NG d¯ ui,j 1 1 =− (ωF · n Δ)i,j,l,m + φi,j (x, y, u)dv = Ri,j , (2) dt Ai,j Ai,j Vi,j m=1 l=1

where each cell has Nf = 4 faces and a NG -point Gaussian quadrature numerical integration procedure is used to evaluate the solution flux, the sum of the advective and diffusive fluxes F = Fa + Fd = Vu − κ ∇u, through each face. The variable u ¯i,j is the average solution scalar, Ai,j is the cell area, ω is the quadrature weighting coefficient, Δ and n are the length of the cell face and unit vector normal to the cell face, respectively. Analytical or numerical integration can be performed to calculate the integral of the particular nonlinear source term, φi,j , for cell (i, j). In this work, either a two- or four-stage standard Runge-Kutta scheme is used to integrate the system of ordinary differential equations given by (2), depending on the desired accuracy. The hyperbolic fluxes, Fa · n = u V · n, at each quadrature point are determined from the left and right solution values, ul and ur , as follows:

ul (V · n) if V · n ≥ 0, Fa · n = (3) ur (V · n) if V · n < 0. The solution states, ul and ur , are determined by performing piecewise k-order polynomial solution reconstruction within each computational cell. Herein, the k-order CENO reconstruction proposed by Ivan and Groth is used [IG07]. The high-order central ENO (CENO) method is not based on either selecting or weighting reconstructions from multiple stencils. Instead, a hybrid solution reconstruction procedure is used that combines the high-order k-exact least-squares reconstruction technique of Barth [Bar93] based on a fixed central stencil with a monotonicity preserving limited piecewise linear least-squares reconstruction algorithm [Bar93]. The limited reconstruction procedure is applied to computational cells with under-resolved solution content and the unlimited k-exact reconstruction scheme is used for cells in which the solution is fully resolved. Switching in the hybrid procedure is determined by a solution smoothness indicator. Note that for smooth and hyperbolic problems, a k-order reconstruction produces a k + 1-order accurate spatial discretization. A detailed description of the CENO reconstruction

High-Order Central ENO Finite-Volume Scheme

445

procedure and of the smoothness indicator is given in the paper by Ivan and Groth [IG07]. As previously mentioned, the proposed discretization seeks to obtain a global k-order accurate scheme on arbitrary meshes. This implies the use of a k-order accurate gradient for the diffusive flux evaluation, which here is derived from a k+1-order accurate reconstruction. For this reason, even if a quartic (k = 4) reconstruction has the potential to generate a 5th-order discretization of the hyperbolic flux in combination with at least three quadrature points (NG = 3), it will only be 4th-order accurate for the elliptic flux discretization. Therefore, two quadrature points (NG = 2) are sufficient to obtain a global 4th-order accurate scheme with a piecewise quartic reconstruction. In a similar manner to hyperbolic fluxes, numerical elliptic fluxes, Fd · n = −κ ∇u · n , must be evaluated at each quadrature point. Having determined the left and right solution reconstructions, ukl (r) and ukr (r), the solution gradient at the inter-cellular face is obtained as the arithmetic mean of the reconstruction gradients and thus, the flux at the calculation point, r, is evaluated as * 1 ) k k ∇ul (r) + ∇ur (r) · n . Fd · n = −κ (4) 2 The accuracy of (4) can be easily observed to be k-order accurate. To infer other properties such as positivity (i.e. local satisfaction of a discrete maximum principle) or non-existence of odd-even solution decoupling, it is convenient to apply the elliptic discretization to the Laplace operator and analyse the influence coefficient of each entry in the supporting stencil [Coi94]. Note that for a given discretization the influence coefficients depend only on the mesh geometry and not on the actual solution. As proposed by Coirier [Coi94], the positivity of the scheme can be characterized in terms of α0 and α ˜ min coefficients. Ideally, α0 < 0 for stability and α ˜min = 0 for positivity [DABLP99]. In the current work, different mesh geometries were analysed including Cartesian, stretched, and randomly disturbed quadrilateral grids. The analysis has shown that odd-even solution decoupling does not occur. In terms of the stability and positivity, it was found that α0 < 0 but α ˜ min < 0 for discretizations of all order, unfortunately implying that, while stable, none of the discretizations satisfy a discrete maximum principle. Note that for square Cartesian meshes, values for α ˜ min are found to be -0.823 for k = 2, -0.362 for k = 3 and -0.854 for k = 4 when inverse distance geometric weighting is used in the k-exact reconstruction. However, the positivity can be improved by using an inverse distance squared geometric weighting, for which α ˜ min was found to be -0.051 for k = 2, -0.247 for k = 3 and -0.324 for k = 4. For non-Cartesian meshes, large variations in the value of α ˜ min are possible (−5 < α ˜ min < 0), depending on the regularity and topology of the mesh.

446

L. Ivan and C.P.T. Groth

High-order treatment of boundary conditions (BCs) is crucial for developing high-order accurate schemes. Herein, high-order BCs have been imposed by making use of extra rows of ghost cells or by constraining the least-squares reconstruction in cells adjacent to boundaries as described in [OV02]. A flexible block-based hierarchical data structure is used in conjunction with the CENO scheme described above to facilitate automatic solutiondirected mesh adaptation on body-fitted multi-block quadrilateral mesh. In this work, an h-refinement criterion based on the solution smoothness indicator is used to control the refinement of AMR mesh. The AMR CENO algorithm and the refinement criterion are described in details in [IG07].

3 Numerical Results The accuracy of the hybrid CENO algorithm was first assessed based on the circular advection at constant angular velocity of the smooth inflow variation u(x, 0) = ed sin6 (π d) if d ∈ [0 : 1], otherwise 0, where d = x − 0.3. The predicted solution distribution for this problem using the quartic CENO reconstruction on a 80×80 Cartesian mesh is shown in Fig. 1. The error norms for both 4th-order versions (k = 3 and k = 4) of the proposed CENO scheme are also depicted in Fig. 1. As the mesh is refined, the slopes of the L1 -, L2 and L∞ -norms approach in the asymptotic limit -4.53, -4.55 and -4.58 for the cubic reconstruction (k = 3), and -4.94 in all error norms for the quartic reconstruction (k = 4), respectively, indicating that the expected theoretical order of accuracy has been achieved in each case. The high-order CENO scheme was also applied in conjunction with AMR to a problem similar to the previous one. In this case the inflow function was u(x, 0) = e2 d sin6 (2 π d) if d ∈ [0 : 0.8], otherwise 0, where d = x − 0.4. This function exhibits two smooth extrema and a discontinuity close to the second peak. Figure 2 shows the final solution profile along the cross-section

Outflow

Farfield

y

Farfield

1

u 1.63 1.40 1.17 0.95 0.72 0.49 0.26 0.03

10

⎪Ε⎪1, ⎪Ε⎪2, ⎪Ε⎪∞

1.5

0

10

-1

10

-2

10

-3

10

-4

10

-5

10

-6

⎪Ε⎪1 (k=3) ⎪Ε⎪2 (k=3) ⎪Ε⎪∞ (k=3) ⎪Ε⎪1 (k=4) ⎪Ε⎪2 (k=4) ⎪Ε⎪∞ (k=4)

0.5

4.0

5.0

Inflow 0

0

0.5

1

x

1.5

50

100

150

200

250

N1/2

Fig. 1. Predicted 4th-order (k = 4) CENO solution (left) and L1 , L2 and L∞ error norms for k = 3 and k = 4 CENO reconstructions (right)


447

1.5

A

Exact Soln u

y

1

A

0.5

0

0.5

1

Distance

1.5

0

0

0.5

1

1.5

x

Fig. 2. Comparison of 4th-order (k = 3) AMR CENO solution on the final mesh and exact solution along A-A (left) and final mesh with 2,911 10×10 blocks (right)

A-A compared against the exact solution and the final multi-block AMR mesh and the location of line A-A. The results clearly show that the proposed AMR scheme in conjunction with the h-refinement criteria based on the smoothness indicator of the hybrid CENO reconstruction is capable of refining both under-resolved (inaccurate) and non-smooth regions of the solution and will not unnecessarily refine resolved solution content. The smooth peaks are all well captured whereas the solution discontinuity is well identified by the smoothness indicator and well resolved by the AMR procedure. The numerical scheme was also tested for solutions to the Laplace equation for the curved boundary domain shown in Fig. 3. Dirichlet BCs were implemented based on the exact solution u(x, y) = eμ x (cos(μ y) + 2 sin(μ y)), μ = 1.5. A 4th-order solution to this problem and the L1 , L2 , and L∞ error norms for cubic and quartic interpolants are also shown in Fig. 3. The slopes of the L1 - and L2 -norms reach in the asymptotic limit -3.86 and -3.85 for k = 3 and -3.86 and -3.81 for k = 4, respectively. Even if the orders of accuracy for these two interpolants are essentially the same and close to the theoretical value, there is about one order of magnitude difference between the absolute error values, demonstrating the benefits of using a quartic interpolant. The proposed scheme is now applied to problems involving different ratios of advection and diffusion terms as determined by the Péclet number (Pe). Solution of (1) with V = (v0 , 0) and κ(x, y) = 0.01, for the BCs shown in Fig. 4 is considered for Pe = 0.1, Pe = 1, and Pe = 10, depending on v0 . Figure 4 shows the numerical solution obtained for Pe = 10 on an 80 × 40 Cartesian mesh and the error norms for the three Péclet numbers. The results show that the errors generated by the quartic polynomial are consistently lower than those of the cubic interpolant by at least one order of magnitude for all cases and obtain the theoretical accuracy in all error norms. Thus, the L1 -norm for k = 4 is -4.02 (Pe = 0.1), -4.30 (Pe = 1.0) and -3.92 (Pe= 10.0), respectively. In the case of cubic interpolant, the L1 -norm is -3.92 (Pe = 0.1), -3.88 (Pe = 1.0)

448

L. Ivan and C.P.T. Groth 8

Quartic Soln (k=4), 40x40

u

6

5

4

0

1

2

3

4

10

⎪Ε⎪1, ⎪Ε⎪2, ⎪Ε⎪∞

y

7

⎪Ε⎪1 (k=3) ⎪Ε⎪2 (k=3) ⎪Ε⎪∞ (k=3) ⎪Ε⎪1 (k=4) ⎪Ε⎪2 (k=4) ⎪Ε⎪∞ (k=4)

101

728.20 671.87 615.55 559.22 502.90 446.57 390.25 333.92 277.60 221.27 168.61 127.40 71.07 14.75 -4.03 -41.58 -97.91 -154.23

0

10

-1

10

-2

3.0

10-3

10-4

10

4.0

-5

50

5

N

x

100

150

200

1/2

Fig. 3. Fourth-order (k = 4) solution to the Laplace equation on a 40 × 40 grid (left) and L1 -, L2 - and L∞ error norms for cubic and quartic interpolants (right) 10-3

0.5

1

1.5

10

10

-6

3.0

10-7

10-8

10

4.0

-9

10-10

2

2.5

100

3

x

N

Pe = 1.0 Cartesian Mesh

⎪Ε⎪1 (k=3) ⎪Ε⎪2 (k=3) ⎪Ε⎪∞ (k=3) ⎪Ε⎪1 (k=4) ⎪Ε⎪2 (k=4) ⎪Ε⎪∞ (k=4)

10-3

⎪Ε⎪1, ⎪Ε⎪2, ⎪Ε⎪∞

10

-4

10-5 10-6 10

3.0

-7

10-8

10

200

300

400 500

1/2


10-1

⎪Ε⎪1, ⎪Ε⎪2, ⎪Ε⎪∞

10

-2

⎪Ε⎪1 (k=3) ⎪Ε⎪2 (k=3) ⎪Ε⎪∞ (k=3) ⎪Ε⎪1 (k=4) ⎪Ε⎪2 (k=4) ⎪Ε⎪∞ (k=4)

-4

10-5

⎪Ε⎪1, ⎪Ε⎪2, ⎪Ε⎪∞

Outflow, ∂u(x,y)/∂y=0

Inflow, u(0,y)=sin( π y)

0.96 0.90 0.84 0.78 0.73 0.67 0.61 0.55 0.49 0.43 0.37 0.31 0.25 0.20 0.14 0.08 0.02

Dirichlet, u(x,0)=0 0


u

Dirichlet, u(x,1)=0

⎪Ε⎪1 (k=3) ⎪Ε⎪2 (k=3) ⎪Ε⎪∞ (k=3) ⎪Ε⎪1 (k=4) ⎪Ε⎪2 (k=4) ⎪Ε⎪∞ (k=4)

-3

10-5 3.0

10-7

10-9 4.0

10

4.0

10-9

-10

10-11 100

N1/2

200

300

400 500

50

100

N

150

200

250

1/2

Fig. 4. Fourth-order solution for a channel flow problem with Pe = 10 (upper-left) and L1 -, L2 - and L∞ error norms obtained with cubic (k = 3) and quartic (k = 4) interpolants for Pe = 0.1 (upper-right), Pe = 1 (lower-left) and Pe = 10 (lower-right)

and -3.53 (Pe= 10.0), respectively. It can be also seen in the error plots that, for the same accuracy level, the cubic interpolant requires almost twice as many cells as the quartic one. Taking into account that both reconstructions have


449

identical stencils, the only extra cost associated with quartic reconstruction is to determine and store five additional derivatives.

4 Concluding Remarks A high-order finite-volume scheme with AMR for the advection-diffusion equation has been proposed with a number of desirable features. Future research will extend the algorithm to the solution of the full Navier-Stokes equations.

References [Abg94] [BC05] [Bar93] [CS89] [Coi94] [DABLP99] [HEOC87] [HS99] [IG07] [MJ06] [OV02] [WA06]

Abgrall, R.: J. Comput. Phys. 114, 45–58 (1994) Barad, M., Colella, P.: J. Comput. Phys. 209, 1–18 (2005) Barth, T.J.: AIAA, Paper 93-0668 (1993) Cockburn, B., Shu, C.-W.: Math. Comput. 52, 411 (1989) Coirier, W.J.: PhD thesis, University of Michigan (1994) Delanaye, M., Aftomis, M.J., Berger, M.J., Liu, Y., Pulliam, T.H.: AIAA, Paper 99-0777 (1999) Harten, A., Enquist, B., Osher, S., Chakravarthy, S.R.: J. Comput. Phys. 71, 231–303 (1987) Hu, C., Shu, C.-W.: J. Comput. Phys. 150, 97–127 (1999) Ivan, L., Groth, C.P.T.: AIAA, Paper 2007-4323 (2007) May, G., Jameson, A.: AIAA, Paper 2006-304 (2006) Ollivier-Gooch, C.F., Van Altena, M.: J. Comput. Phys. 181, 729 (2002) Wolf, W.R., Azevedo, J.: AIAAJ, 44(10), 2295–2310 (2006)