# Tensor formulation of 3-D mimetic finite differences and applications to elliptic problems.

Abstract. The mimetic discretization of a boundary value problem (BVP) seeks to reproduce the same underlying properties that are satisfied by the continuous solution. In particular, the Castillo-Grone mimetic finite difference gradient and divergence fulfill a discrete version of the integration-by-parts theorem on 1-D staggered grids. For the approximation to this integral principle, a boundary flux operator is introduced that also intervenes with the discretization of the given BVP. In this work, we present a tensor formulation of these three mimetic operators on three-dimensional rectangular grids. These operators are used in the formulation of new mimetic schemes for second-order elliptic equations under general Robin boundary conditions. We formally discuss the consistency of these numerical schemes in the case of second-order discretizations and also bound the eigenvalue spectrum of the corresponding linear system. This analysis guarantees the non-singularity of the associated system matrix for a wide range of model parameters and proves the convergence of the proposed mimetic discretizations. In addition, we easily construct fourth-order accurate mimetic operators and extend these discretizations to rectangular grids with a local refinement in any direction. Both of these numerical capabilities are inherited from the original tensor formulation. As a numerical assessment, we solve a boundary-layer test problem with increasing difficulty as a sensitivity parameter is gradually adjusted. Results on uniform grids show optimal convergence rates while the solutions computed after a smooth grid clustering exhibit a significant gain in accuracy for the same number of grid cells.Key words. mimetic finite differences, tensor products, locally refined grids, elliptic equations

AMS subject classifications. 65H17, 65N06, 40A30

1. Introduction. Mimetic finite difference (MFD) approximations to differential vector operators along with compatible inner products satisfy a discrete analog of the Gauss divergence theorem. On nodal grids, an important MFD family comprises the summation-by-part (SBP) difference operators, whose development and initial applications were focused on the conservative discretization of wave propagation problems [22, 23, 24, 30]. The SBP discretization with suitable boundary treatments mimics the energy dissipation properties of the continuous wave model, and therefore numerical methods are provably stable. This attractive property has encouraged researchers to apply SBP on more general wave propagation problems with complex geometries by considering structured curvilinear grids [2, 10, 21]. By construction, SBP discretization operators are less accurate with respect to differentiation at the boundary nodes than for the interior mesh points. On staggered grids, we find two established MFD discretization methods. The first staggered MFD family is based on the discrete support operators (SO), the Divergence D, and the Gradient G that preserve the conservative property of negative adjointness satisfied by the continuous counterparts [18]. The SO method has been extensively applied to diffusion and hyperbolic problems even on non-uniform meshes [11, 15, 17, 19]. The operators that were reformulated in [15] on polygonal meshes provide second-order accuracy on the whole discretization grid, but in several of the previous SO applications, a first-order approximation of the boundary fluxes is performed. The second staggered family presents discrete MFD D and G operators that exhibit second-, fourth-, and even sixth-order accuracy for all grid locations including boundaries. These conservative operators are introduced in [5], and later Castillo and Grone presented an original algebraic construction of D and G in [4]. Recently, this methodology has been improved in [26] by using linear constrained optimization to enforce the mimetic conditions, and new 1-D operators with eighth-order accuracy and higher have been obtained. For a numerical treatment of Neumann and Robin boundary conditions, the authors in [7] developed a new operator B to approximate the boundary fluxes with the same accuracy as for the interior grid using the discrete Green identity as a construction basis. This Castillo-Grone MFD set of D, G, and B operators are used in several applications for 2-D diffusion and wave propagation problems as we detail below, however, its formulation and use on 3-D meshes is still under development. This paper contributes to this field.

The sound evolution of SO and Castillo-Grone MFD during the last decades has led to a mature status among the numerical methods for partial differential equations. In fact, three new survey books have been published recently, cf. [6, 8, 20], and [6] is fully dedicated to 1-D and 2-D Castillo-Grone operators. These textbooks provide a complete and updated review of construction methodologies, application areas, and related approaches up to the year 2014. Tensor formulations of higher-dimensional mimetic operators are discussed in [6, 20] for rectangular or logically rectangular grids. In a two-dimensional context, some authors have developed second-order mimetic schemes for elliptic problems [13, 14] based on the 1-D discretization proposed in [7], but none of them explicitly mention the use of tensor constructions. A tensor-based mimetic scheme for solving transient Maxwell equations under Dirichlet conditions is introduced in [25] on 2-D rectangular meshes. Recently, this formulation has been used in [27] as a basis for the object-oriented mimetic toolkit, which is an ongoing software project with current applications to 2-D Poisson equations. An initial formulation of three-dimensional mimetic schemes for elliptic problems based on Castillo-Grone operators without an explicit tensor formulation has been reported in [1]. However, the numerical results for several tests there are reduced to the case of 2-D logically rectangular grids. Similarly, these MFD operators have been expressed in a 3-D tensor-like notation in [28] for the modeling of acoustic wave propagation, but the numerical implementations are limited to 2-D Cartesian meshes. To the best of our knowledge, a complete 3-D tensor formulation of Castillo-Grone mimetic operators accompanied by satisfactory numerical testing has not yet been reported in the literature.

Tensor-based extensions of finite difference (FD) operators to higher spatial dimensions from their standard 1-D formulas is a well-known technique [9, 29]. Nonetheless, a similar application to the Castillo-Grone MFD operators has not been fully achieved in three dimensions as we have mentioned above. The main goal of this article is to provide a consistent 3-D formulation of these operators starting from the second-order mimetic discretizations described in [7, 14], which are unique on 1-D and 2-D equidistributed grids. These new tensor formulations can be easily used to construct 3-D versions of fourth- or higher-order accurate discretizations based on the 1-D operators given in [4, 5] with similar accuracy. Next, we apply these new mimetic operators to the discretization of a special family of elliptic problems. The convergence of our new schemes is formally studied and later observed for various numerical experiments. It is well known that high-order mimetic operators are not unique, and their performance is parameter dependent, but here we only use some particular cases to test our implementations.

The content of this article has been divided into seven sections. The first section is this introduction with our review of related works. In Section 2 we present a new tensor mimetic discretization of the continuous three-dimensional gradient and divergence by using Kronecker products of the 1-D operators described in [7, 14]. Next, we use these tensor formulations in Section 3 to obtain a novel second-order accurate mimetic discretization of a 3-D family of elliptic equations and develop a consistency analysis of this discretization. Section 4 presents a simple eigenvalue analysis to justify the solvability and stability of this discretization. In Section 5, we discuss the straightforward generalization of mimetic operators to fourth-order accuracy and grids with variable spacing. In Section 6, we apply the new mimetic scheme to a synthetic test problem with an adjustable boundary layer and evaluate the results in terms of numerical convergence. Finally, we discuss the conclusions from this work in Section 7.

2. Tensor formulation of mimetic D, G, and B. In a multidimensional domain [OMEGA], a generalized version of the Gauss divergence theorem can be written for differentiable fields v and u as

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

where n represents the continuous outer normal vector to the domain's boundary surface [partial derivative][OMEGA]. According to [4, 7], a FD method for a boundary value problem (BVP) is called mimetic if the discretizations of the continuous gradient[??], the divergence [??]., and the boundary flux [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is performed by discrete operators G, D, and BG, respectively, such that all together satisfy the following approximation to (2.1)

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Here, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], Q and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] P denote weighted inner products defined by square and positive definite matrices Q and P. Notice that the role of B is well suited for BVPs where v = [??]u, but in general contexts, this operator incorporates the projection of v on the outer normal vector at discrete boundary elements. The matrix structure of the operators G, D, and B, depend on a particular discretization of the domain [OMEGA] and the grid distribution of the vector evaluations [??] and [??] of the continuous fields v and u, respectively.

2.1. 1-D operators and local truncation errors. In the 1-D case, the authors in [4, 5, 7, 14] define a uniform partition of the interval [OMEGA] = [0, 1] given by the (N + 1) nodes [x.sub.i] = ih (0 [less than or equal to] i [less than or equal to] N), for h = 1/N, in addition to the cell centers [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The grid location of the components of the vectors [??] and [??] are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The matrices G and D collect a set of staggered differentiation stencils such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a vector of approximations to du/dx at all grid nodes, while [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] estimates dv/dx at every cell center. To make the dimensions consistent in the leftmost vector inner product in (2.2), the row vector [0,..., 0] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is considered as the first and last rows of the matrix D. Thus, the operator D belongs to the matrix space [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] whose non-zero entries at the inner rows correspond to the central stencil for the staggered FD stencil [h.sup.-1] [-1, 1]. This second-order FD stencil also represents the non-zero elements of the inner rows of the operator G with the matrix dimension (N + 1) x (N + 2). Now, the first row in G is simply the lateral stencil [h.sup.-1] [- 8/3, 3, - 1/3] that keeps the second-order accuracy for the differentiation process. An approximation to du/dx at [x.sub.N] can be obtained by the application of this latter stencil after reversing the order and changing the sign of each coefficient. This last backward FD formula defines the last row in G. For these choices of the G and D operators, the weights Q and P are diagonal matrices whose coefficients correspond to the rectangle and Newton-Cotes 3/8-quadrature rules, respectively. That is, Q = h diag(1,..., 1), and P = h diag(3/8, 9/8, 1,..., 1, 9/8, 3/8). The boundary operator B is derived from (2.2) and given by

(2.3) B = QD + [G.sup.T] P.

Thus, B corresponds to the following (N +2)x(N +1) matrix in this second-order formulation

(2.4) B = [[-e.sub.1]; [b.sub.0]; -[b.sub.0]; 0,..., 0; -[b.sub.N]; [b.sub.N]; [e.sub.N]].

In this equation, the non-zero rows of B are the canonical vectors [e.sub.1] and [e.sub.N], and the special rows are [b.sub.0] = (1/8, - 1/8, 0,..., 0) and [b.sub.N] = (0,..., 0, 1/8, - 1/8). The first and last rows of B incorporate the outward normal vector at both boundaries and permit the approximation of [partial derivative]u/[partial derivative]n at these locations similarly by the first and last components of the vector [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The availability of separate operators G and D leads to a natural approximation to the Laplacian operator by DG. However, the mimetic discretization of the Laplacian of the field u is actually given at the cell centers by the interior components of (BG + DG)[??]. For clarity, below we list the whole set of FD stencils comprised into this approximation vector along with the associated local truncation errors

(2.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

(2.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The lateral approximations to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and the central stencil for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with i = 5, 7,..., N - 5/2 correspond to standard second-order Taylor formulas. On the other hand, the one-sided stencils for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] result from adding the special contribution of the non-zero interior entries of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] to the standard Laplacian approximation DG at those grid centers. This process reduces the accuracy to first order, but it is necessary to enforce the conservation property (2.2) for the discrete solutions for second-order differential equations. Previous mimetic applications of this kind have discussed this nominal accuracy deficiency in the discretization of

boundary value problems [7] and two-dimensional elliptic equations [14] but at the same time have shown consistent second-order convergence in numerical tests. First-order local truncation errors are only present at a few near-boundary grid layers, and this does not degrade the global convergence. In this work, we also provide a set of 3-D numerical experiments where the mimetic solutions show global second-order convergence. The local truncation errors given by equation (2.6) are the basis of the consistency analysis of our 3-D mimetic numerical scheme stated in Section 3.

2.2. Extension to 3-D rectangular grids. Generalizations of the above mimetic gradient, divergence, and boundary flux operators to 3-D rectangular meshes can be achieved by a subsequent application of these 1-D operators along each grid coordinate in a natural order. This process leads to an elegant formulation of 3-D mimetic operators as the tensor product of their 1-D counterparts. With no loss of generality, we consider the unit cube [OMEGA] = [0, 1] x [0, 1] x [0, 1] as the discretization domain. Again, we define N equal-sized cells in the x direction, where the grid nodes [x.sub.i] and cell centers [x.sub.i+1/2] are now collected into vectors [X.sub.n] and [X.sub.c], respectively. In addition, we introduce the vector [X.sub.b] = ([x.sub.0], [x.sub.N]) that only contains boundary nodes, and these grid edges joined with the cell centers form the vector [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Let us denote by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the discrete 1-D gradient and divergence operators that perform a mimetic differentiation along an x gridline. Note that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] represents a linear operator that maps values of a function u evaluated at [X.sub.cb] to approximations of [partial derivative]u/[partial derivative]x at [X.sub.n]. Similarly, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] linearly maps evaluations of a function v at the nodes [X.sub.n] to approximations of [partial derivative]v/[partial derivative]x at the mid-cell points [X.sub.c].

To proceed with the staggered gridding of [OMEGA], we next consider uniform partitions of M cells along the y direction and L cells in the z direction as well as coordinate-dependent grid vectors [Y.sub.n], [Y.sub.c], [Y.sub.b], [Y.sub.cb], [Z.sub.n], [Z.sub.c], [Z.sub.b], and [Z.sub.cb], and finally a pair of 1-D gradient and divergence operators on each of these coordinate lines, i.e., [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The grid locations of a given scalar field f and a vector field v = ([v.sub.1], [v.sub.2], [v.sub.3]) are presented in the second column of Table 2.1. The operator x denotes the Cartesian product used in set theory, and its application to this context renders tuples of grid sites. Evaluations of u, [v.sub.1], [v.sub.2], and [v.sub.3] are orderly arranged into vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], respectively, in such a way that the index associated to the x coordinate increases at first, followed by the index of the y coordinate, and then the z-index. That is, the discrete values follow a lexicographical or natural ordering of these computational vectors. Note that this 3-D grid distribution of the u-values results from the tensor product of the three 1-D distributions along each coordinate axis for the same field and incorporates the four grid corners of each grid plane at a constant value z. This simple fact brings a tremendous regularity into the grid geometry that highly simplifies the construction of the discrete 3-D gradient operator as we describe below. This corner-aided grid formulation has not been explored in previous mimetic applications. The conventional mimetic u gridding omits corner nodes as well as all boundary values of u at the bottom and top grid planes with the heights z = 0 and z = 1, respectively. In Table 2.1, the set of u-grid points in the conventional formulation is given in the last column, which excludes a total of 4(N + M + L + 2) corner and boundary points. On dense grids, this extra amount of evaluations is negligible compared to the total number of nodes for u which is O[(N + 2) * (M + 2) * (L + 2)]. Away from these exceptions, the grid cells for both formulations are the same. Figures 2.1 and 2.2 illustrate the generalization of such cells from 1-D to 2-D and then to 3-D dimensions in the case of either an interior cell or a boundary cell. The components of the vector field v are distributed in a staggered way along the grid. Each of these components takes a separate location coincident with the central point of a cell side. Figures 2.1 and 2.2 also depict this grid arrangement for the v components as the grid dimension increases. The grid distribution of the values for [v.sub.1], [v.sub.2], and [v.sub.3] in our formulation replicate the one used in the traditional mimetic framework as shown in Table 2.1.

The staggered grid distribution of the discrete field u leads to a straightforward formulation of the 3-D mimetic gradient [G.sup.3D] in terms of three block components

(2.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Here, I represents the identity matrix with the dimension corresponding to its integer subindex. The operator [cross product] corresponds to the matrix Kronecker product, and we follow the most-right priority ordering when defining the block components of [G.sup.3D]. The product in the first block, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], yields an approximation to the first component [partial derivative]u/[partial derivative]x of the continuous gradient at the center point of every cell face with a constant x grid coordinate. Due to the inclusion of boundary corner evaluations of u in our formulation, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] also produces approximations to this gradient component at the eight boundaries of the grid cube. This feature is not offered by the conventional mimetic formulation, and Table 2.2 compares the image of both discrete first-component gradient operators. Similar observations can be made for the two remaining block components of [G.sup.3D] designed to approximate [partial derivative]u/[partial derivative]y and [partial derivative]u/[partial derivative]z on the corresponding cell faces. Table 2.2 also shows the grid image of these approximations. The availability of all discrete gradient components at the grid boundaries might allow a precise approximation of the generalized fluxes K [??] u for a material-dependent full tensor K even though the staggered location of the individual components demands certain interpolations. This kind of boundary conditions is common in diffusion problems on heterogeneous media where similar interpolations must be required to discretize [??].(K [??] u) at the interior grid points.

We next define a new discrete vector [??] that subsequently collects all values of [??] followed by [??] and finally [??], keeping the natural ordering already imposed on those vectors. Thus, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The mimetic approximation to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be computed at the central points of every grid cells by the product [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where the 3-D divergence operator [D.sup.3D] also follows a block-componentwise structure

(2.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Above, a matrix [??] is derived from an identity matrix after adding zeros to the first and last rows as well as to the first and last columns. For instance,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

in case of the x coordinate. Table 2.2 shows that the grid locations of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are similar to those defined in the conventional formulation. These locations coincide with the interior grid sites of the scalar function u shown in Figures 2.1 and 2.2.

The mimetic approximation to the continuous boundary flux [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where [??] represents the outward normal vector at the boundary grid faces, can be computed by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for [B.sup.3D] given by

(2.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Again, our rectangular gridding leads to a compact definition of the three-dimensional B operator as the tensor product of the individual one-dimensional versions.

3. A mimetic scheme for elliptic problems. The mimetic operators presented in the previous sections can be used to formulate a three-dimensional, second-order, mimetic scheme for elliptic equations. In this process, we consider a general second-order elliptic equation of the form

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

under a Robin boundary condition

(3.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In these equations, k, [omega], F, a, [gamma], and f are given functions, and u is the unknown function. We assume that [gamma], k, and [omega] are bounded from below by a positive constant. In this paper, we also consider the situations when k is a scalar function or a diagonal spatially variable tensor. A complete anisotropic formulation for these diffusion-type equations is outside the scope of this paper, and it is a topic of current research; see [6] (and the additional references therein) for a possible formulation in that direction in a two-dimensional context. Let us assume that the domain [OMEGA] is a Cartesian parallelepiped with n as the normal unitary exterior vector at its boundary [partial derivative][OMEGA]. The combination of equations (3.1) and (3.2) is a well-posed problem whose solution will be approximated by the mimetic scheme proposed in this section.

A formal mimetic discretization of the above elliptic problem can be accomplished by a direct substitution of the continuous differential operators by their three-dimensional mimetic approximations. In the case of equation (3.1), the corresponding mimetic discretization obeys

(3.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In this equation, [??] is the mimetic approximation to the function u at the grid sites given in Table 2.1, [??] is the projection of the function F at the block centers, K is a diagonal matrix representing the values of the function k at the blocks faces, and W is also a diagonal matrix whose entries are the values of the function [omega] at the block centers. In a similar way, the mimetic discretizations of the boundary condition (3.2) is given by

(3.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In this relation, A and [GAMMA] are diagonal matrices with the values of a and [gamma] at the exterior faces of the boundary blocks. At the same boundary locations, the evaluations of f are collected into the vector [??]. Similar to the 1-D mimetic setting, the discrete operator [B.sup.3D] [G.sup.3D] allows an approximation of the directional derivative at the boundary faces. The discretizations (3.3) and (3.4) can be combined by the superposition principle into a single equation

(3.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

which represents the full three-dimensional mimetic scheme associated to the elliptic problem (3.1) and (3.2). Equation (3.5) is very neat, correct, and has the intrinsic computational advantage of avoiding an explicit enumeration of the discrete values of u along the 3-D grid. The grid locations of the unknown u automatically follow a lexicographic order resulting from the tensor product formulation, which greatly simplifies the coding of the linear system (3.5). However, an explicit description of the equations defined by (3.5) is important for a consistency analysis. In the light of the 1-D mimetic approximations given by (2.5), one may expect unconventional formulas for the 3D Laplacian in the case of (3.5), which are different from the standard seven-point FD stencil. Actually, the only non-standard Laplacian formulas are those defined at the center points of the boundary grid blocks as well as similar approximations designed at its adjacent block(s). In other words, equations associated to grid blocks that are two blocks away from the boundary reduce to the standard seven-point Laplacian stencil. By taking this into account, the minimum grid size that allows describing all the possible stencils is a 5 x 5 x 5 grid. By symmetry, the new mimetic scheme only presents six non-standard stencils, which are shown in Figure 3.1. In these stencil graphs, the grid point (i, j, k) corresponds to a block center and therefore does not intervene in the discretization of the boundary conditions (3.2). The stencils in Figure 3.1 are topologically singular because all of them contain at least one branch with a node on the boundary [partial derivative][OMEGA]. Given that we have used integer tuples to identify the stencil's center points, any boundary node will have a half integer among its indices. Here, we refer to a branch as singular, if it has four nodes, one of which is a boundary node, and the other nodes correspond to block centers. Since the general equations in (3.5) are quite involved, we assume below that k, [gamma], and a are constant and equal to one and the grid is uniform with a common grid step h. These simplifications allow us to present the unconventional discrete equations in (3.5) with a minimum of notation while keeping their non-standard features.

The display a) in Figure 3.1 shows a mimetic Laplacian stencil with only one singular branch and corresponds to the following formula

(3.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This stencil represents a typical center node of a grid block that is one block away from a boundary face in only one direction. The stencil in display b) of the same figure shows two singular branches and follows the equation

(3.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This type of Laplacian discretization is associated to a center node whose block is simultaneously one block away from only two adjacent grid boundary faces. Display c) depicts the less standard of these six mimetic stencils with three singular branches with the corresponding formula

(3.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This type of stencils arise at the center node of a grid block displaced by one block from three adjacent grid boundary faces. These blocks are the ones located close to the eight corners of the grid parallelepiped. Notice that all these 3-D stencils in (3.6), (3.7), and (3.8) result from a cross combination of the standard centered FD [h.sup.-2][1, -2, 1] and the 1-D mimetic stencil placed on the third row of (2.5) with a first-order truncation error as given by (2.6). Thus, the accuracy of the Laplacian approximations in (3.6), (3.7), and (3.8) is only of first order, but it could improve to second order in the case of harmonic functions.

Each of the next three stencils in Figure 3.1 (panels d), e), and f)), represent a non-standard Laplacian discretization at the center of a boundary block, i.e., a grid block with at least one face on the boundary [partial derivative][OMEGA]. The stencil in graph d) has one singular branch and corresponds to the approximating equation

(3.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This stencil is applied at the center of a grid block with exactly two boundary faces, which explains the two short branches in its configuration. The stencil in display e) shows two singular branches and satisfies the discrete equation

(3.10) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This stencil corresponds to a grid block with only one face on the grid boundary. The node on that boundary face belongs to the short branch depicted in the stencil graph, and the other branches are singular with the terminal nodes on the boundary faces. The three boundary faces reached by this stencil are orthogonal and assemble a corner of the grid parallelepiped. Finally, the stencil in display f) is also associated to the boundary blocks with just one face on [partial derivative][OMEGA], but it has only one singular branch. This stencil is given by

(3.11) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The configuration of the above 3-D stencils in (3.9), (3.10), and (3.11) combine the central FD [h.sup.-2][1, -2, 1] and both first-order 1-D mimetic stencils given by the second and third rows of (2.5). As a result, the local truncation error of these 3-D Laplacian approximations is formally O(h). All other stencils in the new mimetic scheme have seven points which exhibit second-order truncation errors at the grid interior. Since the new scheme is developed on an equidistributed Cartesian grid, the equations discretizing the boundary conditions are identical to those described in one-dimensional mimetic schemes and yields second-order truncation errors; see [7] and [14].

The above consistency analysis reveals that the mimetic discrete Laplacian approximations are of first order in the vicinity of the grid boundaries. However, this low nominal convergence is just the result of our simplified convergence analysis, but the global empirical convergence might be higher for several numerical tests and is mainly dominated by the interior second-order accuracy. Computational experiments discussed later give evidence that the new mimetic scheme actually yields optimal second-order convergence rates.

4. Eigenvalue estimation and stability considerations. An application of the discrete integration by parts formula (2.2) extended to our 3-D formulation shows that any eigenvalue associated to the system matrix in (3.5) is real and positive. Let us denote this matrix by [PHI], which for positive and constant parameters [kappa], [omega], a, and [gamma] can be written as

(4.1) [PHI] = -[D.sup.3D] [G.sup.3D] + [LAMBDA] + [GAMMA][B.sup.3D][G.sup.3D].

The matrix [LAMBDA] is diagonal with entries following the same lexicographical order as the discrete field [??], and it is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In the following and for the sake of clarity of notation, we omit the dimension superscript in the mimetic differentiation and quadrature operators in (4.1) and the 3-D generalization of (2.2).

PROPOSITION 4.1. The eigenvalues of the matrix [PHI] in (4.1) are real and positive.

Proof. In the case of a non-zero complex vector [??], where [??] = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the 3-D analog of (2.2) satisfies the identity

(4.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Here, we focus on second-order discretizations, where Q = hI. Thus, (4.2) is equivalent to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

where the right-hand side is always non-negative given that the diagonal weight matrix P is positive definite. Next, we suppose that the vector u corresponds to an eigenvector of the matrix [PHI] with the associated eigenvalue [lambda]. The Rayleigh quotient for this pair results in

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Thus, [lambda] is real and positive given that the matrix [LAMBDA] is also diagonal with entries larger than zero.

Proposition 4.1 proves the non-singularity of the system matrix [PHI] for the case of general Robin boundary conditions (3.2) coupled to the elliptic model problem (3.1). Actually, the above proof easily extends to the case of positive and spatially variable coefficients [kappa], [omega], a, and [gamma]. Moreover, in the particular case of Dirichlet boundary conditions, a simple application of Gershgorin's theorem [9] renders lower bounds for the eigenvalues of the matrix [PHI]. In the following result, we consider the same additional constraints of [kappa] = a = [gamma] = 1 imposed on the simplified versions of the mimetic Laplacian stencils (3.6)-(3.11).

PROPOSITION 4.2. Any eigenvalue [lambda] of the matrix [PHI] in (4.1) satisfies

[omega] [less than or equal to] [lambda] [less than or equal to] [omega] + O([h.sup.-2]).

Proof. The proof is a simple application of Gershgorin circles. We first consider an interior stencil of the system (3.5) that corresponds to the standard seven-point FD Laplacian approximation that is far enough away from the boundary contribution of the mimetic B in (2.9),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

where the magnitude of the off-diagonal entries add up to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

and the diagonal entry corresponds to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Thus, according to the Gershgorin circles,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

We next consider the unconventional stencil in equation (3.6). In this case,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

After an appropriate spatial scaling of the original elliptic problem, it is always the case that 0 < h < 1. In this interval, we have that ([h.sup.-2] - [h.sup.-1]/2) > 0, and the Gershgorin circles for the stencil (3.6) leads to the estimate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Notice that the strict inequalities hold for the above eigenvalue bounds because the Dirichlet boundary data actually add to the right-hand side vector in the system (3.5). The Gershgorin bounds for the remaining equations of this system are obtained in a similar way.

Proposition 4.2 reinforces the fact that the matrix [PHI] is invertible as the minimum eigenvalue remains bounded from below by [omega] > 0 independent of how h [right arrow] 0 under any grid refinement. In this sense, the mimetic scheme is stable and must converge to the exact solution of the given problem (3.1) and (3.2). On the other hand, if the Gershgorin upper bound is tight, then the maximum eigenvalue may grow as O([h.sup.-2]), and we also expect a similar increase of the condition number of the matrix [PHI]. This is the case for standard FD schemes with a symmetric discretization matrix for this type of elliptic problems. Now, a similar application of Gershgorin's theorem to our scheme matrix in (3.5) under Neumann or Robin boundary conditions does not lead to simple bounds for the eigenvalues. Thus, we empirically extend this eigenvalue analysis to Robin conditions by using MATLAB built-in functions and plot the smallest and the largest eigenvalue of the matrix [PHI] in Figure 4.1 as well as the [L.sub.2] condition number as the number of grid cells N increases. These results correspond to parameter values [kappa] = a = [gamma] = 1 and also [omega] = 1 with Robin-type boundary conditions. We scale the extreme eigenvalues [[lambda].sup.MIN] and [[lambda].sup.MAX] by the parameters [omega] and [N.sup.2], respectively, in order to make a correspondence to the above theoretical bounds for the Dirichlet problems. Notice the perfect scaling followed by both eigenvalues in these experiments. In addition, Figure 4.1 displays the [L.sup.2] condition number of [PHI] weighted by the ratio [lambda]MAX/[lambda]MIN even though this matrix is far from symmetric, and we have denoted this new term by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Note that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is directly proportional to [N.sup.2], so one may expect a linear increase of the number of iterations required to solve the system (3.5) by a standard iterative solver for sparse matrices.

We end this section with a brief comment on the properties of the system matrix [PHI] under a more general parametrization of the BVP (3.1)-(3.2). This matrix is non-symmetric for Robin boundary conditions [gamma] [not equal to] 0 because of the mimetic lateral approximation of the fluxes by the operator [B.sup.3D][G.sup.3D]. It may become indefinite according to the sign variation of the parameters [kappa] and w in the problem domain. Thus, the implementation of the mimetic scheme (3.5) requires specialized iterative solvers for non-symmetric and indefinite matrices. The authors in [16] present a 2-D discretization of the Poisson equation based on the mimetic divergence and gradient operators. This study carefully explores the application of the Generalized Minimal Residual (GMRES) and the BiConjugate Gradient (BiCGstab) methods to these problems, and recommends using the latter method with an incomplete LU preconditioning. Even though the FD discretizations in [16] omit the contribution of the operator [B.sup.3D][G.sup.3D] for Laplacian approximations, we here consider BiCGstab as one choice for solving the system.

5. Higher-order operators and locally refined grids. Two advantages offered by the tensor formulation of the 3-D gradient and divergence operators in equations (2.7) and (2.8) is the effortless construction of higher-order accurate counterparts as well as their extension to more general rectangular grids with a local refinement in any coordinate direction. In [5] and [4], Castillo and collaborators propose a three-parametric fourth-order family for each [G.sup.1D] and [D.sup.1D] operator along with consistently accurate quadratures [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Using a particular choice of [G.sup.1D] and [D.sup.1D] in combination with P and Q in the mimetic identity (2.3) yields the corresponding 1-D boundary operator. A direct substitution of these 1-D operators into the tensor formulae (2.7), (2.8), and (2.9) lead to 3-D fourth-order accurate operators and allows us to build a new higher-order formulation of the mimetic scheme (3.5). In the next section, we use the following fourth-order discretizations of [G.sup.1D] and [D.sup.1D] to solve some test problems numerically:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The second feature of the 3-D tensor formulation is an easy inclusion of locally refined rectangular grids. In [3], uniform [G.sup.1D], [D.sup.1D], P, and Q operators are adapted to 1-D refined meshes by using the standard mapping technique. The authors consider that a given arbitrary-spaced 1-D grid smoothly transforms to the uniform staggered mesh defined above in Section 2.1. The Jacobian of the transformation of the nodal grid is given by the diagonal matrix [J.sub.D], while the Jacobian of the complementary mesh comprising centers and boundary points is given by a different diagonal matrix [J.sub.G], both of them are readily defined as

[J.sub.D] = [D.sup.1D] [([x.sub.0], [x.sub.1],..., [x.sub.N]).sup.T], [J.sub.G] = [G.sup.1D] [([x.sub.0], [x.sub.1/2],..., [x.sub.N-1/2], [x.sub.N]).sup.T].

Thus, the mimetic operators on a 1-D non-uniform grid can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In [3], the authors also show that the 1-D operator B remains unchanged under any smooth coordinate mapping by employing a simple substitution of these non-uniform operators into (2.3). Therefore, the 3-D operators defined by the equations (2.7) and (2.8) can easily accommodate a variable grid stepping in any coordinate direction by replacing the uniform [G.sup.1D] and [D.sup.1D] operators by their non-uniform versions [G.sup.1D] and [D.sup.1D] as long as the 3-D grid remains rectangular. In the numerical experiments given below, we also exploit this virtue of our 3-D tensor operators to gain accuracy with almost the same computational cost as uniformly-spaced simulations.

Finally, we would like to note that a simplified boundary operator [B.sup.3D] can be defined by using B = [-[e.sub.1], 0,..., 0, [e.sub.N]] as the 1-D operator in equation (2.9) instead of the mimetic version given by (2.4). This simpler [B.sup.3D] can also couple Neumann or Robin boundary conditions to the scheme in (3.5) but leads to a more conventional FD elliptic scheme with a local truncation errors only depending on the nominal construction accuracy of the [G.sup.3D] and [D.sup.3D] differentiators.

6. Numerical tests. To assess the accuracy and experimental convergence of the proposed mimetic discretizations in (2.9), we use the following diffusion test problem with a diagonal variable tensor [kappa]I:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The exact solution is available and corresponds to

u = [x.sup.m] + [y.sup.m] + [z.sup.m].

The difficulty in this numerical test increases with the degree m of the polynomial solution as the boundary layers become harder to be resolved in the vicinity of the domain corners (with the exception of the origin point). First, we perform second-order simulations for m = 5, 6, 8, and 10 using square grids with an increasing number of cells N = 10, 20, 40, 50, and 60 along each coordinate direction. On the same grids, we also compute numerical solutions using the fourth-order mimetic operators detailed in the previous section. Finally, we use the common grid clustering T([epsilon]) = arctan([epsilon])/arctan(1) along each physical coordinate x, y, and z of the unit cube, where [epsilon] follows a uniform discretization of the computational interval [0, 1]. Thus, the grid spacing reduces away from the point (0, 0, 0) to allow for a better resolution of the steep gradients of the solution at the other corners. In this set of non-uniform tests, our choice for the accuracy was of second order. A similar arctan mesh adaptation was used in [3] to efficiently solve 1-D elliptic problems in the presence of boundary layers, and we here apply it to our 3-D tests. As the numerical results below show, this grid clustering allow reaching a higher accuracy than for the simulations on uniform grids with the same cell density. However, we want to note that our selection of this arctan adaptation has not been guided by any optimal or automatic strategy of grid refinement, and our goal is simply to exercise the easy adaptability of our mimetic tensor formulation to orthogonal locally refined grids.

All numerical experiments were carried out with MATLAB R2102 on an Intel i7-4600U with 8 GB RAM. To solve the linear systems given by (3.5), we use the built-in implementation of the Quasi-Minimal Residual Method qmr.m with SSOR preconditioning, whose successful performance in solving difficult elliptic equations after FD discretizations has been recently reported [12]. In all our experiments, qmr converges with a relative tolerance of [10.sup.-4] in less than a thousand iterations. In addition, we follow [16] and apply the MATLAB implementation of the BiConjugate Gradient method to these experiments, but it converges slower than qmr.m on our dense grids, and the computation of the ILU preconditioner also increases the computational demand.

Figure 6.1 depicts the [L.sub.[infinity]] relative errors of the solutions for the second-order approximation that use a constant grid step (panel a) and locally refined grids (panel b) in addition to the similar misfits given by the fourth-order simulations on uniform grids (panel c). In all cases, a perfect second-order convergence is also displayed as reference. One can see that the solution accuracy effectively decays as m grows in response to the problem difficulty. For a given number of grid cells, the errors are larger in the case of uniform second-order discretizations, while the fourth-order solutions give the most accurate approximations. A complementary accuracy information is given in Table 6.1, where we list the [L.sub.[infinity]] errors for each computed solution in the case of m = 10. Relative to the misfits of the uniform second-order simulations, the error present in the second-order non-uniform solution corresponds to 31%, while the uniform fourth-order misfits is only 13% for the available coarsest grid 10 x 10 x 10. Similar relative comparisons when using 60 x 60 x 60 grid cells, indicate that the former accuracy advantage stay in 32%, while the latter ratio drops significantly to 0.02%. The accuracy gain offered by any of these numerical devices, grid clustering or fourth-order discretizations, can be also assessed from the information in Table 6.1 and actually leads to a reduction of the computing time as shown by Figure 6.2. The elapsed times in this figure only count the construction and the qmr solution of the linear system (3.5) for each experiment where m = 10. If we set the [L.sub.[infinity]]-accuracy target to 2%, then we need to perform second-order simulations either on a 60 x 60 x 60 uniform grid or on a 40 x 40 x 40 locally refined grid, while the fourth-order scheme exceeds this tolerance on a 20 x 20 x 20 uniform grid. The computation times of these simulations were 40.4, 7.13, and 1.28 seconds, respectively, as depicted in Figure 6.2. This figure clearly illustrates the benefits of these aforementioned numerical enhancements on the grids with 30 x 30 x 30 cells and more. Note that higher precision is accompanied by a slight increase of computing times that remain in the same order of magnitude as the CPU time spent by the second-order uniform simulations. These advantages are better realized by the fourth-order scheme on rectangular grids.

Table 6.2 reports the least-square regression slopes to the error data, and the results on the uniform grids show the excellent agreement with respect to the nominal convergence rates for different degrees m. In addition, the second-order solutions on the locally refined grids present the typical convergence deficiency associated to variable spacing, but it is well compensated by its higher accuracy with respect to the results on uniform meshes.

7. Conclusions. In this work, we first construct second-order accurate mimetic gradient G, divergence D, and boundary B operators on 3-D rectangular meshes by using tensor products of 1-D operators defined along each grid coordinate. The mimetic Laplacian approximation is embedded into the matrix operator BG + DG, which here is used for the discretization of certain families of elliptic differential equations with Robin boundary conditions. Then, we carry out a full consistency analysis that leads to standard second-order Taylor truncation errors at the interior grid points but also reveals that non-standard mimetic stencils are formally first-order accurate near the boundaries. Linear systems arising from such discretizations lead to non-symmetric matrices with positive eigenvalues under a wide range of model parameters, and therefore the numerical stability and convergence are guaranteed. Next, the flexibility of the tensor formulation is demonstrated by an easy construction of 3-D fourth-order accurate operators starting from their 1-D counterparts. We also reformulate all proposed 3-D mimetic operators to include variable spacing along any coordinate direction but where the grids remain orthogonal. Finally, we assess the impact of both numerical features, fourth-order accuracy and local grid refinement, on the performance of mimetic discretizations by solving an elliptic problem with a parametrized boundary layer. This parameter is used to gradually increase the problem difficulty. On uniform grids, global experimental convergence rates are consistent with the nominal discretization as the grid spacing reduces, and the accuracy slightly degrades as the problem becomes harder. Thus, the accuracy of fourth-order uniform simulations rapidly exceeds the precision for second-order results with a slight increase of computing times. On the other hand, local grid refinements allow an efficient resolution of the exact steep gradients, and the resulting errors are approximately one third of those achieved on uniform equally dense grids.

Acknowledgement. The first and second authors are grateful to Banco Central de Venezuela (BCV) and Universidad Central de Venezuela (CDCH-UCV) for travel funding. We would also like to thank Professors Beatriz Otero and Eliecer Correa for their unreserved assistance throughout this work.

REFERENCES

[1] M. ABOUALI AND J. CASTILLO, Solving Poisson equation with Robin boundary condition on a curvilinear mesh using high order mimetic discretization methods, Math. Comput. Simulation, in press, online 2014.

[2] D. APPELO AND N. A. PETERSSON, A stable finite difference method for the elastic wave equation on complex geometries with free surfaces, Commun. Comput. Phys., 5 (2009), pp. 84-107.

[3] E. D. BATISTA AND J. E. CASTILLO, Mimetic schemes on non-uniform structured meshes, Electron. Trans. Numer. Anal., 34 (2008/09), pp. 152-162. http://etna.ricam.oeaw.ac.at/vol.34.2008-2009/pp152-162.dir/pp152-162.pdf

[4] J. E. CASTILLO AND R. D. GRONE, A matrix analysis approach to higher-order approximations for divergence and gradients satisfying a global conservation law, SIAM J. Matrix Anal. Appl., 25 (2003), pp. 128-142.

[5] J. E. CASTILLO, J. M. HYMAN, M. SHASHKOV, AND S. STEINBERG, Fourth- and sixth-order conservative finite difference approximations of the divergence and gradient, Appl. Numer. Math., 37 (2001), pp. 171-187.

[6] J. E. CASTILLO AND G. MIRANDA, Mimetic Discretization Methods, CRC Press, Boca Raton, 2013.

[7] J. E. CASTILLO AND M. YASUDA, Linear systems arising for second-order mimetic divergence and gradient discretizations, J. Math. Model. Algorithms, 4 (2005), pp. 67-82.

[8] L. BEIRAO DA VEIGA, K. LIPNIKOV, AND G. MANZINI, The Mimetic Finite Difference Method for Elliptic Problems, Springer, Cham 2014.

[9] J. W. DEMMEL, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997.

[10] K. DURU AND E. M. DUNHAM, Dynamic earthquake rupture simulations on nonplanar faults embedded in 3D geometrically complex, heterogeneous elastic solids, J. Comput. Phys., 305 (2016), pp. 185-207.

[11] G. ELY, S. DAY, AND J.-B. MINSTER, A support-operator method for viscoelastic wave modelling in 3-D heterogeneous media, Geophys. J. Internat., 172 (2008), pp. 331-344.

[12] A. GRAYVER, Three-Dimensional Controlled-Source Electromagnetic Inversion Using Modern Computational Concepts, PhD. Thesis, Department of Geosciences, Freie Universitat Berlin, Berlin, 2013.

[13] J. GUEVARA-JORDAN AND J. ARTEAGA-ARISPE, A second order mimetic approach for tracer flow in oil reservoirs, Proceedings Latin American and Caribbean Petroleum Engineering Conference X, SPE-107366-MS, Society of Petroleum Engineers, Richardson, US, 2007, pp. 626-633.

[14] J. M. GUEVARA-JORDAN, S. ROJAS, M. FREITES-VILLEGAS, AND J. E. CASTILLO, Convergence of a mimetic finite difference method for static diffusion equation, Adv. Difference Equ., (2007), Art. ID 12303 (12 pages).

[15] V. GYRYA AND K. LIPNIKOV, High-order mimetic finite difference method for diffusion problems on polygonal meshes, J. Comput. Phys., 227 (2008), pp. 8841-8854.

[16] F. F. HERNANDEZ, J. E. CASTILLO, AND G. A. LARRAZABAL, Large sparse linear systems arising from mimetic discretization, Comput. Math. Appl., 53 (2007), pp. 1-11.

[17] J. HYMAN, J. MOREL, M. SHASHKOV, AND S. STEINBERG, Mimetic finite difference methods for diffusion equations, Comput. Geosci., 6 (2002), pp. 333-352.

[18] J. M. HYMAN AND M. SHASHKOV, Adjoint operators for the natural discretizations of the divergence, gradient and curl on logically rectangular grids, Appl. Numer. Math., 25 (1997), pp. 413-442.

[19]--, Mimetic finite difference methods for Maxwell's equations and the equations of magnetic diffusion, Progr. Electromagn. Res., 32 (2001), pp. 89-121.

[20] D. JUSTO, High Order Mimetic Methods and Absorbing Boundary Conditions, VDM Verlag, Saarbrucken, 2009.

[21] J. E. KOZDON, E. M. DUNHAM, AND J. NORDSTROM, Simulation of dynamic earthquake ruptures in complex geometries using high-order finite difference methods, J. Sci. Comput., 55 (2013), pp. 92-124.

[22] H. KREISS AND G. SCHERER, Finite element and finite difference methods for hyperbolic partial differential equations, in Mathematical Aspects of Finite Elements in Partial Differential Equations, C. de Boor, ed., Academic Press, New York, 1974, pp. 195 - 212.

[23] P. OLSSON, Summation by parts, projections, and stability. I, Math. Comp., 64 (1995), pp. 1035-1065.

[24]--, Summation by parts, projections, and stability. II, Math. Comp., 64 (1995), pp. 1473-1493.

[25] J. RUNYAN, A Novel Higher Order Finite Difference Time Domain Method Based On The Castillo-Grone Mimetic Curl Operator With Applications Concerning The Time-Dependent Maxwell Equations, Master's Thesis, Computational Science Research Center, San Diego State University, San Diego 2011.

[26] E. SANCHEZ, C. PAOLINI, P. BLOMGREN, AND J. E. CASTILLO, Algorithms for higher-order mimetic operators, in Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2014, R. M. Kirby, M. Berzins, and J. S. Hesthaven, eds., Lecture Notes in Computational Science and Engineering 106, Springer, Cham, 2015, pp. 425-434.

[27] E. J. SANCHEZ, C. P. PAOLINI, AND J. E. CASTILLO, The mimetic methods toolkit: an object-oriented API for mimetic finite differences, J. Comput. Appl. Math., 270 (2014), pp. 308-322.

[28] F. SOLANO-FEO, J. M. GUEVARA-JORDAN, O. ROJAS, B. OTERO, AND R. RODRIGUEZ, A new mimetic scheme for the acoustic wave equation, J. Comput. Appl. Math., 295 (2016), pp. 2-12.

[29] M. SPIEGEL, Calculus of Finite Differences and Differences Equations, McGraw Hill, New York, 1994.

[30] B. STRAND, Summation by parts for finite difference approximations for d/dx, J. Comput. Phys., 110 (1994), pp. 47-67.

J. BLANCO ([dagger]), O. ROJAS ([dagger][parallel]), C. CHACON ([double dagger]), J. M. GUEVARA-JORDAN ([section]), AND J. CASTILLO([paragraph])

(*) Received May 29, 2016. Accepted October 17, 2016. Published online on December 9, 2016. Recommended by Andy Wathen.

([dagger]) Centro de Investigacion de Operaciones y Modelos Matematicos, Facultad de Ciencias, Universidad Central de Venezuela, Caracas 1040 (jaime.blanco@ciens.ucv.ve, rojasotilio@gmail.com).

([double dagger]) Area de Matematicas, Universidad Nacional Abierta, Caracas (ccsuescun@gmail.com).

([section]) Escuela de Matematicas, Facultad de Ciencias, Universidad Central de Venezuela, Caracas 1040 (jmguevarajordan@gmail.com).

([paragraph]) Computational Science Research Center, San Diego State University, San Diego, USA (jcastillo@mail.sdsu.edu).

([parallel]) Barcelona Supercomputing Center (BSC), Jordi Girona 29, 08034-Barcelona, Spain (otilio.rojas@bsc.es).

TABLE 2.2 Grid image of mimetic gradient and divergence approximations. Operator Corner-aided [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.n] * [Y.sub.cb] * [.sub.Zcb] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.cb] * [Y.sub.n] * [Z.sub.cb] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.cb] * [Y.sub.cb] * [Z.sub.n] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.c] * [Y.sub.c] * [Z.sub.c] Operator Conventional [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.n] x [Y.sub.c] x [Z.sub.c] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.c] x [Y.sub.n] x [Z.sub.c] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.c] x [Y.sub.c] x [Z.sub.n] [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [X.sub.c] x [Y.sub.c] x [Z.sub.c] TABLE 6.1 Relative [L.[infinity]] errors of numerical solutions in the case of m = 10. Grid size 2nd-order (uniform) 2nd-order (refined) 4th-order (uniform) 10 x 10 x 10 0.752379 0.232773 0.09797807 20 x 20 x 20 0.187049 0.058873 0.00390037 30 x 30 x 30 0.082925 0.026355 0.00044975 40 x 40 x 40 0.046690 0.014881 0.00007690 50 x 50 x 50 0.029914 0.009546 0.00001322 60 x 60 x 60 0.020789 0.006640 0.00000426 TABLE 6.2 [L.[infinity]] convergence rates of the 2nd- and 4th-order numerical solutions. m 2nd-order (uniform) 2nd-order (refined) 4th-order (uniform) 5 2.005 1.944 4.365 6 2.006 1.939 4.590 8 2.005 1.974 5.253 10 2.003 1.985 5.618

Printer friendly Cite/link Email Feedback | |

Author: | Blanco, J.; Rojas, O.; Chacon, C.; Guevara-Jordan, J.M.; Castillo, J. |
---|---|

Publication: | Electronic Transactions on Numerical Analysis |

Article Type: | Report |

Date: | Jan 1, 2016 |

Words: | 9247 |

Previous Article: | Dirichlet-Neumann and Neumann-Neumann waveform relaxation algorithms for parabolic problems. |

Next Article: | Weighted Hermite quadrature rules. |

Topics: |