Algebraic Construction of a Strongly Consistent, Permutationally Symmetric and Conservative Difference Scheme for 3D Steady Stokes Flow †

: By using symbolic algebraic computation, we construct a strongly-consistent second-order ﬁnite difference scheme for steady three-dimensional Stokes ﬂow and a Cartesian solution grid. The scheme has the second order of accuracy and incorporates the pressure Poisson equation. This equation is the integrability condition for the discrete momentum and continuity equations. Our algebraic approach to the construction of difference schemes suggested by the second and the third authors combines the ﬁnite volume method, numerical integration, and difference elimination. We make use of the techniques of the differential and difference Janet/Gröbner bases for performing related computations. To prove the strong consistency of the generated scheme, we use these bases to correlate the differential ideal generated by the polynomials in the Stokes equations with the difference ideal generated by the polynomials in the constructed difference scheme. As this takes place, our difference scheme is conservative and inherits permutation symmetry of the differential Stokes ﬂow. For the obtained scheme, we compute the modiﬁed differential system and use it to analyze the scheme’s accuracy.


Introduction
In this paper, we extend to the three-dimensional case the results of paper [1], where we generated and studied the strong consistent finite difference scheme for two-dimensional (2D) steady Stokes flow of an incompressible fluid. For the 3D steady Stokes flow, the governing system of partial differential equations (PDE) reads: Here, x, y, z are the independent variables; the velocities u, v, and w, the pressure p, and the external forces f (1) , f (2) , and f (3) are the dependent variables; the constant Re is the Reynolds number; and ∆ := ∂ xx + ∂ yy + ∂ zz is the Laplace operator. Equation (1) approximate the incompressible Navier-Stokes equations: when Re 1. Under the last condition, the nonlinear inertia terms in Equation (2) can be neglected (cf. [2], Secttion 22.11). We refer to the book [3] and to the references therein for the fundamental mathematical theory of the Stokes flow.
Equations (1) and (2) possess the permutational symmetry: In our finite difference approximation (FDA) to Equation (1), we want to preserve this symmetry, and for this purpose, we choose a Cartesian solution grid. Furthermore, as shown by the research results presented in [4], mimetic discretizations, i.e., such discrete approximations to PDE that mimic their basic algebraic properties, are more likely to produce highly accurate and stable numerical results. In addition to the symmetry (3), among the algebraic properties of the PDE (1) to be preserved at the discrete level, we make sure of the conservativity and strong consistency or s-consistency of FDA. Conservativity means inheritance at the discrete level of the underlying conservation laws of PDE (1). The conventional notion of consistency (cf. [5], Chapter 7) provides reduction of the FDA to the original PDE when the grid spacings tend to zero. In other words, a consistent discretization approximates PDE. s-consistency is the novel concept introduced in [6,7]. For Stokes flow, it means (see Definition 2) not only approximation of the equations in PDE by difference equations, but also approximation of elements in the differential ideal generated by the polynomials in (1) by elements in the difference ideal (cf. [8], Chapter 2) generated by the polynomials in the FDA to (1). In particular, if the s-consistency holds, then among the consequences of FDA, i.e., among elements in the ideal it generates, there are approximants of all integrability conditions (see [9], Chapter 2) for the original PDE.
Thus, our aim is to construct, given a Cartesian grid, a conservative FDA (difference scheme) to the governing Stokes system (1). The scheme contains a discrete version of the pressure Poisson equation (integrability condition) and is s-consistent with (1). For this purpose, we use the approach proposed in [10]. The approach is comprised of the finite volume method, numerical integration, and difference elimination. For the generated scheme, we apply the algorithmic criterion for verification of its s-consistency. This criterion was designed in [6] and relates Janet/Gröbner bases of the differential and difference ideals generated by PDE and FDA. The computational experiments performed in papers [11,12] with the Navier-Stokes equations revealed a substantially better numerical behavior of s-consistent schemes than that of s-inconsistent ones.
Equation (1) are linear. It makes construction and analysis of their solutions much easier than in the case of the Navier-Stokes equations. Besides, one can fully algorithmically generate difference approximations to Equation (1) and to analyze their s-consistency. To perform computation of Janet/Gröbner bases, we use two Maple packages implementing the involutive algorithm (cf. [13]): the package JANET [14] for linear differential systems and the package LDA [15] (linear difference algebra) for linear difference systems.
Having the difference scheme constructed, we compute a modified differential system of this FDA, which we refer to as modified Stokes flow, and use it to analyze the order of approximation of the scheme and its consistency. The method of modified equation suggested in [16] is widely used (see [17], Chapter 8 and [18], Section 5.5) in studying difference schemes. The method provides a natural and unified platform to study such basic properties of the scheme as consistency and stability, order of approximation, convergence, dissipativity, dispersion, and invariance. However, to our knowledge, the methods for the computation of modified equations have not been extended yet to non-evolutionary PDE systems, and in [1], by applying the technique of differential Janet/Gröbner bases to the 2D steady Stokes flow for the first time, a modified equation was constructed for non-evolutionary system of PDE.
The present paper is organized as follows. In Section 2, we generate a difference scheme for Equation (1) by applying the approach of paper [10]. In Section 3, we show that our scheme is s-consistent and demonstrate the s-inconsistency of another scheme obtained by a tempting compactification of the constructed scheme. The computation of a 3D modified Stokes system is performed in Section 4. Here, by the example of the s-inconsistent scheme of Section 3, we show how the modified Stokes system detects the s-inconsistency. Finally, some concluding remarks are given in Section 5.

Algorithmic Generation of the Difference Scheme for Stokes Flow
For this purpose, we choose the Cartesian solution grid with the grid spacing h and apply the approach of paper [10] to generate a difference scheme for Equation (1).
Step 1. Completion of to involution (we refer to [9] and to the references therein for the theory of involution). We fix the orderly lexicographic ranking of partial derivatives related to Equation (1) with: For this ranking, the Maple package JANET [14] outputs the Janet involutive form of Equation (1), which is the minimal reduced differential Gröbner basis form: Here, we underlined the leaders, i.e., the highest ranking partial derivatives occurring in Equation (5). The equation F 5 is the pressure Poisson equation. This equation is the integrability condition for system (1). It can be expressed in terms of the left-hand sides of the equations in (1) as: Remark 1. The differential polynomial F (2) in Equation (5) is that in Equation (1) reduced modulo the continuity equation F (1) .
Step 2. Conversion into the integral form. We choose the cubic control volume Ω of Figure 1 as a stencil for difference approximation of the partial derivatives and rewrite equations F (1) , F (2) , F (3) , and F (4) into the equivalent integral conservation law form.
where ∂Ω is the boundary of Ω, and we apply the Gauss-Ostrogradsky formula (8).
To preserve at the discrete level the symmetry of System (1) under the permutation Equation (3), we use in Equation (7) the original form of F (2) given in Equation (1) (see Remark 1).
Step 3. Addition of integral relations for derivatives. We add to System (7) the following exact integral relations: Here, x j , y k , z l are the points of boundary ∂Ω shown in Figure 1.
Step 4. Numerical evaluation of integrals. Now, to evaluate the integrals in Equation (7), we apply the midpoint rule for the integrals over ∂Ω, the trapezoidal rule for the integrals in (9), and approximate the triple integrals over Ω as: As a result, we obtain the system (10) of difference equations for the grid functions: approximating functions: and for the grid functions approximating their partial derivatives: j, k, l ∈ Z w x j, k, l ≈ w x (jh, kh, lh) , w y j, k, l ≈ w y (jh, kh, lh) , w z j, k, l ≈ w z (jh, kh, lh) .
Step 5. Difference elimination of partial derivatives. To eliminate the grid functions u x , u y , u z , v x , v y , v z , w x , w y , w z , we compute a difference Janet/Gröbner basis form of the set of linear difference polynomials occurring in the left-hand sides of Equation (10). The computation can be done with the Maple package LDA [6] if one selects the orderly lexicographic ranking, which is a difference analogue of the differential ranking (4) used on Step 1: with the right-shift operators σ 1 , σ 2 , and σ 3 acting as translations, for instance, In the output of the LDA, there are four difference polynomials that do not contain the grid functions u x , u y , u z , v x , v y , v z , w x , w y , w z . These polynomials comprise a difference scheme. Being interreduced, the obtained scheme does not possess symmetry under the transformation (3). For this reason, we prefer the following redundant, but conservative and symmetric form of the scheme: where ∆ 1 and ∆ 2 are finite difference discretizations of the Laplace operator acting on a grid function g j, k, l as: ∆ 2 g j, k, l := g j+4, k+2, l+2 + g j+2, k+4, l+2 + g j+2, k+2, l+4 − 6g j+2, k+2, l+2 4h 2 Remark 2. It should be noted that EquationF (5) in the system (13) can also be obtained from the integral form of F 5 in Equation (5) with the stencil of Figure 1 if one uses the midpoint rule for the surface integration of the p x , p y , and p z and for evaluation of the additional integral relations for these derivatives: x j+2 whereas the volume integrals of f (1) , f (2) , and f (3) are evaluated by the trapezoidal rule.
The difference scheme (13) approximates the involutive system (5) and demonstrates the full correspondence between the differential and difference Janet/Gröbner bases forms of these systems. Such correspondence is a consequence of our choice of the differential (4) and difference (11) rankings.

Consistency Analysis
] be the ring of differential polynomials over the differential field of rational functions in parameter Re. We consider the dependent variables of the Stokes flow (1) as differential indeterminates and their grid approximants as difference indeterminates (cf. [8], Chapter 2). Respectively, we denote byR = Q(Re, h){u, v, p, f (1) , f (2) , f (3) } the difference polynomial ring whose elements are polynomials in the grid functions with the right-shift operators σ 1 , σ 2 , and σ 3 (see (12)).
Every element in I vanishes on solutions of the Stokes flow (1), and every element inĨ vanishes on solutions of (13). We refer to elements in I (respectively, inĨ) as the consequences of Equation (1) (respectively, of Equation (13)).
It is clear that to approximate Equation (5), the scheme (13) must be pairwise consistent with the involutive differential form (5). We call this sort of consistency weak consistency or w-consistency (cf. [6]). The difference polynomial set {F (1) ,F (2) ,F (3) ,F (4) ,F (5) } in Equation (13) is w-consistent with the differential system (5) since: Remark 3. It is a universally-adopted notion of consistency for a finite difference discretization of PDE systems (cf. [5], Chapter 7) and means the reduction of Equation (13) to Equation (5) when the mesh step h goes to zero.

Definition 1. [6]
We shall say that a difference equationf = 0,f ∈R implies the differential equation f = 0, f ∈ R, and writef f if the Taylor expansion about a grid point, after clearing denominators containing h, yields:f The following definition establishes the consistency interrelation between the differential and difference ideals generated by Equations (1) and (13), respectively. If such a consistency holds, then it provides a certain inheritance by the difference scheme of the algebraic properties of Stokes flow.
Proof. The set of difference polynomials in Equation (13), by its construction, is a Janet/Gröbner basis of the elimination idealĨ 0 ∩R whereĨ 0 is the difference ideal generated by the polynomials in Equation (10) (cf. [19], Thm. 2.3.4). The same set is also a Janet/Gröbner basis for the idealĨ and ranking (11). It is readily verified with the LDA package. Hence, the w-consistency (16) implies s-consistency.

Remark 5.
Generally, given a linear PDE system and its FDA, the s-consistency conditions (20) can be algorithmically verified by using relevant built-in routines of the Maple packages JANET [14] and LDA [15].
It is clear that w-consistency (16) follows from s-consistency of Equation (13) with Equation (5). However, the converse is not true, i.e., w-consistency does not imply s-consistency, in general. Thus, for numerical simulation of the Stokes flow, it is tempting to replaceF (5) in Equation (13) with a more compact discretization: This substitution preserves w-consistency since: However, as the following Proposition shows, the obtained scheme is not s-consistent.
Since the differential polynomials F (6) and F (7) are irreducible modulo the differential ideal generated by the differential polynomials in Equation (1), it follows that F (5) and F (6) are not consequences of the Stokes equations. Therefore, there are solutions to the Stokes equations that do not satisfy Equation (23).

Remark 6.
As a result of s-inconsistency, Equation (23) impose constraints on the external forces, which are not consequences of the governing differential equations (1).

Modified Stokes Flow
To construct a modified equation (cf. [18], Sect. 5.5) for the discrete Stokes flow (13), a numerical solution of the governing differential system (1), given external forces f (1) , f (2) , and f (3) , should be considered as a set of continuous differentiable functions {u, v, w, p} whose values at the grid points satisfy the difference system (13). Since this system describes the differential system (5) only approximately, one cannot expect that a continuous solution interpolating the grid values exactly satisfies Equation (5). In fact, it satisfies another set of differential equations, which we call the modified steady Stokes flow or modified flow for short.
The method of modified differential equation treats the difference equations comprising the scheme as infinite order differential equations obtained by replacing the various shift operators in the difference equations by the Taylor series about a grid point. For equations of an evolutionary type, the next step is to eliminate all derivatives with respect to the evolutionary variable of order greater than one by their expressions via partial derivatives with respect to the other independent variables. This step allows obtaining a kind of canonical form of the modified equation. Then, truncation of the order of the differential representations in the grid steps gives a family of modified equations ("differential approximations") of the difference scheme.
As we show here, the fact that the polynomial parts of both differential and difference systems are Gröbner bases of the ideals, they generate and satisfy the condition (16), which implies the s-consistency condition (20) and allows developing a constructive procedure for the computation of the modified flow. Due to the fact that the finite differences in the scheme (13) approximate the partial derivatives occurring in Equation (5) with accuracy O(h 2 ), it is natural to expect the scheme to have the second order of accuracy. For this reason, we restrict ourselves to the computation of the second order modified flow.
The Taylor expansions of the difference polynomials in Equation (13) at the grid point: (j + 1, k + 1, l + 1) forF (1) ,F (2) ,F (3) ,F (4) and at the point: (j + 2, k + 2, l + 2) forF (5) yield: where we explicitly write the terms of order h 2 . The calculation of the right-hand sides in Equation (24), as well as the computation of the expressions given below were performed with the help of the freely-available Python library SYMPY (http://www.sympy.org/) for symbolic mathematics.

Remark 7.
The Taylor expansions of the elements in difference scheme (13) and also those in the scheme 1 } over the chosen grid points contain only the even powers of h. The reason is that all the finite differences occurring in the equations of both schemes are the central difference approximations of the partial derivatives occurring in (5).
Furthermore, let us reduce the terms of order h 2 in the right-hand sides of (24) modulo the differential Janet/Gröbner basis (13). Such reduction gives us a canonical form of the second order modified flow. This is because the normal form of a polynomial computed modulo a Gröbner basis basis is uniquely defined (cf. [19], Section 2.1). The normal form can be computed by means of the routine InvReduce in the Maple package JANET.
Thereby, the Taylor expansion of the difference polynomials produces the second order modified Stokes flow in the form of System (25). As we know, solutions of Stokes flow (1) satisfy the integrability condition (6), which we rewrite as: Substitution of the expansions (25) into the equality (26) yields that its left-hand side is equal to zero for the terms explicitly written in (25). As we show below, this is an implication of the s-consistency of the scheme.

Proposition 2.
A w-consistent difference scheme for involutive Stokes flow (5) is s-consistent only if its Taylor expansion based on the central-difference formulas for derivatives and reduced modulo system (5), after its substitution into the left-hand side of Equation (26), vanishes for every order in h 2 .
The implication (30) is a consequence of the fact that the normal form of the differential polynomial (30) modulo Equation (5), when it is nonzero, is not an element of the differential ideal generated by the polynomials in (5), which contradicts the s-consistency ofG. Because of the same argument, the equality (31) holds for any k.

Corollary 2.
Given ranking (11), a w-consistent difference scheme for involutive Stokes flow (5) is s-consistent if and only if its set of polynomials is a difference Janet/Gröbner basis.
Proof. "⇐" From our choice of the ranking and the corresponding structure (5) of the differential Janet/Gröbner basis whose leaders are underlined, it follows that a w-consistent difference scheme composed of five difference polynomials: has the only difference S-polynomial of the form (28) approximating the left-hand side of the differential integrability condition (26). Together with the Taylor expansion (27), the relations (29)-(31) imply the reduction of S-polynomial (28) to zero modulo the set (32). Thus, the scheme is a Janet/Gröbner basis. "⇒" If a w-consistent set (32) is a Janet/Gröbner basis, then by Theorem 1, it is s-consistent.
It is instructive to illustrate Proposition 2 and Corollary 2 by the s-inconsistent difference scheme: 1 ,F 1 ,F 1 ,F 1 } analyzed in Section 3. It has the first four difference equations equal to those in the system (13), 1 given by Equation (21). Because of the distinction ofF Expression (34) contains terms of second order in h. Up to the factor h 2 , the sum of these terms coincides with the differential polynomial F (7) in Equation (23). Thus, the presence of the second order terms in (34) is intimately related to the s-inconsistency of (33) with the Stokes flow (1). Apparently, the PDE system (33) cannot be considered as a modified Stokes flow.

Conclusions
For the three-dimensional incompressible steady Stokes flow (1) and a Cartesian solution grid, we applied a computer algebra-based approach suggested in [10] to derive the s-consistent, conservative and permutationally-invariant difference scheme (13) for which we constructed the modified Stokes flow (25). The structure of the latter shows that the derived scheme has order O(h 2 ).
Our constructive approach to derivation of the scheme combines the methods of differential and difference Gröbner bases. The first method applied to the governing system of Stokes equations (1) to complete this system to the involution form (5) that incorporates the pressure Poisson equation F (4) . The second method allows deriving the difference scheme on the chosen grid by means of difference elimination and verifying its s-consistency by applying the criterion of Theorem 1. This criterion is fully algorithmic for linear PDE systems [6].
In addition, we exploited both methods to construct the modified Stokes flow (25). The structure of the constructed modified flow is determined by the choice of the differential and difference rankings. We performed computational experiments with several rankings and finally preferred the ranking (4) for the differential case and (11) for the difference case as the best suited. All the related computations were done by means of the Maple packages JANET [14] and LDA [15] and the freely-available Python library SYMPY (http://www.sympy.org/) for symbolic mathematics.
Since our difference scheme (13) is obtained from the four governing equations (1) by constructing the difference Janet/Gröbner basis (see also Remark 2), we checked with the help of Gröbner bases whether the difference ideal generated byF := {F (2) ,F (3) ,F (4) } contains a difference approximation to the continuity equation. In the case of the existence of such an approximation, it might be used for the numerical study of Stokes flow in the velocity-pressure formulation. However, the straightforward computation with LDA shows that the discrete version of F (1) is not a consequence ofF. Therefore, to use the velocity-pressure formulation, one has to add information on the continuity equation toF via the corresponding boundary condition (cf. [20]).