Higher-Order Accurate Finite Volume Discretization of the Three-Dimensional Poisson Equation Based on An Equation Error Method

Efficient higher-order accurate finite volume schemes are developed for the threedimensional Poisson’s equation based on optimizations of an equation error expansion on local control volumes. A weighted quadrature of local compact fluxes and the flux integral form of the equation are utilized to formulate the local equation error expansions. Efficient quadrature weights for the schemes are then determined through a minimization of the error expansion for higher-order accurate discretizations of the equation. Consequently, the leading numerical viscosity coefficients are more accurately and completely determined to optimize the weight parameters for uniform higher-order convergence suitable for effective numerical modeling of physical phenomena. Effectiveness of the schemes are evaluated through the solution of the associated eigenvalue problem. Numerical results and analysis of the schemes demonstrate the effectiveness of the methodology.


Introduction
The problems of interest in this work have been extensively investigated and we describe some of the relevant research. Higher order compact schemes for elliptic equations have been well-investigated, [5,23,[25][26][27][28], since they achieve high-order accuracies without significant increases in the bandwidths of the coefficient matrices. In [19,25,27], and other application problems [6,17,31], the univariate Taylor series expansion is used to derive the finite difference approximations of individual terms of the differential equation and then coupled to obtain the schemes for multiple spatial dimensions. Subsequently, truncation errors are formulated to assess the accuracies of the schemes. In [12,13], the multivariate Taylor expansion is used in developing higher-order finite volume schemes for elliptic equations in two spatial dimensions. First, local approximations for the unknown and the source are utilized to formulate an equation error expansion for the integral form of the equation. Generalized weighted quadratures of the local approximations are utilized for the equation error expansions based on flux integral formulations of the equation in order to capture all local compact fluxes and preserve operator properties of the equation within the computational domain. The weight parameters are subsequently determined to minimize this wellbalanced error by eliminating the leading terms of the error. As a result, the leading coefficients of the resulting residuals which are the numerical viscosity coefficients are consequently more accurately and completely determined for the computational domain. These viscosity coefficients which describe the growth rates of local residual errors are then optimized by right choices of quadrature weights for the schemes to ensure uniformly diminishing residual errors. In [11], the approach was extended to develop a space-time finite volume differencing framework for effective higher-order accurate discretizations for parabolic equations.
In this paper, we utilize the local equation error expansion approach [11][12][13] to develop higher-order accurate discretizations for the three-dimensional Poisson's equation based on approximations of the flux integral formulations of the equation. This ensures that the resulting residual errors are more accurately and completely determined for optimal quadrature weights for the discrete equations and further guarantee uniformly converging residual errors. That is, the focus is on the efficient representation of both the solution and the source term in the discrete representations of the equations. Our work in [11,12] shows that the rates of higher-order accurate convergence of the solution depends on efficient discretization of the source term in the discrete representation of the equation. While the 19-point stencil and the 27point stencil may be used for the solution [14,25,27], the associated stencils to be used for the source term should be emphasized as either the 7-point stencil or the 19-point stencil or the 27-point stencil. By using the equation error expansion which combines all local compact fluxes of the solution and their corresponding source term approximations in describing the resulting residual errors, the errors associated with all the various stencils for the source term may be optimally assessed. The weighted quadrature descriptions express the approximation of the divergence of the flow about each point within the computational domain and offers the right framework to allow for effective local flux improvements to ensure more uniform higher-order accuracies. Hence, there is effective representation of fluxes to neighboring grid points which ensures conservation of flow within the computational domain.
The equation error expansions allow for flexible configurations of local grid-point clouds [7,10,16] to be adopted as desired in order to effectively approximate fluxes in all compact directions of neighboring grid points needed to account for high frequencies in data distributions [15,30]. The structured distribution of grid points [1] ensures that accuracy estimates of the local residual error can be guaranteed. Thus, a local equation error [2] is formulated for the integral formulation of the equation on the control volume instead of just the coordinate directions as in traditional finite difference formulations [19,25]. This error is then locally constrained by the derivatives of the equation through the Cauchy-Kovalevskaya procedure [20] and then minimized by eliminating the leading terms iteratively to determine the weights to approximate the equation. These weights are further optimized to control the growth of residual errors and ensure uniformly diminishing computational errors and robust higher-order accurate convergence of the schemes.
The paper is organized as follows: In Section 2, we present the discretization framework of the method for a general elliptic equation in flux divergence form in 2 . In Section 3, we apply the method to develop new efficient higher-order schemes for the three-dimensional Poisson equation. We discuss quadrature-weight optimizations of the residual errors for consistent fourth-order convergence for the resulting schemes in Section 4. In Section 5, we apply the resulting schemes to discuss the solution of the associated eigenvalue problem for the unit cube to demonstrate the efficiency of the numerical schemes. Numerical errors for the eigenvalue problem associated with different local stencil supports for the discretizations are demonstrated in Section 6. Sixth-order accurate discretizations for the three-dimensional Poisson equation is discussed in Section 7 and conclusions are presented in Section 8.

Finite Volume Differencing Discretization Framework
We describe the finite volume differencing discretization method for the elliptic boundary value problem where Ω is a bounded domain in 2 or 3 with a smooth boundary Γ. We assume that κ ∈ L ∞ (Ω) is positive and the source function q ∈ L 2 (Ω).
To develop a higher-order accurate discretization for (2.1) with a robust computational convergence, the discretization framework must effectively be able to represent all local fluxes to as many neighboring grid points as possible within the computational domain. That is, the framework must be conservative of all local compact fluxes [4] required for consistent and robust higher-order accuracies. We therefore formulate the equation (2.1) over local control volumes which can support all possible local fluxes to neighboring mesh points rather than independently in univariate coordinate directions [24] as in traditional finite difference formulations. We first write the integral formulation of the equation (2.1) as where dv = dxdy in 2 or dv = dxdydz in 3 . By the divergence theorem, the equation (2.3) is rewritten into the flux integral balance form as (2.4) where ν is the unit outward normal to the boundary S of the domain Ω. Now, consider the domain Ω partitioned into control volumes where each control volume is identified by its centroid mesh point and a distribution of neighboring mesh points. A uniform distribution of grid points is utilized for this work but the approach is applicable for a non-uniform distribution as well. A combination of uniform distribution for regular grid points and non-uniform distribution for irregular grid-points may be adopted [6]. The grid-point clouds for neighboring control volumes overlap [29]  where νh is the unit outward normal vector to Sh which is the boundary of Qh. The equation (2.5) represents the conservation of u about the centroid grid point of the control volume such that variations in the local source term distribution within the control volume are compensated for by the radial fluxes through the boundary Sh [9]. That is, the distribution of u within the control volume Qh is completely determined where ni is number of quadrature points, wi is the collocation weight for the local directional flux ∇u·νhi = (ui−u0) along the radial direction νhi toward location of ui. Based on the adopted set of compact quadrature points, the quadrature approximation of (2.4) about the centroid becomes , (2.8) where κ(u0) is an averaging value of κ(u) about the centroid. The resulting residual is described as (2.10) ∂Qh Qh Clearly, the framework as described in (2.8) allows for regular and non-regular distribution of quadrature points locally about each centroid adaptively. For a more accurate modeling of the local equation error (2.9) using a generalized weighted quadrature description in 2 for instance, we adopt a natural local multivariate Taylor expansion for u about each centroid (x0,y0) by (2.11) where φ is sufficiently smooth and locally defined everywhere such that Consequently, we define the local source term expansion by (2.12) where Δκ is the local differential operator description of (2.1) with unique characteristics of the equation [21] such that In this way, any grid functions of φ and the source term f about each centroid may be determined and utilized to describe the approximations the integral fluxes in (2.4) in the form of (2.8). Thus, the grid point spacings may not necessarily need to be uniform and can be adaptively determined locally. As a result, any desired quadrature points about the centroid may be utilized or included in the approximation of the flux integrals to discretize (2.4) locally. The radial fluxes towards the neighboring quadrature points describe their relative dependencies locally on the total flux within the computational domain and therefore weighted accordingly.
To enable effective higher-order accuracies, the Cauchy-Kovalevskayaprocedure [20]   where ni is the desired number of quadrature points, and wi is the quadrature weight for the local directional flux (φi − φ0). The number of quadrature points forming the desired local distribution is part of the available degrees of freedom, which may be increased to improve local accuracy by incorporating sub-grid-scale points and nonlocal fluxes [3]. To obtain the specifics the discretization for (2.1), the discrete minimax approach is utilized to determine the quadrature weights to annihilate the leading terms of the error expansion. That is, the weights wi and vi are determined to annihilate the leading terms of the multivariate error expansion of (2.5) and to further regulate the growth of the residual error through efficient collocations of the local source terms fi and the solutions φi. The eventual order p of the approximation depends on the degrees of freedom available to discretize the equation (2.4) locally.
One advantage here is that for various innovative ways [22,30] to incorporate local micro scale properties into the numerical model, our comprehensive approach is naturally efficient in determining the right sampling and collocations of the source required for effective and robust higher-order accuracy.

Finite Volume Discretization for Poisson's Equation in 3D
We formulate the equation error expansion (2.9) about each grid point using local compact cloud of 27 grid points within the three-dimensional computational domain. The corresponding control volume Qh, for a uniformly partitioned computational domain is illustrated in Figure 3.1 where the centroid grid-point function φ0, is surrounded by 26 mesh points distributed at three radial distances of and from the centroid.

Fig. 3.1. A compact cloud of neighboring grid-point distribution about each Centroid
As discussed above, we use the multivariate Taylor expansion φ(x0 + α,y0 + τ,z0 + ν) defined in the form of  where U2 = 1 − w1 − 4w7 − 4w19 By annihilating the leading coefficients U2, F2 and U4 of the error (3.2), the quadrature weights for the stiffness matrix H of the discretization are transformed as functions of w19 by where w0 is the weight of centroid value φ0. Consequently, the quadrature weights for the mass matrix Q are transformed as functions of β7 and β19 by 1 β0 = + 12β7 + 16β19 6048 12 As indicated in [13], diagonal dominance of H requires W19 to be chosen to satisfy (3.6) for which the solution is found to be . (3.7) Clearly, there are several possible ways to choose {w7,w19,β7,β19} for a fourthorder accuracy. In the next section, we discuss the optimization of such choices for higher-order convergence robustness of the discretization.

Residual Error Optimization for Robust Higher-order Accurate Convergence
For a uniform distribution of quadrature grid points as illustrated in Figure 3 where p = 1,2,3,··· and i + j + k = p. Notice that the coefficients of higher-order derivatives of φ and f have same signs. For the resulting discretization to produce robust higher-order computational convergence rates, the computational errors must diminish monotonically. After incorporating higher-order derivatives of the equation, the residual error terms assume the form of (4.2). Therefore, robust higher-order computational convergence requires that the coefficients of the higher-order derivatives of both φ and f within the leading terms of the residual error (4.1) maintain similar signs as described in (4.2 to render the leading term of the residual error decreasing monotonically as resolution is refined.

Laplacian Eigenvalue Approximations
To computationally evaluate the higher-order accurate convergence of the resulting schemes, consider the associated eigenvalue problem for the three-dimensional Laplacian −Δu = λu, in ∈ Ω (5.1) u = 0, on ∈ Γ (5.2) where Ω is the unit cube [0 1] × [0 1] × [0 1]. As discussed above, the equation error for (5.1) about the centroid of each control volume Qh ∈ Ω as described in (2.6) is The finite dimensional representation of (5.3) by utilizing the local uniform cloud of quadrature points on the control volume shown in Figure 3.1 may be described as where Φ is a column vector of the discrete representation of the eigenfunctions on Ω which are locally supported by the quadrature points on the control volumes and Λ as the matrix of associated eigenvalue estimates. The mass and stiffness for the discretizations are Q and H respectively, and E is a column vector of ones. As described by the discrete approximation (5.4), the robustness of the resulting numerical schemes for higher-order accuracy depends on the convergence properties of the local residual error Rφ0 in convergence analyses. The eigenvalue estimates, Λ and the associated eigenfunctions Φ are a function of the geometry of the domain and the boundary conditions [8]. Therefore, the accuracy of the estimates Φ and Λ are a direct reflection of the level of effectiveness of the discretization method in sufficiently approximating the equation on the geometry of the domain [18]. Thus, the effectiveness of the discretization for higher-order accuracies may be evaluated through the convergence properties of the local residual error RQ0 as described in (3.5).
The parameters w19, β7, and β19 are free for a fourth-order accuracy and so various strategies may be adopted to further optimize the discretization for convergence robustness. Since the local residual error RQ0 serves as the local source term for the discrete error equation associated with the discretization, our strategy for robustness is to effectively ensure that RQ0 is more uniformly diminishing by choosing the parameters to avoid fluctuations about zero. The computational profiles of various choices of w19, β7, and β19 for higherorder accurate convergence are discussed for the solution of (5.4) by the eigs function in Matlab.

Numerical Errors for Different Local Stencil Supports
In this Section, we discuss numerical results for the eigenvalue estimates associated with the different choices of the local computational supports for the stiffness and mass matrix pair in the discrete equation (5.4). That is, we examine the computations errors in estimating the eigenvalues of the unit cube associated with the possible choices of the parameters w19, β7, and β19 in (5.4).  [19,25,27], is recovered with w19 = 0, β7 = 0, and β19 = 0 in (3.5) and uses a 19×7 point local support for H×Q and the associated local residual error is given as Clearly, condition has a leading error expansion term with same signs for the higher-order derivatives as described in (4.2) ensuring a monotonic convergence than for the 7−point stencil mass matrix with β7 = 0. Furthermore, the mass matrix Q with provides an improved uniform fourth-order accuracies than with as demonstrated between Figures 6.2 and 6.3. With , the mixed partial derivatives of the source term f is eliminated from the O(h 4 ) term and the size of the leading error term gets polluted by O(h 6 ) term for some eigenfunctions.
A convergence analysis through resolution refinement for the 19−point stencil mass matrix Q with is shown from Figure 6 = 0 with = 0 with 6.3. Twenty-Seven Point Stencil. Next, we examine the case for a 27−point stencil local support for the discretization. Since the parameters w19, β7, and β19 are free for a fourth-order accuracy as indicated by the residual error (3.5), our strategy for robustness and consistent higher-order accurate convergence is to ensure that the residual error is uniformly diminishing by reducing the size and fluctuations of the residual error as much as possible as resolution is refined. Therefore, we exercise the choices for w19, β7, and β19 to control leading coefficients in both the O(h 4 ) and O(h 6 ) terms to minimize associated local error pollution. First, we utilize w19 to eliminate ∂x2 ∂ ∂y 6f 2∂z2 in the 6  As evident from Figure 6.5, resolution refinement fails to improve on local accuracies for larger eigenvalues with the 19 × 19 discretization. On the other hand, the 27 × 27 discretization shows much improvement in local accuracy as illustrated in in Figure 6.6. Clearly, the 27−point compact local support for both the stiffness and mass matrices provide for a uniform higher-order accurate convergence rates for the eigenvalue calculations as well as consistent local improvements in accuracies as resolution is refined which is suitable for effective numerical modeling of physical phenomena.

Sixth-Order Accuracy
Since the O(h 4 ) term in the local residual error only involves derivatives of the source term, the right sampling of the numerical source term at sub-grid scale points on the computational stencil may be utilized to achieve sixthorder accuracy. If the fourth-order derivatives of the source term are readily available then the sixth-order accuracy may be achieved by determining β7 to annihilate the O(h 4 ) term in the residual error by .
(7.1) On the other hand, since F fxxxxxx + fyyyyyy + fzzzzzz 6  and then adding (7.2) to the discretization such that Thus, the sixth-order discrete representation of the equations may be described as HΦ = QF − FQh, (7.5) where the stiffness matrix H is composed by (6.2) and the mass matrix is composed by (6.3) with and FQh given in (7.2). Similar sixth-order computational experiments have been demonstrated in [13].

Conclusion
We have demonstrated the effectiveness of the finite volume discretization approach in developing efficient higher-order accurate schemes for the three-dimensional Poisson's equation based on the equation error method. Using a more balanced integral form of the equation to formulate the local error expansion provides for an efficient framework for a more complete local flux modeling on the computational domain that ensures local accuracies as well as uniform higher-order accurate convergence of the resulting schemes. In particular, using the 27−point stencil local support for both the mass and stiffness matrices in the discrete equation achieves robust higher-order accurate convergence, and very suitable for large scale eigenvalue problem in 3D due to local accuracy improvements for finer resolutions.