Abstract
Elliptic grid generation equations based on the Laplacian operator have the well-known property of clustering the mesh near convex boundaries and declustering it near concave boundaries. In prior work, a new differential operator was derived and presented to address this issue. This new operator retains the strong smoothing properties of the Laplacian without the latter’s adverse curvature effects. However, the new operator exhibits slower convergence properties than the Laplacian, which can lead to increased turnaround times and in some cases preclude the achievement of convergence to machine accuracy. In the work presented here, a Newton linearization of the new operator is presented, with the objective of achieving more robust convergence properties. Sample solutions are presented by evaluating a number of solvers and preconditioners and assessing the convergence properties of the solution process. The efficiency of each solution method is demonstrated with applications to two-dimensional airfoil meshes.
1. Introduction
Elliptic grid generation is a powerful tool for optimizing the quality of a mesh in a complex domain. The solution of elliptic equations based on the Laplace operator can be used to control virtually all the characteristics of a given mesh, through the proper specification of control functions. Various forms of these equations have been proposed, based on the early works of Winslow [1], Barfield [2], Chu [3], Godunov and Prokopov [4], Amsden and Hirt [5], and Thompson et al. [6–10].
However, elliptic grid generation equations based on the Laplacian operator have the well-known property of clustering the mesh near convex boundaries and declustering it near concave boundaries. The result of this effect is that the mesh spacing near curved boundaries may not reflect the spacing prescribed by the user. In a prior work [11, 12], a new differential operator was derived to address this issue. This new operator retains the strong smoothing properties of the Laplacian operator without the latter’s adverse curvature effects. The capability of this new operator, herein referred to as the curvature operator, is illustrated in Figures 1(a) and 1(b) for a simple mesh over a cylinder with equally-spaced boundary points. It can be seen from Figure 1(b) that the mesh generated with curvature control yields the desired cell sizes near the concave corners at the top and bottom of the cylinder, whereas the mesh generated without curvature control increases the mesh spacing in these areas. The capabilities of this new operator have been demonstrated in numerous applications involving the generation of meshes in both two and three dimensions on complex geometry such as complete aircraft configurations [11].
However, it has been observed that the convergence of the numerical solution of the curvature operator is more problematic and slower than the solution of the basic Laplace operator. The main reason for this is that the curvature operator is inherently more sensitive to mesh curvature, and thus, any spurious or discontinuous curvature distribution arising during the solution process may cause a slowdown in convergence. The approach employed thus far to deal with this issue has been to simply under-relax the solution process. However, this approach leads to increased turnaround times, and in some cases precludes the achievement of full convergence to machine accuracy.
There is extensive literature on the numerical solution of similar systems of nonlinear partial differential equations. Davis and McCammon [13] showed that the incomplete Cholesky conjugate gradient (ICCG) method is superior to relaxation methods for solving a system of linear equations arising from the finite-difference form of the linearized Poisson–Boltzmann equation. Nicholls and Honig [14] developed an efficient algorithm for the numerical solution of the Poisson–Boltzmann equation by the finite-difference method utilizing successive over-relaxation through the rapid estimation of the optimum relaxation parameter. Cai et al. [15] implemented and optimized seven finite-difference solvers for the full nonlinear Poisson–Boltzmann equation, including four relaxation methods, one conjugate gradient method, and two inexact Newton methods. Wang and Luo [16] assessed five commonly used finite-difference solvers for the Poisson–Boltzmann equation and found that the modified incomplete Cholesky conjugate gradient and geometric multigrid are the most efficient in the diversified test set.
Knoll et al. [17] developed a nonlinear solution method for the nonequilibrium radiation diffusion problem by combining an outer Newton-based iteration with the inner conjugate gradient-like (Krylov) iteration.
Rider et al. [18] developed a robust and scalable solution of the equilibrium radiation diffusion problem using the GMRES Newton–Krylov method preconditioned by a multigrid method resulting from a Picard-type linearization of the governing equation.
Berndta et al. [19] developed two nonlinear solvers based on the Jacobian-Free Newton–Krylov (JFNK) methodology to solve the Laplace–Beltrami system of nonlinear, elliptic, and partial differential equations discretized using a finite-element method. They also developed two different preconditioners, both of which employ existing algebraic multigrid (AMG) methods. These solvers have been demonstrated to be significantly faster than a standard Newton–Krylov approach.
Gupta et al. [20] combined a compact high-order finite-difference approximation with a multigrid V-cycle algorithm to solve the two-dimensional Poisson elliptic differential equations and showed that the nine-point discretization formula, combined with a full-weighting projection operator, is much more accurate than the five-point discretization formula. Zhang [21] employed a fourth-order compact finite-difference scheme (FOS) with the multigrid algorithm to solve the three-dimensional Poisson equation and showed that it yields a fast and high-accuracy 3D Poisson solver.
Ilić et al. [22] showed that the fractional Poisson equation can be approximately solved using a finite-difference discretization of the Laplacian to produce an appropriate matrix representation of the operator and proposed an algorithm based on a Krylov subspace method that could be used efficiently to solve the resulting linear system.
Froese and Oberman [23] developed a fast and accurate finite-difference solver for the elliptic Monge–Ampère equation using a hybrid finite-difference discretization which selects between an accurate standard finite-difference discretization and a stable monotone discretization.
In the work presented herein, the grid generation equations with and without curvature control functions are discretized using finite-differences and a Newton linearization of the discretized equations is then implemented. To solve the resulting linear system, a number of solvers and preconditioners drawn from the Portable Extensible Toolkit for Scientific Computation (PETSc) library are implemented and evaluated.
In what follows, Section 2 provides an overview of the curvature operator and its derivation, followed by a detailed description of the methodology employed to linearize the equations. Section 3 presents the numerical implementation process and setup of the matrix equation, and Section 4 presents applications of the solution process to C-meshes around a two-dimensional airfoil.
2. Mathematical Derivation
2.1. Overview of the Curvature Operator
In this section, a brief overview of the curvature operator and its derivation will be presented followed by a detailed description of the methodology employed to linearize the equations.
2.1.1. Mesh Generation Equations
Denoting as the curvilinear coordinates, it was shown in [11, 24] that for a mesh with uniform spacing, the curvature operator can be expressed as follows:where represents the curvature control functions and is given by the following:whereand where are cyclic, with no sum on and .
In the above expressions, the are the contravariant metric tensor components, and the are the Christoffel symbols of the second kind. The physical interpretation of equation (3) can be made by noting that, for a surface, represents the sum of the local principal curvatures of the surface.
As shown in [11], in the general case of a nonuniform mesh, the curvature operator must be expressed in terms of the functions defining the grid clustering, denoted here as . The functions are defined by a mapping of the spacing distributions of the physical mesh onto a computational space scaled to a unit square, as shown in Figure 2. The points along the boundaries of the space represent the normalized arc length distributions along the boundaries of the physical space, while the points in the interior of the space are interpolated algebraically from the boundaries [25].
Thus, for a nonuniform mesh, equation (1) is applied to the functions, as follows:whereandand where the bar over the terms in equations (4) to (6) denotes that they are functions of , not functions of the curvilinear coordinates [11].
The abovementioned relations can be recast in terms of the curvilinear coordinates by expressing the Laplacian of in equation (4) as a function of the curvilinear coordinates, which is given by the following expression:
It is to be noted that in equation (7) and similar expressions in the remainder of this work, the summation convention is used whenever repeated upper and lower dummy indices are present, unless otherwise noted. Isolating in equation (7) then leads to the following:where
In equation (9), the second term on the right-hand side represents the functions controlling the mesh spacing, as derived by [25]. Equation (10) gives the general form of the curvature control functions, , for a mesh with nonuniform spacing. For a mesh with uniform spacing, equation (9) reduces to equation (1).
The mesh generation equations are then built by using the identity (where ) and expressing it in terms of the curvilinear coordinates:
The above system of equations is solved to yield the mesh coordinates, with the terms given by equations (9) and (10).
2.1.2. Evaluation of Curvature Control Functions
The numerical generation of the mesh using equation (11) requires the evaluation of the terms using equations (9) and (10) and the evaluation of the terms using equations (5) and (6).
The solution process begins by defining the functions, following the methodology described in detail in [25]. These functions are initially calculated based on the mesh clustering on the boundaries and are interpolated into the interior of the domain. During the iterative solution process, the functions can be updated using the Neumann boundary conditions to enforce orthogonality conditions where desired, thereby modifying the functions both on the boundaries and in the domain at every iteration.
The functions are expressed in terms of the curvilinear coordinates , where . Once the functions have been defined, their first and second derivatives with respect to can be computed directly at every point in the domain via finite differences of the coordinate transformation. To evaluate the second term on the right-hand side of equation (9), it is also required to compute the derivatives of the curvilinear coordinates with respect to the functions, and these can be evaluated using the following well-known relations [26]:where and are cyclic, and where is the determinant of the Jacobian matrix of the -to- coordinate transformation:
To evaluate the terms in equation (10), the control functions must be evaluated with respect to the coordinates. To do this, we must compute the metric tensor components and Christoffel symbols of the Cartesian -to- coordinate transformation. The metric tensors of the transformation can be evaluated using the following well-known relations [27]:where and are cyclic. To evaluate equation (14), the derivatives of the Cartesian coordinates with respect to the coordinates are required, and these can be computed using the chain rule involving the coordinates:where the are evaluated using equation (12).
To evaluate the Christoffel symbols with respect to the coordinates in (6), denoted as , the following expression can be used [27]:
In turn, to evaluate the second derivative of the physical coordinates with respect to the coordinates, it is useful to express the second term on the right-hand side of equation (9) in the following compact form:
Substituting (19) back into equation (9) then leads to the following:
From equation (19) and the transformation laws for the Christoffel symbols [27], it can be shown that:
Expressions (12) to (21) can be used to evaluate all terms required to construct the curvature functions appearing in equations (9) and (10), and the grid generation equation (11) can then be solved numerically to generate the mesh.
2.2. Two-Dimensional Form of the Equations
In this work, the applications of the grid generation equations are limited to two-dimensional domains. To this end, if we let and in equation (11), we can write
In two dimensions, the transformation relations between the covariant and contravariant metric tensors reduce to the following:
Substituting relations (23) into equation (22), we then obtain the following:
If we now combine equation (23) with equation (20) and expand, we obtain
If we define the functions and as follows, and insert equations (25) and (26) into equation (24), we obtain
In the above equations, the components of the curvature control functions are represented by the and terms. If and are set to zero, the standard grid generation equations without curvature control functions are recovered.
2.2.1. Expansion of Curvature Control Functions in Two Dimensions
The general expression for the curvature control functions are given by equations (5) and (10). In two-dimensions, , , and the expansion of the Christoffel symbols given by equation (18), in combination with the relations given in (23), lead to the following expressions for the curvature control functions:where .
It can also be shown thatwhere is the determinant of the metric tensor of the transformation, is the determinant of the metric tensor of the transformation, is the determinant of the metric tensor of the transformation, and is the determinant of the Jacobian matrix of the transformation (13).
Inserting relations (29) into (28a)-(28b) then yields:where the first and second derivatives of r with respect to the coordinates can be computed using equations (17) and (21), respectively. The latter equation can be simplified further by defining the function as follows:
Equation (21) can then be rewritten in a more compact form:
The curvature functions in equation (27) are then evaluated using (10):
In this work, numerical solutions of the grid generation equation (27) will be evaluated both with and without the curvature control functions to quantify the impact of the latter on the convergence of the numerical solution.
2.3. Discretization Scheme
The discretization of equation (27) is accomplished using second-order finite-difference expressions for all first and second derivatives given by the following expressions:
The updated values of the coordinates at iteration , denoted as , are expressed in terms of the old values at iteration plus the increments from iteration to , denoted as follows:
Substituting finite-difference expressions (34) and (35) into (27), and linearizing the equations by neglecting all second-order terms, , leads to a linear system of equations that must be solved at every iteration. Although (27) is a vector equation, because the grid coordinates are coupled through the product of the metric tensor components and the Christoffel symbols, the coordinates must be solved simultaneously through the solution of one matrix equation encompassing the entire system of equations.
The methodology employed to linearize the discretized form of the grid generation equations is described in detail in the next section.
2.4. Linearization of Grid Generation Equations
In this section, the linearization of the grid generation equations is developed in two steps, first without the curvature control functions and then with the latter included. The linearization process involves the discretization of the equations using second-order finite differences followed by a Newton linearization of the discretized equations. As will be shown in the following sections, the addition of the curvature functions significantly increases the complexity of the derivations.
2.4.1. Linearization of Equations without Curvature Functions
The grid generation equations without the curvature functions are obtained by setting and in (27). Employing the notation described in Section 2.3, the updated value of equation (27) at iteration can be written as follows:
Expanding the above, we obtain the following:
Multiplying all terms in (37) and neglecting all higher-order terms, the equation reduces to
To simplify the subsequent derivations, the left-hand side of equation (38) is expressed as follows:where
In equation (40c), all terms must be expressed as functions of the grid coordinates, . To this end, expressions for and are derived from equations (26a) and (26b) as follows:
In the above expressions, the terms represent the grid stretching functions, which are only functions of the coordinates. Since the latter are held fixed during the solution process, they do not require the calculation of increments.
Substituting equations (41a) and (41b) into equation (40c) and collecting terms, we obtain the following:
If we now substitute equation (31) in equation (42), we can write,
We also need to express the change in the metric tensor components in terms of the grid coordinates. The following expressions can be used for this purpose:
Differentiating equation (44) then leads to the following:
Substituting equations (45a)-(45c) into equation (43) and collecting terms, we get an expression for in terms of the increments in the derivatives of the grid coordinates:
Inserting equations (46), (40a) and (40b) into (39), we then obtain the following:
Equation (47) is a vector equation in which , , and . To separate this equation into its and components, we note that from equation (31), we have the following:
Therefore, the and components of equation (47) can be written as
Equations (49) and (50) are the fully linearized versions of the grid generation equations without curvature control functions. In these equations, all the derivatives are computed using finite differences involving and , using the expressions given in (34). However, only the updates in the derivatives in these equations, i.e., the terms, involve and , while all other terms involve old values of and from the previous iteration. Thus, the discretized form of these equations leads to a system of linear equations that can be solved for and at every grid point. The numerical procedure to assemble the linear system of equations based on (49) and (50) is described in detail in Section 3.
2.4.2. Linearization of Curvature Control Functions
The curvature control functions in equation (27) are given by equation (33). The updated values of these terms at iteration , denoted as , are expressed in terms of the old values at iteration plus the increments from iteration to :
In equation (33), the stretching function derivative terms are held constant during the solution process, and thus is given by the following:where the terms are given by equations (30a) and (30b). In the latter equations, is the Jacobian of the coordinate transformation, and is held constant during the solution process. However, all the other terms, involving derivatives of with respect to the coordinates, or metric tensor components of the coordinate transformation, will vary at every iteration. Therefore, using the chain rule for differentiation, can be expressed asand is given by the following:
To evaluate equations (53) and (54), every term involving an update in the grid coordinates, i.e., every term preceded by a “,” must be expanded further. We begin by evaluating . From equation (17), can be written as follows:
Assuming that the stretching function derivative terms are held constant, can thus be written as
Therefore, expanding the vector dot product and collecting terms, can be expressed as follows:
Similarly, can be written as
To evaluate the terms, we assume again that the stretching function derivatives are held constant. Therefore, using equation (32), can be expressed as follows:
Similarly, evaluating using equation (31) and substituting it into equation (59), we obtain the following:
Equating the vector components on both sides of equation (60), we can write
Similarly, using equation (55), the components of can be expressed as follows:
Therefore, taking the dot product of and and substituting the expressions for and given by equations (61) and (62), can be expressed as
From the expression for the metric tensor in equation (14), can be expressed as follows:
Therefore, the change in the covariant metric tensor components for the coordinates can be expressed as follows:
If we then expand the dot products in equation (64) and rearrange, can be expressed as follows:
From the chain rule applied to the differentiation of quotients, we can write the following:
Substituting equation (66) into equations (67a)-(67b), we then have the following
Now, combining equations (58), (63), (67a) and (68a) into equation (53) and rearranging, we get the following expression for :
After rearranging and collecting terms, we can finally express as follows:whereand
Similarly, if we combine equations (58), (63), (67b) and (68b) into equation (54) and rearrange, we get the following expression for :
Rearranging and collecting terms in (73), we can then write as follows:whereand
Using expressions for and from equations (70) and (74), respectively, in equation (52) can now be expressed as
We now have all the terms necessary to evaluate in equation (51).
In Section 2.4.1, the derivation of the linearized grid generation equations ((49) and (50)) excluded the curvature terms . Adding these terms back into equation (27) and applying the linearization process described in Section 2.4.1, the following additional terms must be added to the left-hand side of equation (38):
Separating equation (78) into its and components then yields the following expressions:
and
If we now add expression (79a) to the left-hand side of equation (49) and expression (79b) to left-hand side of equation (50), and then substitute the expressions for and given by equation (77), we obtain the following complete linearized grid generation equations with curvature control functions:and
Equations (80) and (81) are the fully linearized equations that must be solved numerically at each grid point to compute the and coordinates. To facilitate the discretization of these equations, we now recast them in terms of the coefficients multiplying each of the increments in the derivative terms, i.e., the coefficients of and , leading to equations of the following form:whereandand where
In equations (84a)–(84j) to (86a)–(86b), the curvature control functions can be deactivated by setting all the and terms to zero. The numerical implementation of equations (82) and (83) is described in detail in the next section.
3. Numerical Implementation
3.1. Setup of the Matrix Equation
Writing the finite-difference expressions in (34) in terms of the change in and coordinates from iteration to , we have
In the mesh generation process, the curvilinear coordinates, and , are the independent variables, and for convenience, they are set equal to the mesh indices and , respectively. This implies that and vary from 1 to and , respectively, and thus everywhere in the mesh. Therefore, substituting the finite-difference expressions from equations (87a)–(87j) into equations (82) and (83), yields the following:and
For the ease of presentation, we introduce the parameters