# Newton Linearization of the Curvature Operator in Structured Grid Generation with Sample Solutions

Authors: Pat Piperni; Maahi M. Talukder
Journal: International Journal of Mathematics and Mathematical Sciences
 (2022)

---

## 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.


---

## Body

## 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].Figure 1
Grids generated without and with curvature control, (a) with no curvature control and (b) with curvature control.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, Section2 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ξi 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:(1)∇2ξi=Ci;i=1,2,3,where Ci represents the curvature control functions and is given by the following:(2)Ci=−K1i+K2igii,where(3)K1i+K2i=1ggii3/2Γjjigkk−2Γjkigjk+Γkkigjj,and where i,j,k are cyclic, with no sum on j and k.In the above expressions, thegij are the contravariant metric tensor components, and the Γjki are the Christoffel symbols of the second kind. The physical interpretation of equation (3) can be made by noting that, for a ξi=constant surface, K1i+K2i 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 sii=1,2,3. The si 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 si space represent the normalized arc length distributions along the boundaries of the physical space, while the points in the interior of the si space are interpolated algebraically from the boundaries [25].Figure 2
Computational domain defining thesi stretching functions.Thus, for a nonuniform mesh, equation (1) is applied to the si functions, as follows:(4)∇2si=C¯i;i=1,2,3,where(5)C¯i=−K¯1i+K¯2ig¯ii,and(6)K¯1i+K¯2i=1g¯g¯ii3/2Γ¯jjig¯kk−2Γ¯jkig¯jk+Γ¯kkig¯jj,and where the bar over the terms in equations (4) to (6) denotes that they are functions of si, not functions of the curvilinear coordinates ξi [11].The abovementioned relations can be recast in terms of the curvilinear coordinates by expressing the Laplacian ofsi in equation (4) as a function of the curvilinear coordinates, which is given by the following expression:(7)∇2si=gjk∂2si∂ξj∂ξk+∇2ξp∂si∂ξp,(8)=C¯i.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 ∇2ξp in equation (7) then leads to the following:(9)∇2ξp=Rp−gjk∂ξp∂si∂2si∂ξj∂ξk;p=1,2,3,where(10)Rp=C¯i∂ξp∂si.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, Rp, 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∇2r=0 (where r=x1x2x3T=xyzT) and expressing it in terms of the curvilinear coordinates:(11)∇2r=gjk∂2r∂ξj∂ξk+∇2ξp∂r∂ξp=0.The above system of equations is solved to yield the mesh coordinates, with the∇2ξp 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 ∇2ξp terms using equations (9) and (10) and the evaluation of the C¯i terms using equations (5) and (6).The solution process begins by defining thesi 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 si functions can be updated using the Neumann boundary conditions to enforce orthogonality conditions where desired, thereby modifying the si functions both on the boundaries and in the domain at every iteration.Thesi functions are expressed in terms of the curvilinear coordinates ξi, where s=s1s2s3T. Once the si functions have been defined, their first and second derivatives with respect to ξi 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 si functions, and these can be evaluated using the following well-known relations [26]:(12)∂ξi∂sj=1J^∂ss∂ξr∂sn∂ξm−∂ss∂ξm∂sn∂ξr,where i,r,m and j,s,n are cyclic, and where J^ is the determinant of the Jacobian matrix of the si-to-ξi coordinate transformation:(13)J^=det∂si∂ξj.To evaluate theRp terms in equation (10), the C¯i control functions must be evaluated with respect to the si coordinates. To do this, we must compute the metric tensor components and Christoffel symbols of the Cartesian xi-to-si coordinate transformation. The metric tensors of the transformation can be evaluated using the following well-known relations [27]:(14)g¯ij=rsi·rsj=∂xp∂si∂xp∂sj.(15)g¯=detg¯ij.(16)g¯ij=g¯rsg¯mn−g¯rng¯msg¯.where i,r,m and j,s,n are cyclic. To evaluate equation (14), the derivatives of the Cartesian coordinates with respect to the si coordinates are required, and these can be computed using the chain rule involving the ξi coordinates:(17)∂r∂sj=∂r∂ξp∂ξp∂sj,where the ∂ξp/∂sj are evaluated using equation (12).To evaluate the Christoffel symbols with respect to thesi coordinates in (6), denoted as Γ¯jki, the following expression can be used [27]:(18)Γ¯jki=g¯ilrsl·rsjsk,In turn, to evaluate the second derivative of the physical coordinates with respect to thesi coordinates, it is useful to express the second term on the right-hand side of equation (9) in the following compact form:(19)pjkp=−∂ξp∂si∂2si∂ξj∂ξk;p=1,2,3.Substituting (19) back into equation (9) then leads to the following:(20)∇2ξp=Rp+gjkpjkp;p=1,2,3.From equation (19) and the transformation laws for the Christoffel symbols [27], it can be shown that:(21)∂2r∂sj∂sk=∂2r∂ξm∂ξn+pmnp∂r∂ξp∂ξm∂sj∂ξn∂sk.Expressions (12) to (21) can be used to evaluate all terms required to construct the Rp 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ξ=ξ1 and η=ξ2 in equation (11), we can write(22)∇2r=g11rξξ+2g12rξη+g22rηη+∇2ξrξ+∇2ηrη=0.In two dimensions, the transformation relations between the covariant and contravariant metric tensors reduce to the following:(23)g11=g22g,g12=−g12g,g22=g11g.Substituting relations (23) into equation (22), we then obtain the following:(24)g22rξξ−2g12rξη+g11rηη+g∇2ξrξ+∇2ηrη=0.If we now combine equation (23) with equation (20) and expand, we obtain(25a)∇2ξ=R1+g22p111−2g12p121+g11p221g.(25b)∇2η=R2+g22p112−2g12p122+g11p222g.If we define the functionsP1 and P2 as follows, (26a)P1=g22p111−2g12p121+g11p221,(26b)P2=g22p112−2g12p122+g11p222,and insert equations (25) and (26) into equation (24), we obtain(27)g22rξξ−2g12rξη+g11rηη+P1+gR1rξ+P2+gR2rη=0.In the above equations, the components of the curvature control functions are represented by theR1 and R2 terms. If R1 and R2 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, g23=g13=g23=g13=0, g33=g33=1, and the expansion of the Christoffel symbols Γ¯jki given by equation (18), in combination with the relations given in (23), lead to the following expressions for the curvature control functions:(28a)C¯1=−K¯11+K¯21g¯11=−1g¯rs·rtt−g¯12g¯22rt·rtt,(28b)C¯2=−K¯12+K¯22g¯22=−1g¯rt·rss−g¯12g¯11rs·rss,where s,t=s1,s2.It can also be shown that(29)1g¯=g^g=J^2g.where g is the determinant of the metric tensor of the x,y⟶ξ,η transformation, g¯ is the determinant of the metric tensor of the x,y⟶s,t transformation, g^ is the determinant of the metric tensor of the s,t⟶ξ,η transformation, and J^ is the determinant of the Jacobian matrix of the s,t⟶ξ,η transformation (13).Inserting relations (29) into (28a)-(28b) then yields:(30a)gC¯1=−J^2rs·rtt−g¯12g¯22rt·rtt,(30b)gC¯2=−J^2rt·rss−g¯12g¯11rs·rss,where the first and second derivatives of r with respect to the s,t coordinates can be computed using equations (17) and (21), respectively. The latter equation can be simplified further by defining the function R as follows:(31)Rξiξj=rξiξj+pijprξp.Equation (21) can then be rewritten in a more compact form:(32)rsjsk=Rξmξnξsjmξskn.The curvature functions in equation (27) are then evaluated using (10):(33)gRp=gC¯iξsip.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:(34)rξ=ri+1,j−ri−1,j2Δξrη=ri,j+1−ri,j−12Δηrξξ=ri+1,j−2ri,j+ri−1,jΔξ2rηη=ri,j+1−2ri,j+ri,j−1Δη2rξη=ri+1,j+1−ri+1,j−1−ri−1,j+1+ri−1,j−14ΔξΔη.The updated values of the coordinates at iterationn+1, denoted as ri,j+, are expressed in terms of the old values at iteration n plus the increments from iteration n to n+1, denoted as follows:(35)ri,j+=ri,j+Δri,j.Substituting finite-difference expressions (34) and (35) into (27), and linearizing the equations by neglecting all second-order terms, OΔri,j2, leads to a linear system of equations that must be solved at every iteration. Although (27) is a vector equation, because the grid coordinates xi are coupled through the product of the metric tensor components and the Christoffel symbols, the xi 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 settingR1=0 and R2=0 in (27). Employing the notation described in Section 2.3, the updated value of equation (27) at iteration n+1 can be written as follows:(36)g22+rξξ+−2g12+rξη++g11+rηη++P1+rξ++P2+rη+=0.Expanding the above, we obtain the following:(37)g22+Δg22rξξ+Δrξξ−2g12+Δg12rξη+Δrξη+g11+Δg11rηη+Δrηη+P1+ΔP1rξ+Δrξ+P2+ΔP2rη+Δrη=0.Multiplying all terms in (37) and neglecting all higher-order terms, the equation reduces to(38)g22rξξ−2g12rξη+g22rηη+P1rξ+P2rη+g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη+Δg22rξξ−2Δg12rξη+Δg11rηη+ΔP1rξ+ΔP2rη=0.To simplify the subsequent derivations, the left-hand side of equation (38) is expressed as follows:(39)G+=G+ΔG1+ΔG2=0,where(40a)G=g22rξξ−2g12rξη+g11rηη+P1rξ+P2rη.(40b)ΔG1=g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη.(40c)ΔG2=Δg22rξξ−2Δg12rξη+Δg11rηη+ΔP1rξ+ΔP2rη.In equation (40c), all terms must be expressed as functions of the grid coordinates, rij. To this end, expressions for ΔP1 and ΔP2 are derived from equations (26a) and (26b) as follows:(41a)ΔP1=Δg22p111−2Δg12p121+Δg11p221.(41b)ΔP2=Δg22p112−2Δg12p122+Δg11p222.In the above expressions, thepijp terms represent the grid stretching functions, which are only functions of the s,t 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:(42)ΔG2=Δg22rξξ+p111rξ+p112rη−2Δg12rξη+p121rξ+p122rη+Δg11rηη+p221rξ+p222rη.If we now substitute equation (31) in equation (42), we can write,(43)ΔG2=Δg22Rξξ−2Δg12Rξη+Δg11Rηη.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:(44)g11=rξ·rξ,g12=rξ·rη,g22=rη·rη.Differentiating equation (44) then leads to the following:(45a)Δg11=2rξ·Δrξ=2xξΔxξ+yξΔyξΔg12=rξ·Δrη+Δrξ·rη.(45b)=xξΔxη+yξΔyη+xηΔxξ+yηΔyξ.(45c)Δg22=2rη·Δrη=2xηΔxη+yηΔyη.Substituting equations (45a)-(45c) into equation (43) and collecting terms, we get an expression for ΔG2 in terms of the increments in the derivatives of the grid coordinates:(46)ΔG2=2Rηηxξ−RξηxηΔxξ+Rξξxη−RξηxξΔxη+Rηηyξ−RξηyηΔyξ+Rξξyη−RξηyξΔyη.Inserting equations (46), (40a) and (40b) into (39), we then obtain the following:(47)g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη+2Rηηxξ−RξηxηΔxξ+2Rξξxη−RξηxξΔxη+2Rηηyξ−RξηyηΔyξ+2Rξξyη−RξηyξΔyη=−G.Equation (47) is a vector equation in which r=xyT, R=XYT, and G=GxGyT. To separate this equation into its x and y components, we note that from equation (31), we have the following:(48a)Xξiξj=xξiξj+pijpxξp.(48b)Yξiξj=yξiξj+pijpyξp.Therefore, thex and y components of equation (47) can be written as(49)g22Δxξξ−2g12Δxξη+g11Δxηη+P1+2Xηηxξ−XξηxηΔxξ+P2+2Xξξxη−XξηxξΔxη+2Xηηyξ−XξηyηΔyξ+2Xξξyη−XξηyξΔyη=−Gx.(50)g22Δyξξ−2g12Δyξη+g11Δyηη+2Yηηxξ−YξηxηΔxξ+2Yξξxη−YξηxξΔxη+P1+2Yηηyξ−YξηyηΔyξ+P2+2Yξξyη−YξηyξΔyη=−Gy.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 xi,j and yi,j, using the expressions given in (34). However, only the updates in the derivatives in these equations, i.e., the Δ terms, involve xi,j+ and yi,j+, while all other terms involve old values of xi,j and yi,j from the previous iteration. Thus, the discretized form of these equations leads to a system of linear equations that can be solved for xi,j+ and yi,j+ 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 n+1, denoted as gRp+, are expressed in terms of the old values at iteration n plus the increments from iteration n to n+1:(51)gRp+=gRp+ΔgRp.In equation (33), the stretching function derivative terms ξsip are held constant during the solution process, and thus ΔgRp is given by the following:(52)ΔgRp=ΔgC¯iξsip,where the gC¯i terms are given by equations (30a) and (30b). In the latter equations, J^ is the Jacobian of the s,t⟶ξ,η coordinate transformation, and is held constant during the solution process. However, all the other terms, involving derivatives of r with respect to the s,t coordinates, or metric tensor components of the x,y⟶s,t coordinate transformation, will vary at every iteration. Therefore, using the chain rule for differentiation, ΔgC¯1 can be expressed as(53)ΔgC¯1=−J^2Δrs·rtt+rs·Δrtt−Δg¯12g¯22rt·rtt−g¯12g¯22Δrt·rtt+rt·Δrtt,and ΔgC¯2 is given by the following:(54)ΔgC¯2=−J^2Δrt·rss+rt·Δrss−Δg¯12g¯11rs·rss−g¯12g¯11Δrs·rss+rs·Δrss.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 Δrsi. From equation (17), rsk can be written as follows:(55)rsi=rξmξsim.Assuming that the stretching function derivative termsξsim are held constant, Δrsi can thus be written as(56)Δrsi=Δrξmξsim.Therefore, expanding the vector dot product and collecting terms,Δrsi·rsj can be expressed as follows:(57)Δrsi·rsj=Δrξmξsim·rsj=Δxξmlξsimxsjl=ξsimxsjlΔxξml.Similarly,Δrsi·rsisj can be written as(58)Δrsi·rsisj=ξsimxsisjlΔxξml.To evaluate thersi·Δrsisj terms, we assume again that the stretching function derivatives are held constant. Therefore, using equation (32), Δrsisj can be expressed as follows:(59)Δrsisj=ΔRξmξnξsimξsjn.Similarly, evaluatingΔR using equation (31) and substituting it into equation (59), we obtain the following:(60)Δrsisj=Δrξmξn+pmnkΔrξkξsimξsjn.Equating the vector components on both sides of equation (60), we can write(61)Δxsisjl=Δxξmξnl+pmnkΔxξklξsimξsjn.Similarly, using equation (55), the components of rsi can be expressed as follows:(62)xsil=xξmlξsim.Therefore, taking the dot product ofrsi and Δrsisj and substituting the expressions for Δxsisjl and xsil given by equations (61) and (62), rsi·Δrsisj can be expressed as(63)rsi·Δrsisj=xsilΔxsisjl=xsilΔxξmξnl+pmnkΔxξklξsimξsjn=xsilξsimξsjnΔxξmξnl+xsilξsimξsjnpmnkΔxξkl.From the expression for the metric tensor in equation (14), Δg¯ij can be expressed as follows:(64)Δg¯ij=Δrsi·rsj+Δrsj·rsi.Therefore, the change in the covariant metric tensor components for thesi coordinates can be expressed as follows:(65a)Δg¯12=Δrs·rt+rs·Δrt.(65b)Δg¯11=2rs·Δrs.(65c)Δg¯22=2rt·Δrt.If we then expand the dot products in equation (64) and rearrange, Δg¯ij can be expressed as follows:(66)Δg¯ij=Δxsilxsjl+xsilΔxsjl=Δxξmlξsimxsjl+xsilΔxξmlξsjm=ξsimxsjl+ξsjmxsilΔxξml.From the chain rule applied to the differentiation of quotients, we can write the following:(67a)Δg¯12g¯22=Δg¯12g¯22−g¯12g¯22Δg¯22g¯22.(67b)Δg¯12g¯11=Δg¯12g¯11−g¯12g¯11Δg¯11g¯11.Substituting equation (66) into equations (67a)-(67b), we then have the following(68a)Δg¯12g¯22=1g¯22ξsmxtl+ξtmxslΔxξml−2g¯12g¯222ξtmxtlΔxξml.(68b)Δg¯12g¯11=1g¯11ξsmxtl+ξtmxslΔxξml−2g¯12g¯112ξsmxslΔxξml.Now, combining equations (58), (63), (67a) and (68a) into equation (53) and rearranging, we get the following expression for ΔgC¯1:(69)ΔgC¯1=−J2^ξsmxttlΔxξml+xslξtmξtnΔxξmξnl+xslξtmξtnpmnkΔxξkl−rt·rttg¯22ξsmxtl+ξtmxsl−2g¯12g¯22ξtmxtlΔxξml−g¯12g¯22ξtmxttlΔxξml−g¯12g¯22xtlξtmξtnΔxξmξnl+xtlξtmξtnpmnkΔxξkl.After rearranging and collecting terms, we can finally expressΔgC¯1 as follows:(70)ΔgC¯1=Almn1Δxξmξnl+Blk1Δxξkl,where(71)Almn1=−J2^ξtmξtnxsl−g¯12g¯22xtl,and(72)Blk1=−J2^xttlξsk−g¯12g¯22ξtk+ξtmξtnpmnkxsl−g¯12g¯22xtl−rt·rttg¯22ξtkxsl−2g¯12g¯22xtl+ξskxtl.Similarly, if we combine equations (58), (63), (67b) and (68b) into equation (54) and rearrange, we get the following expression for ΔgC¯2:(73)ΔgC¯2=−J2^ξtkxsslΔxξkl+xtlξsmξsnΔxξmξnl+xtlξsmξsnpmnkΔxξkl−rs·rssg¯11ξskxtl+ξtkxsl−2g¯12g¯11ξskxslΔxξkl−g¯12g¯11ξskxsslΔxξkl−g¯12g¯11xslξsmξsnΔxξmξnl+xslξsmξsnpmnkΔxξkl.Rearranging and collecting terms in (73), we can then write ΔgC¯2 as follows:(74)ΔgC¯2=Almn2Δxξmξnl+Blk2Δxξkl,where(75)Almn2=−J2^ξsmξsnxtl−g¯12g¯11xsl,and(76)Blk2=−J2^xsslξtk−g¯12g¯11ξsk+ξsmξsnpmnkxtl−g¯12g¯11xsl−rs·rssg¯11ξskxtl−2g¯12g¯11xsl+ξtkxsl.Using expressions forΔgC¯1 and ΔgC¯2 from equations (70) and (74), respectively, ΔgRp in equation (52) can now be expressed as(77)ΔgRp=ΔgC¯1ξsp+ΔgC¯2ξtp=Almn1ξsp+Almn2ξtpΔxξmξnl+Blk1ξsp+Blk2ξtpΔxξkl.We now have all the terms necessary to evaluategRp+ in equation (51).In Section2.4.1, the derivation of the linearized grid generation equations ((49) and (50)) excluded the curvature terms gRp. 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):(78)gR1rξ+R2rη+R1Δrξ+R2Δrη+ΔgR1rξ+ΔgR2rη.Separating equation (78) into its x and y components then yields the following expressions:(79a)gR1xξ+R2xη+R1Δxξ+R2Δxη+ΔgR1xξ+ΔgR2xη,and(79b)gR1yξ+R2yη+R1Δyξ+R2Δyη+ΔgR1yξ+ΔgR2yη.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 ΔgC¯1 and ΔgC¯2 given by equation (77), we obtain the following complete linearized grid generation equations with curvature control functions:(80)g22Δxξξ−2g12Δxξη+g11Δxξξ+P1+2Xηηxξ−Xξηxη+gR1Δxξ+P2+2Xξξxη−Xξηxξ+gR2Δxη+2Xηηyξ−XξηyηΔyξ+2Xξξyη−XξηyξΔyη+Almn1ξs+Almn2ξtxξΔxξmξnl+Blk1ξs+Blk2ξtxξΔxξkl+Almn1ηs+Almn2ηtxηΔxξmξnl+Blk1ηs+Blk2ηtxηΔxξkl=−Gx−gR1xξ−gR2xη.and(81)g22Δyξξ−2g12Δyξη+g11Δyξξ+2Yηηxξ−YξηxηΔxξ+2Yξξxη−YξηxξΔxη+P1+2Yηηyξ−Yξηyη+gR1Δyξ+P2+2Yξξyη−Yξηyξ+gR2Δyη+Almn1ξs+Almn2ξtyξΔxξmξnl+Blk1ξs+Blk2ξtyξΔxξkl+Almn1ηs+Almn2ηtyηΔxξmξnl+Blk1ηs+Blk2ηtyηΔxξkl=−Gy−gR1yξ−gR2yη.Equations (80) and (81) are the fully linearized equations that must be solved numerically at each grid point to compute the x and y 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 Δxξml and Δxξmξnl, leading to equations of the following form:(82)Cxij1Δxξ+Cxij2Δxη+Cxij3Δyξ+Cxij4Δyη+Cxij5Δxξξ+Cxij6Δxξη+Cxij7Δxηη+Cxij8Δyξξ+Cxij9Δyξη+Cxij10Δyηη=Rhxij.(83)Cyij1Δxξ+Cyij2Δxη+Cyij3Δyξ+Cyij4Δyη+Cyij5Δxξξ+Cyij6Δxξη+Cyij7Δxηη+Cyij8Δyξξ+Cyij9Δyξη+Cyij10Δyηη=Rhyij.where(84a)Cxij1=P1+2Xηηxξ−Xξηxη+gR1+B111xs+B112xt.(84b)Cxij2=P2+2Xξξxη−Xξηxξ+gR2+B121xs+B122xt.(84c)Cxij3=2Xηηyξ−Xξηyη+B211xs+B212xt.(84d)Cxij4=2Xξξyη−Xξηyξ+B221xs+B222xt.(84e)Cxij5=g22+A1111xs+A1112xt.(84f)Cxij6=−2g12+A1121xs+A1122xt.(84g)Cxij7=g11+A1221xs+A1222xt.(84h)Cxij8=A2111xs+A2112xt.(84i)Cxij9=A2121xs+A2122xt.(84j)Cxij10=A2221xs+A2222xt,and(85a)Cyij1=2Yηηxξ−Yξηxη+B111ys+B112yt.(85b)Cyij2=2Yξξxη−Yξηxξ+B121ys+B122yt.(85c)Cyij3=P1+2Yηηyξ−Yξηyη+gR1+B211ys+B212yt.(85d)Cyij4=P2+2Yξξyη−Yξηyξ+gR2+B221ys+B222yt.(85e)Cyij5=A1111ys+A1112yt.(85f)Cyij6=A1121ys+A1122yt.(85g)Cyij7=A1221ys+A1222yt.(85h)Cyij8=g22+A2111ys+A2112yt.(85i)Cyij9=−2g12+A2121ys+A2122yt.(85j)Cyij10=g11+A2221ys+A2222yt,and where(86a)Rhxij=−Gx−gR1xξ−gR2xη.(86b)Rhyij=−Gy−gR1yξ−gR2yη.In equations (84a)–(84j) to (86a)–(86b), the curvature control functions can be deactivated by setting all the Rp,Almnp and Blkp terms to zero. The numerical implementation of equations (82) and (83) is described in detail in the next section.


## 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ξi 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:(1)∇2ξi=Ci;i=1,2,3,where Ci represents the curvature control functions and is given by the following:(2)Ci=−K1i+K2igii,where(3)K1i+K2i=1ggii3/2Γjjigkk−2Γjkigjk+Γkkigjj,and where i,j,k are cyclic, with no sum on j and k.In the above expressions, thegij are the contravariant metric tensor components, and the Γjki are the Christoffel symbols of the second kind. The physical interpretation of equation (3) can be made by noting that, for a ξi=constant surface, K1i+K2i 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 sii=1,2,3. The si 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 si space represent the normalized arc length distributions along the boundaries of the physical space, while the points in the interior of the si space are interpolated algebraically from the boundaries [25].Figure 2
Computational domain defining thesi stretching functions.Thus, for a nonuniform mesh, equation (1) is applied to the si functions, as follows:(4)∇2si=C¯i;i=1,2,3,where(5)C¯i=−K¯1i+K¯2ig¯ii,and(6)K¯1i+K¯2i=1g¯g¯ii3/2Γ¯jjig¯kk−2Γ¯jkig¯jk+Γ¯kkig¯jj,and where the bar over the terms in equations (4) to (6) denotes that they are functions of si, not functions of the curvilinear coordinates ξi [11].The abovementioned relations can be recast in terms of the curvilinear coordinates by expressing the Laplacian ofsi in equation (4) as a function of the curvilinear coordinates, which is given by the following expression:(7)∇2si=gjk∂2si∂ξj∂ξk+∇2ξp∂si∂ξp,(8)=C¯i.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 ∇2ξp in equation (7) then leads to the following:(9)∇2ξp=Rp−gjk∂ξp∂si∂2si∂ξj∂ξk;p=1,2,3,where(10)Rp=C¯i∂ξp∂si.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, Rp, 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∇2r=0 (where r=x1x2x3T=xyzT) and expressing it in terms of the curvilinear coordinates:(11)∇2r=gjk∂2r∂ξj∂ξk+∇2ξp∂r∂ξp=0.The above system of equations is solved to yield the mesh coordinates, with the∇2ξp 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 ∇2ξp terms using equations (9) and (10) and the evaluation of the C¯i terms using equations (5) and (6).The solution process begins by defining thesi 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 si functions can be updated using the Neumann boundary conditions to enforce orthogonality conditions where desired, thereby modifying the si functions both on the boundaries and in the domain at every iteration.Thesi functions are expressed in terms of the curvilinear coordinates ξi, where s=s1s2s3T. Once the si functions have been defined, their first and second derivatives with respect to ξi 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 si functions, and these can be evaluated using the following well-known relations [26]:(12)∂ξi∂sj=1J^∂ss∂ξr∂sn∂ξm−∂ss∂ξm∂sn∂ξr,where i,r,m and j,s,n are cyclic, and where J^ is the determinant of the Jacobian matrix of the si-to-ξi coordinate transformation:(13)J^=det∂si∂ξj.To evaluate theRp terms in equation (10), the C¯i control functions must be evaluated with respect to the si coordinates. To do this, we must compute the metric tensor components and Christoffel symbols of the Cartesian xi-to-si coordinate transformation. The metric tensors of the transformation can be evaluated using the following well-known relations [27]:(14)g¯ij=rsi·rsj=∂xp∂si∂xp∂sj.(15)g¯=detg¯ij.(16)g¯ij=g¯rsg¯mn−g¯rng¯msg¯.where i,r,m and j,s,n are cyclic. To evaluate equation (14), the derivatives of the Cartesian coordinates with respect to the si coordinates are required, and these can be computed using the chain rule involving the ξi coordinates:(17)∂r∂sj=∂r∂ξp∂ξp∂sj,where the ∂ξp/∂sj are evaluated using equation (12).To evaluate the Christoffel symbols with respect to thesi coordinates in (6), denoted as Γ¯jki, the following expression can be used [27]:(18)Γ¯jki=g¯ilrsl·rsjsk,In turn, to evaluate the second derivative of the physical coordinates with respect to thesi coordinates, it is useful to express the second term on the right-hand side of equation (9) in the following compact form:(19)pjkp=−∂ξp∂si∂2si∂ξj∂ξk;p=1,2,3.Substituting (19) back into equation (9) then leads to the following:(20)∇2ξp=Rp+gjkpjkp;p=1,2,3.From equation (19) and the transformation laws for the Christoffel symbols [27], it can be shown that:(21)∂2r∂sj∂sk=∂2r∂ξm∂ξn+pmnp∂r∂ξp∂ξm∂sj∂ξn∂sk.Expressions (12) to (21) can be used to evaluate all terms required to construct the Rp curvature functions appearing in equations (9) and (10), and the grid generation equation (11) can then be solved numerically to generate the mesh.


## 2.1.1. Mesh Generation Equations

Denotingξi 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:(1)∇2ξi=Ci;i=1,2,3,where Ci represents the curvature control functions and is given by the following:(2)Ci=−K1i+K2igii,where(3)K1i+K2i=1ggii3/2Γjjigkk−2Γjkigjk+Γkkigjj,and where i,j,k are cyclic, with no sum on j and k.In the above expressions, thegij are the contravariant metric tensor components, and the Γjki are the Christoffel symbols of the second kind. The physical interpretation of equation (3) can be made by noting that, for a ξi=constant surface, K1i+K2i 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 sii=1,2,3. The si 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 si space represent the normalized arc length distributions along the boundaries of the physical space, while the points in the interior of the si space are interpolated algebraically from the boundaries [25].Figure 2
Computational domain defining thesi stretching functions.Thus, for a nonuniform mesh, equation (1) is applied to the si functions, as follows:(4)∇2si=C¯i;i=1,2,3,where(5)C¯i=−K¯1i+K¯2ig¯ii,and(6)K¯1i+K¯2i=1g¯g¯ii3/2Γ¯jjig¯kk−2Γ¯jkig¯jk+Γ¯kkig¯jj,and where the bar over the terms in equations (4) to (6) denotes that they are functions of si, not functions of the curvilinear coordinates ξi [11].The abovementioned relations can be recast in terms of the curvilinear coordinates by expressing the Laplacian ofsi in equation (4) as a function of the curvilinear coordinates, which is given by the following expression:(7)∇2si=gjk∂2si∂ξj∂ξk+∇2ξp∂si∂ξp,(8)=C¯i.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 ∇2ξp in equation (7) then leads to the following:(9)∇2ξp=Rp−gjk∂ξp∂si∂2si∂ξj∂ξk;p=1,2,3,where(10)Rp=C¯i∂ξp∂si.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, Rp, 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∇2r=0 (where r=x1x2x3T=xyzT) and expressing it in terms of the curvilinear coordinates:(11)∇2r=gjk∂2r∂ξj∂ξk+∇2ξp∂r∂ξp=0.The above system of equations is solved to yield the mesh coordinates, with the∇2ξp 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 ∇2ξp terms using equations (9) and (10) and the evaluation of the C¯i terms using equations (5) and (6).The solution process begins by defining thesi 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 si functions can be updated using the Neumann boundary conditions to enforce orthogonality conditions where desired, thereby modifying the si functions both on the boundaries and in the domain at every iteration.Thesi functions are expressed in terms of the curvilinear coordinates ξi, where s=s1s2s3T. Once the si functions have been defined, their first and second derivatives with respect to ξi 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 si functions, and these can be evaluated using the following well-known relations [26]:(12)∂ξi∂sj=1J^∂ss∂ξr∂sn∂ξm−∂ss∂ξm∂sn∂ξr,where i,r,m and j,s,n are cyclic, and where J^ is the determinant of the Jacobian matrix of the si-to-ξi coordinate transformation:(13)J^=det∂si∂ξj.To evaluate theRp terms in equation (10), the C¯i control functions must be evaluated with respect to the si coordinates. To do this, we must compute the metric tensor components and Christoffel symbols of the Cartesian xi-to-si coordinate transformation. The metric tensors of the transformation can be evaluated using the following well-known relations [27]:(14)g¯ij=rsi·rsj=∂xp∂si∂xp∂sj.(15)g¯=detg¯ij.(16)g¯ij=g¯rsg¯mn−g¯rng¯msg¯.where i,r,m and j,s,n are cyclic. To evaluate equation (14), the derivatives of the Cartesian coordinates with respect to the si coordinates are required, and these can be computed using the chain rule involving the ξi coordinates:(17)∂r∂sj=∂r∂ξp∂ξp∂sj,where the ∂ξp/∂sj are evaluated using equation (12).To evaluate the Christoffel symbols with respect to thesi coordinates in (6), denoted as Γ¯jki, the following expression can be used [27]:(18)Γ¯jki=g¯ilrsl·rsjsk,In turn, to evaluate the second derivative of the physical coordinates with respect to thesi coordinates, it is useful to express the second term on the right-hand side of equation (9) in the following compact form:(19)pjkp=−∂ξp∂si∂2si∂ξj∂ξk;p=1,2,3.Substituting (19) back into equation (9) then leads to the following:(20)∇2ξp=Rp+gjkpjkp;p=1,2,3.From equation (19) and the transformation laws for the Christoffel symbols [27], it can be shown that:(21)∂2r∂sj∂sk=∂2r∂ξm∂ξn+pmnp∂r∂ξp∂ξm∂sj∂ξn∂sk.Expressions (12) to (21) can be used to evaluate all terms required to construct the Rp 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ξ=ξ1 and η=ξ2 in equation (11), we can write(22)∇2r=g11rξξ+2g12rξη+g22rηη+∇2ξrξ+∇2ηrη=0.In two dimensions, the transformation relations between the covariant and contravariant metric tensors reduce to the following:(23)g11=g22g,g12=−g12g,g22=g11g.Substituting relations (23) into equation (22), we then obtain the following:(24)g22rξξ−2g12rξη+g11rηη+g∇2ξrξ+∇2ηrη=0.If we now combine equation (23) with equation (20) and expand, we obtain(25a)∇2ξ=R1+g22p111−2g12p121+g11p221g.(25b)∇2η=R2+g22p112−2g12p122+g11p222g.If we define the functionsP1 and P2 as follows, (26a)P1=g22p111−2g12p121+g11p221,(26b)P2=g22p112−2g12p122+g11p222,and insert equations (25) and (26) into equation (24), we obtain(27)g22rξξ−2g12rξη+g11rηη+P1+gR1rξ+P2+gR2rη=0.In the above equations, the components of the curvature control functions are represented by theR1 and R2 terms. If R1 and R2 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, g23=g13=g23=g13=0, g33=g33=1, and the expansion of the Christoffel symbols Γ¯jki given by equation (18), in combination with the relations given in (23), lead to the following expressions for the curvature control functions:(28a)C¯1=−K¯11+K¯21g¯11=−1g¯rs·rtt−g¯12g¯22rt·rtt,(28b)C¯2=−K¯12+K¯22g¯22=−1g¯rt·rss−g¯12g¯11rs·rss,where s,t=s1,s2.It can also be shown that(29)1g¯=g^g=J^2g.where g is the determinant of the metric tensor of the x,y⟶ξ,η transformation, g¯ is the determinant of the metric tensor of the x,y⟶s,t transformation, g^ is the determinant of the metric tensor of the s,t⟶ξ,η transformation, and J^ is the determinant of the Jacobian matrix of the s,t⟶ξ,η transformation (13).Inserting relations (29) into (28a)-(28b) then yields:(30a)gC¯1=−J^2rs·rtt−g¯12g¯22rt·rtt,(30b)gC¯2=−J^2rt·rss−g¯12g¯11rs·rss,where the first and second derivatives of r with respect to the s,t coordinates can be computed using equations (17) and (21), respectively. The latter equation can be simplified further by defining the function R as follows:(31)Rξiξj=rξiξj+pijprξp.Equation (21) can then be rewritten in a more compact form:(32)rsjsk=Rξmξnξsjmξskn.The curvature functions in equation (27) are then evaluated using (10):(33)gRp=gC¯iξsip.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.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, g23=g13=g23=g13=0, g33=g33=1, and the expansion of the Christoffel symbols Γ¯jki given by equation (18), in combination with the relations given in (23), lead to the following expressions for the curvature control functions:(28a)C¯1=−K¯11+K¯21g¯11=−1g¯rs·rtt−g¯12g¯22rt·rtt,(28b)C¯2=−K¯12+K¯22g¯22=−1g¯rt·rss−g¯12g¯11rs·rss,where s,t=s1,s2.It can also be shown that(29)1g¯=g^g=J^2g.where g is the determinant of the metric tensor of the x,y⟶ξ,η transformation, g¯ is the determinant of the metric tensor of the x,y⟶s,t transformation, g^ is the determinant of the metric tensor of the s,t⟶ξ,η transformation, and J^ is the determinant of the Jacobian matrix of the s,t⟶ξ,η transformation (13).Inserting relations (29) into (28a)-(28b) then yields:(30a)gC¯1=−J^2rs·rtt−g¯12g¯22rt·rtt,(30b)gC¯2=−J^2rt·rss−g¯12g¯11rs·rss,where the first and second derivatives of r with respect to the s,t coordinates can be computed using equations (17) and (21), respectively. The latter equation can be simplified further by defining the function R as follows:(31)Rξiξj=rξiξj+pijprξp.Equation (21) can then be rewritten in a more compact form:(32)rsjsk=Rξmξnξsjmξskn.The curvature functions in equation (27) are then evaluated using (10):(33)gRp=gC¯iξsip.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:(34)rξ=ri+1,j−ri−1,j2Δξrη=ri,j+1−ri,j−12Δηrξξ=ri+1,j−2ri,j+ri−1,jΔξ2rηη=ri,j+1−2ri,j+ri,j−1Δη2rξη=ri+1,j+1−ri+1,j−1−ri−1,j+1+ri−1,j−14ΔξΔη.The updated values of the coordinates at iterationn+1, denoted as ri,j+, are expressed in terms of the old values at iteration n plus the increments from iteration n to n+1, denoted as follows:(35)ri,j+=ri,j+Δri,j.Substituting finite-difference expressions (34) and (35) into (27), and linearizing the equations by neglecting all second-order terms, OΔri,j2, leads to a linear system of equations that must be solved at every iteration. Although (27) is a vector equation, because the grid coordinates xi are coupled through the product of the metric tensor components and the Christoffel symbols, the xi 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 settingR1=0 and R2=0 in (27). Employing the notation described in Section 2.3, the updated value of equation (27) at iteration n+1 can be written as follows:(36)g22+rξξ+−2g12+rξη++g11+rηη++P1+rξ++P2+rη+=0.Expanding the above, we obtain the following:(37)g22+Δg22rξξ+Δrξξ−2g12+Δg12rξη+Δrξη+g11+Δg11rηη+Δrηη+P1+ΔP1rξ+Δrξ+P2+ΔP2rη+Δrη=0.Multiplying all terms in (37) and neglecting all higher-order terms, the equation reduces to(38)g22rξξ−2g12rξη+g22rηη+P1rξ+P2rη+g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη+Δg22rξξ−2Δg12rξη+Δg11rηη+ΔP1rξ+ΔP2rη=0.To simplify the subsequent derivations, the left-hand side of equation (38) is expressed as follows:(39)G+=G+ΔG1+ΔG2=0,where(40a)G=g22rξξ−2g12rξη+g11rηη+P1rξ+P2rη.(40b)ΔG1=g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη.(40c)ΔG2=Δg22rξξ−2Δg12rξη+Δg11rηη+ΔP1rξ+ΔP2rη.In equation (40c), all terms must be expressed as functions of the grid coordinates, rij. To this end, expressions for ΔP1 and ΔP2 are derived from equations (26a) and (26b) as follows:(41a)ΔP1=Δg22p111−2Δg12p121+Δg11p221.(41b)ΔP2=Δg22p112−2Δg12p122+Δg11p222.In the above expressions, thepijp terms represent the grid stretching functions, which are only functions of the s,t 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:(42)ΔG2=Δg22rξξ+p111rξ+p112rη−2Δg12rξη+p121rξ+p122rη+Δg11rηη+p221rξ+p222rη.If we now substitute equation (31) in equation (42), we can write,(43)ΔG2=Δg22Rξξ−2Δg12Rξη+Δg11Rηη.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:(44)g11=rξ·rξ,g12=rξ·rη,g22=rη·rη.Differentiating equation (44) then leads to the following:(45a)Δg11=2rξ·Δrξ=2xξΔxξ+yξΔyξΔg12=rξ·Δrη+Δrξ·rη.(45b)=xξΔxη+yξΔyη+xηΔxξ+yηΔyξ.(45c)Δg22=2rη·Δrη=2xηΔxη+yηΔyη.Substituting equations (45a)-(45c) into equation (43) and collecting terms, we get an expression for ΔG2 in terms of the increments in the derivatives of the grid coordinates:(46)ΔG2=2Rηηxξ−RξηxηΔxξ+Rξξxη−RξηxξΔxη+Rηηyξ−RξηyηΔyξ+Rξξyη−RξηyξΔyη.Inserting equations (46), (40a) and (40b) into (39), we then obtain the following:(47)g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη+2Rηηxξ−RξηxηΔxξ+2Rξξxη−RξηxξΔxη+2Rηηyξ−RξηyηΔyξ+2Rξξyη−RξηyξΔyη=−G.Equation (47) is a vector equation in which r=xyT, R=XYT, and G=GxGyT. To separate this equation into its x and y components, we note that from equation (31), we have the following:(48a)Xξiξj=xξiξj+pijpxξp.(48b)Yξiξj=yξiξj+pijpyξp.Therefore, thex and y components of equation (47) can be written as(49)g22Δxξξ−2g12Δxξη+g11Δxηη+P1+2Xηηxξ−XξηxηΔxξ+P2+2Xξξxη−XξηxξΔxη+2Xηηyξ−XξηyηΔyξ+2Xξξyη−XξηyξΔyη=−Gx.(50)g22Δyξξ−2g12Δyξη+g11Δyηη+2Yηηxξ−YξηxηΔxξ+2Yξξxη−YξηxξΔxη+P1+2Yηηyξ−YξηyηΔyξ+P2+2Yξξyη−YξηyξΔyη=−Gy.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 xi,j and yi,j, using the expressions given in (34). However, only the updates in the derivatives in these equations, i.e., the Δ terms, involve xi,j+ and yi,j+, while all other terms involve old values of xi,j and yi,j from the previous iteration. Thus, the discretized form of these equations leads to a system of linear equations that can be solved for xi,j+ and yi,j+ 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 n+1, denoted as gRp+, are expressed in terms of the old values at iteration n plus the increments from iteration n to n+1:(51)gRp+=gRp+ΔgRp.In equation (33), the stretching function derivative terms ξsip are held constant during the solution process, and thus ΔgRp is given by the following:(52)ΔgRp=ΔgC¯iξsip,where the gC¯i terms are given by equations (30a) and (30b). In the latter equations, J^ is the Jacobian of the s,t⟶ξ,η coordinate transformation, and is held constant during the solution process. However, all the other terms, involving derivatives of r with respect to the s,t coordinates, or metric tensor components of the x,y⟶s,t coordinate transformation, will vary at every iteration. Therefore, using the chain rule for differentiation, ΔgC¯1 can be expressed as(53)ΔgC¯1=−J^2Δrs·rtt+rs·Δrtt−Δg¯12g¯22rt·rtt−g¯12g¯22Δrt·rtt+rt·Δrtt,and ΔgC¯2 is given by the following:(54)ΔgC¯2=−J^2Δrt·rss+rt·Δrss−Δg¯12g¯11rs·rss−g¯12g¯11Δrs·rss+rs·Δrss.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 Δrsi. From equation (17), rsk can be written as follows:(55)rsi=rξmξsim.Assuming that the stretching function derivative termsξsim are held constant, Δrsi can thus be written as(56)Δrsi=Δrξmξsim.Therefore, expanding the vector dot product and collecting terms,Δrsi·rsj can be expressed as follows:(57)Δrsi·rsj=Δrξmξsim·rsj=Δxξmlξsimxsjl=ξsimxsjlΔxξml.Similarly,Δrsi·rsisj can be written as(58)Δrsi·rsisj=ξsimxsisjlΔxξml.To evaluate thersi·Δrsisj terms, we assume again that the stretching function derivatives are held constant. Therefore, using equation (32), Δrsisj can be expressed as follows:(59)Δrsisj=ΔRξmξnξsimξsjn.Similarly, evaluatingΔR using equation (31) and substituting it into equation (59), we obtain the following:(60)Δrsisj=Δrξmξn+pmnkΔrξkξsimξsjn.Equating the vector components on both sides of equation (60), we can write(61)Δxsisjl=Δxξmξnl+pmnkΔxξklξsimξsjn.Similarly, using equation (55), the components of rsi can be expressed as follows:(62)xsil=xξmlξsim.Therefore, taking the dot product ofrsi and Δrsisj and substituting the expressions for Δxsisjl and xsil given by equations (61) and (62), rsi·Δrsisj can be expressed as(63)rsi·Δrsisj=xsilΔxsisjl=xsilΔxξmξnl+pmnkΔxξklξsimξsjn=xsilξsimξsjnΔxξmξnl+xsilξsimξsjnpmnkΔxξkl.From the expression for the metric tensor in equation (14), Δg¯ij can be expressed as follows:(64)Δg¯ij=Δrsi·rsj+Δrsj·rsi.Therefore, the change in the covariant metric tensor components for thesi coordinates can be expressed as follows:(65a)Δg¯12=Δrs·rt+rs·Δrt.(65b)Δg¯11=2rs·Δrs.(65c)Δg¯22=2rt·Δrt.If we then expand the dot products in equation (64) and rearrange, Δg¯ij can be expressed as follows:(66)Δg¯ij=Δxsilxsjl+xsilΔxsjl=Δxξmlξsimxsjl+xsilΔxξmlξsjm=ξsimxsjl+ξsjmxsilΔxξml.From the chain rule applied to the differentiation of quotients, we can write the following:(67a)Δg¯12g¯22=Δg¯12g¯22−g¯12g¯22Δg¯22g¯22.(67b)Δg¯12g¯11=Δg¯12g¯11−g¯12g¯11Δg¯11g¯11.Substituting equation (66) into equations (67a)-(67b), we then have the following(68a)Δg¯12g¯22=1g¯22ξsmxtl+ξtmxslΔxξml−2g¯12g¯222ξtmxtlΔxξml.(68b)Δg¯12g¯11=1g¯11ξsmxtl+ξtmxslΔxξml−2g¯12g¯112ξsmxslΔxξml.Now, combining equations (58), (63), (67a) and (68a) into equation (53) and rearranging, we get the following expression for ΔgC¯1:(69)ΔgC¯1=−J2^ξsmxttlΔxξml+xslξtmξtnΔxξmξnl+xslξtmξtnpmnkΔxξkl−rt·rttg¯22ξsmxtl+ξtmxsl−2g¯12g¯22ξtmxtlΔxξml−g¯12g¯22ξtmxttlΔxξml−g¯12g¯22xtlξtmξtnΔxξmξnl+xtlξtmξtnpmnkΔxξkl.After rearranging and collecting terms, we can finally expressΔgC¯1 as follows:(70)ΔgC¯1=Almn1Δxξmξnl+Blk1Δxξkl,where(71)Almn1=−J2^ξtmξtnxsl−g¯12g¯22xtl,and(72)Blk1=−J2^xttlξsk−g¯12g¯22ξtk+ξtmξtnpmnkxsl−g¯12g¯22xtl−rt·rttg¯22ξtkxsl−2g¯12g¯22xtl+ξskxtl.Similarly, if we combine equations (58), (63), (67b) and (68b) into equation (54) and rearrange, we get the following expression for ΔgC¯2:(73)ΔgC¯2=−J2^ξtkxsslΔxξkl+xtlξsmξsnΔxξmξnl+xtlξsmξsnpmnkΔxξkl−rs·rssg¯11ξskxtl+ξtkxsl−2g¯12g¯11ξskxslΔxξkl−g¯12g¯11ξskxsslΔxξkl−g¯12g¯11xslξsmξsnΔxξmξnl+xslξsmξsnpmnkΔxξkl.Rearranging and collecting terms in (73), we can then write ΔgC¯2 as follows:(74)ΔgC¯2=Almn2Δxξmξnl+Blk2Δxξkl,where(75)Almn2=−J2^ξsmξsnxtl−g¯12g¯11xsl,and(76)Blk2=−J2^xsslξtk−g¯12g¯11ξsk+ξsmξsnpmnkxtl−g¯12g¯11xsl−rs·rssg¯11ξskxtl−2g¯12g¯11xsl+ξtkxsl.Using expressions forΔgC¯1 and ΔgC¯2 from equations (70) and (74), respectively, ΔgRp in equation (52) can now be expressed as(77)ΔgRp=ΔgC¯1ξsp+ΔgC¯2ξtp=Almn1ξsp+Almn2ξtpΔxξmξnl+Blk1ξsp+Blk2ξtpΔxξkl.We now have all the terms necessary to evaluategRp+ in equation (51).In Section2.4.1, the derivation of the linearized grid generation equations ((49) and (50)) excluded the curvature terms gRp. 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):(78)gR1rξ+R2rη+R1Δrξ+R2Δrη+ΔgR1rξ+ΔgR2rη.Separating equation (78) into its x and y components then yields the following expressions:(79a)gR1xξ+R2xη+R1Δxξ+R2Δxη+ΔgR1xξ+ΔgR2xη,and(79b)gR1yξ+R2yη+R1Δyξ+R2Δyη+ΔgR1yξ+ΔgR2yη.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 ΔgC¯1 and ΔgC¯2 given by equation (77), we obtain the following complete linearized grid generation equations with curvature control functions:(80)g22Δxξξ−2g12Δxξη+g11Δxξξ+P1+2Xηηxξ−Xξηxη+gR1Δxξ+P2+2Xξξxη−Xξηxξ+gR2Δxη+2Xηηyξ−XξηyηΔyξ+2Xξξyη−XξηyξΔyη+Almn1ξs+Almn2ξtxξΔxξmξnl+Blk1ξs+Blk2ξtxξΔxξkl+Almn1ηs+Almn2ηtxηΔxξmξnl+Blk1ηs+Blk2ηtxηΔxξkl=−Gx−gR1xξ−gR2xη.and(81)g22Δyξξ−2g12Δyξη+g11Δyξξ+2Yηηxξ−YξηxηΔxξ+2Yξξxη−YξηxξΔxη+P1+2Yηηyξ−Yξηyη+gR1Δyξ+P2+2Yξξyη−Yξηyξ+gR2Δyη+Almn1ξs+Almn2ξtyξΔxξmξnl+Blk1ξs+Blk2ξtyξΔxξkl+Almn1ηs+Almn2ηtyηΔxξmξnl+Blk1ηs+Blk2ηtyηΔxξkl=−Gy−gR1yξ−gR2yη.Equations (80) and (81) are the fully linearized equations that must be solved numerically at each grid point to compute the x and y 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 Δxξml and Δxξmξnl, leading to equations of the following form:(82)Cxij1Δxξ+Cxij2Δxη+Cxij3Δyξ+Cxij4Δyη+Cxij5Δxξξ+Cxij6Δxξη+Cxij7Δxηη+Cxij8Δyξξ+Cxij9Δyξη+Cxij10Δyηη=Rhxij.(83)Cyij1Δxξ+Cyij2Δxη+Cyij3Δyξ+Cyij4Δyη+Cyij5Δxξξ+Cyij6Δxξη+Cyij7Δxηη+Cyij8Δyξξ+Cyij9Δyξη+Cyij10Δyηη=Rhyij.where(84a)Cxij1=P1+2Xηηxξ−Xξηxη+gR1+B111xs+B112xt.(84b)Cxij2=P2+2Xξξxη−Xξηxξ+gR2+B121xs+B122xt.(84c)Cxij3=2Xηηyξ−Xξηyη+B211xs+B212xt.(84d)Cxij4=2Xξξyη−Xξηyξ+B221xs+B222xt.(84e)Cxij5=g22+A1111xs+A1112xt.(84f)Cxij6=−2g12+A1121xs+A1122xt.(84g)Cxij7=g11+A1221xs+A1222xt.(84h)Cxij8=A2111xs+A2112xt.(84i)Cxij9=A2121xs+A2122xt.(84j)Cxij10=A2221xs+A2222xt,and(85a)Cyij1=2Yηηxξ−Yξηxη+B111ys+B112yt.(85b)Cyij2=2Yξξxη−Yξηxξ+B121ys+B122yt.(85c)Cyij3=P1+2Yηηyξ−Yξηyη+gR1+B211ys+B212yt.(85d)Cyij4=P2+2Yξξyη−Yξηyξ+gR2+B221ys+B222yt.(85e)Cyij5=A1111ys+A1112yt.(85f)Cyij6=A1121ys+A1122yt.(85g)Cyij7=A1221ys+A1222yt.(85h)Cyij8=g22+A2111ys+A2112yt.(85i)Cyij9=−2g12+A2121ys+A2122yt.(85j)Cyij10=g11+A2221ys+A2222yt,and where(86a)Rhxij=−Gx−gR1xξ−gR2xη.(86b)Rhyij=−Gy−gR1yξ−gR2yη.In equations (84a)–(84j) to (86a)–(86b), the curvature control functions can be deactivated by setting all the Rp,Almnp and Blkp terms to zero. The numerical implementation of equations (82) and (83) is described in detail in the next section.


## 2.4.1. Linearization of Equations without Curvature Functions

The grid generation equations without the curvature functions are obtained by settingR1=0 and R2=0 in (27). Employing the notation described in Section 2.3, the updated value of equation (27) at iteration n+1 can be written as follows:(36)g22+rξξ+−2g12+rξη++g11+rηη++P1+rξ++P2+rη+=0.Expanding the above, we obtain the following:(37)g22+Δg22rξξ+Δrξξ−2g12+Δg12rξη+Δrξη+g11+Δg11rηη+Δrηη+P1+ΔP1rξ+Δrξ+P2+ΔP2rη+Δrη=0.Multiplying all terms in (37) and neglecting all higher-order terms, the equation reduces to(38)g22rξξ−2g12rξη+g22rηη+P1rξ+P2rη+g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη+Δg22rξξ−2Δg12rξη+Δg11rηη+ΔP1rξ+ΔP2rη=0.To simplify the subsequent derivations, the left-hand side of equation (38) is expressed as follows:(39)G+=G+ΔG1+ΔG2=0,where(40a)G=g22rξξ−2g12rξη+g11rηη+P1rξ+P2rη.(40b)ΔG1=g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη.(40c)ΔG2=Δg22rξξ−2Δg12rξη+Δg11rηη+ΔP1rξ+ΔP2rη.In equation (40c), all terms must be expressed as functions of the grid coordinates, rij. To this end, expressions for ΔP1 and ΔP2 are derived from equations (26a) and (26b) as follows:(41a)ΔP1=Δg22p111−2Δg12p121+Δg11p221.(41b)ΔP2=Δg22p112−2Δg12p122+Δg11p222.In the above expressions, thepijp terms represent the grid stretching functions, which are only functions of the s,t 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:(42)ΔG2=Δg22rξξ+p111rξ+p112rη−2Δg12rξη+p121rξ+p122rη+Δg11rηη+p221rξ+p222rη.If we now substitute equation (31) in equation (42), we can write,(43)ΔG2=Δg22Rξξ−2Δg12Rξη+Δg11Rηη.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:(44)g11=rξ·rξ,g12=rξ·rη,g22=rη·rη.Differentiating equation (44) then leads to the following:(45a)Δg11=2rξ·Δrξ=2xξΔxξ+yξΔyξΔg12=rξ·Δrη+Δrξ·rη.(45b)=xξΔxη+yξΔyη+xηΔxξ+yηΔyξ.(45c)Δg22=2rη·Δrη=2xηΔxη+yηΔyη.Substituting equations (45a)-(45c) into equation (43) and collecting terms, we get an expression for ΔG2 in terms of the increments in the derivatives of the grid coordinates:(46)ΔG2=2Rηηxξ−RξηxηΔxξ+Rξξxη−RξηxξΔxη+Rηηyξ−RξηyηΔyξ+Rξξyη−RξηyξΔyη.Inserting equations (46), (40a) and (40b) into (39), we then obtain the following:(47)g22Δrξξ−2g12Δrξη+g11Δrηη+P1Δrξ+P2Δrη+2Rηηxξ−RξηxηΔxξ+2Rξξxη−RξηxξΔxη+2Rηηyξ−RξηyηΔyξ+2Rξξyη−RξηyξΔyη=−G.Equation (47) is a vector equation in which r=xyT, R=XYT, and G=GxGyT. To separate this equation into its x and y components, we note that from equation (31), we have the following:(48a)Xξiξj=xξiξj+pijpxξp.(48b)Yξiξj=yξiξj+pijpyξp.Therefore, thex and y components of equation (47) can be written as(49)g22Δxξξ−2g12Δxξη+g11Δxηη+P1+2Xηηxξ−XξηxηΔxξ+P2+2Xξξxη−XξηxξΔxη+2Xηηyξ−XξηyηΔyξ+2Xξξyη−XξηyξΔyη=−Gx.(50)g22Δyξξ−2g12Δyξη+g11Δyηη+2Yηηxξ−YξηxηΔxξ+2Yξξxη−YξηxξΔxη+P1+2Yηηyξ−YξηyηΔyξ+P2+2Yξξyη−YξηyξΔyη=−Gy.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 xi,j and yi,j, using the expressions given in (34). However, only the updates in the derivatives in these equations, i.e., the Δ terms, involve xi,j+ and yi,j+, while all other terms involve old values of xi,j and yi,j from the previous iteration. Thus, the discretized form of these equations leads to a system of linear equations that can be solved for xi,j+ and yi,j+ 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 n+1, denoted as gRp+, are expressed in terms of the old values at iteration n plus the increments from iteration n to n+1:(51)gRp+=gRp+ΔgRp.In equation (33), the stretching function derivative terms ξsip are held constant during the solution process, and thus ΔgRp is given by the following:(52)ΔgRp=ΔgC¯iξsip,where the gC¯i terms are given by equations (30a) and (30b). In the latter equations, J^ is the Jacobian of the s,t⟶ξ,η coordinate transformation, and is held constant during the solution process. However, all the other terms, involving derivatives of r with respect to the s,t coordinates, or metric tensor components of the x,y⟶s,t coordinate transformation, will vary at every iteration. Therefore, using the chain rule for differentiation, ΔgC¯1 can be expressed as(53)ΔgC¯1=−J^2Δrs·rtt+rs·Δrtt−Δg¯12g¯22rt·rtt−g¯12g¯22Δrt·rtt+rt·Δrtt,and ΔgC¯2 is given by the following:(54)ΔgC¯2=−J^2Δrt·rss+rt·Δrss−Δg¯12g¯11rs·rss−g¯12g¯11Δrs·rss+rs·Δrss.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 Δrsi. From equation (17), rsk can be written as follows:(55)rsi=rξmξsim.Assuming that the stretching function derivative termsξsim are held constant, Δrsi can thus be written as(56)Δrsi=Δrξmξsim.Therefore, expanding the vector dot product and collecting terms,Δrsi·rsj can be expressed as follows:(57)Δrsi·rsj=Δrξmξsim·rsj=Δxξmlξsimxsjl=ξsimxsjlΔxξml.Similarly,Δrsi·rsisj can be written as(58)Δrsi·rsisj=ξsimxsisjlΔxξml.To evaluate thersi·Δrsisj terms, we assume again that the stretching function derivatives are held constant. Therefore, using equation (32), Δrsisj can be expressed as follows:(59)Δrsisj=ΔRξmξnξsimξsjn.Similarly, evaluatingΔR using equation (31) and substituting it into equation (59), we obtain the following:(60)Δrsisj=Δrξmξn+pmnkΔrξkξsimξsjn.Equating the vector components on both sides of equation (60), we can write(61)Δxsisjl=Δxξmξnl+pmnkΔxξklξsimξsjn.Similarly, using equation (55), the components of rsi can be expressed as follows:(62)xsil=xξmlξsim.Therefore, taking the dot product ofrsi and Δrsisj and substituting the expressions for Δxsisjl and xsil given by equations (61) and (62), rsi·Δrsisj can be expressed as(63)rsi·Δrsisj=xsilΔxsisjl=xsilΔxξmξnl+pmnkΔxξklξsimξsjn=xsilξsimξsjnΔxξmξnl+xsilξsimξsjnpmnkΔxξkl.From the expression for the metric tensor in equation (14), Δg¯ij can be expressed as follows:(64)Δg¯ij=Δrsi·rsj+Δrsj·rsi.Therefore, the change in the covariant metric tensor components for thesi coordinates can be expressed as follows:(65a)Δg¯12=Δrs·rt+rs·Δrt.(65b)Δg¯11=2rs·Δrs.(65c)Δg¯22=2rt·Δrt.If we then expand the dot products in equation (64) and rearrange, Δg¯ij can be expressed as follows:(66)Δg¯ij=Δxsilxsjl+xsilΔxsjl=Δxξmlξsimxsjl+xsilΔxξmlξsjm=ξsimxsjl+ξsjmxsilΔxξml.From the chain rule applied to the differentiation of quotients, we can write the following:(67a)Δg¯12g¯22=Δg¯12g¯22−g¯12g¯22Δg¯22g¯22.(67b)Δg¯12g¯11=Δg¯12g¯11−g¯12g¯11Δg¯11g¯11.Substituting equation (66) into equations (67a)-(67b), we then have the following(68a)Δg¯12g¯22=1g¯22ξsmxtl+ξtmxslΔxξml−2g¯12g¯222ξtmxtlΔxξml.(68b)Δg¯12g¯11=1g¯11ξsmxtl+ξtmxslΔxξml−2g¯12g¯112ξsmxslΔxξml.Now, combining equations (58), (63), (67a) and (68a) into equation (53) and rearranging, we get the following expression for ΔgC¯1:(69)ΔgC¯1=−J2^ξsmxttlΔxξml+xslξtmξtnΔxξmξnl+xslξtmξtnpmnkΔxξkl−rt·rttg¯22ξsmxtl+ξtmxsl−2g¯12g¯22ξtmxtlΔxξml−g¯12g¯22ξtmxttlΔxξml−g¯12g¯22xtlξtmξtnΔxξmξnl+xtlξtmξtnpmnkΔxξkl.After rearranging and collecting terms, we can finally expressΔgC¯1 as follows:(70)ΔgC¯1=Almn1Δxξmξnl+Blk1Δxξkl,where(71)Almn1=−J2^ξtmξtnxsl−g¯12g¯22xtl,and(72)Blk1=−J2^xttlξsk−g¯12g¯22ξtk+ξtmξtnpmnkxsl−g¯12g¯22xtl−rt·rttg¯22ξtkxsl−2g¯12g¯22xtl+ξskxtl.Similarly, if we combine equations (58), (63), (67b) and (68b) into equation (54) and rearrange, we get the following expression for ΔgC¯2:(73)ΔgC¯2=−J2^ξtkxsslΔxξkl+xtlξsmξsnΔxξmξnl+xtlξsmξsnpmnkΔxξkl−rs·rssg¯11ξskxtl+ξtkxsl−2g¯12g¯11ξskxslΔxξkl−g¯12g¯11ξskxsslΔxξkl−g¯12g¯11xslξsmξsnΔxξmξnl+xslξsmξsnpmnkΔxξkl.Rearranging and collecting terms in (73), we can then write ΔgC¯2 as follows:(74)ΔgC¯2=Almn2Δxξmξnl+Blk2Δxξkl,where(75)Almn2=−J2^ξsmξsnxtl−g¯12g¯11xsl,and(76)Blk2=−J2^xsslξtk−g¯12g¯11ξsk+ξsmξsnpmnkxtl−g¯12g¯11xsl−rs·rssg¯11ξskxtl−2g¯12g¯11xsl+ξtkxsl.Using expressions forΔgC¯1 and ΔgC¯2 from equations (70) and (74), respectively, ΔgRp in equation (52) can now be expressed as(77)ΔgRp=ΔgC¯1ξsp+ΔgC¯2ξtp=Almn1ξsp+Almn2ξtpΔxξmξnl+Blk1ξsp+Blk2ξtpΔxξkl.We now have all the terms necessary to evaluategRp+ in equation (51).In Section2.4.1, the derivation of the linearized grid generation equations ((49) and (50)) excluded the curvature terms gRp. 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):(78)gR1rξ+R2rη+R1Δrξ+R2Δrη+ΔgR1rξ+ΔgR2rη.Separating equation (78) into its x and y components then yields the following expressions:(79a)gR1xξ+R2xη+R1Δxξ+R2Δxη+ΔgR1xξ+ΔgR2xη,and(79b)gR1yξ+R2yη+R1Δyξ+R2Δyη+ΔgR1yξ+ΔgR2yη.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 ΔgC¯1 and ΔgC¯2 given by equation (77), we obtain the following complete linearized grid generation equations with curvature control functions:(80)g22Δxξξ−2g12Δxξη+g11Δxξξ+P1+2Xηηxξ−Xξηxη+gR1Δxξ+P2+2Xξξxη−Xξηxξ+gR2Δxη+2Xηηyξ−XξηyηΔyξ+2Xξξyη−XξηyξΔyη+Almn1ξs+Almn2ξtxξΔxξmξnl+Blk1ξs+Blk2ξtxξΔxξkl+Almn1ηs+Almn2ηtxηΔxξmξnl+Blk1ηs+Blk2ηtxηΔxξkl=−Gx−gR1xξ−gR2xη.and(81)g22Δyξξ−2g12Δyξη+g11Δyξξ+2Yηηxξ−YξηxηΔxξ+2Yξξxη−YξηxξΔxη+P1+2Yηηyξ−Yξηyη+gR1Δyξ+P2+2Yξξyη−Yξηyξ+gR2Δyη+Almn1ξs+Almn2ξtyξΔxξmξnl+Blk1ξs+Blk2ξtyξΔxξkl+Almn1ηs+Almn2ηtyηΔxξmξnl+Blk1ηs+Blk2ηtyηΔxξkl=−Gy−gR1yξ−gR2yη.Equations (80) and (81) are the fully linearized equations that must be solved numerically at each grid point to compute the x and y 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 Δxξml and Δxξmξnl, leading to equations of the following form:(82)Cxij1Δxξ+Cxij2Δxη+Cxij3Δyξ+Cxij4Δyη+Cxij5Δxξξ+Cxij6Δxξη+Cxij7Δxηη+Cxij8Δyξξ+Cxij9Δyξη+Cxij10Δyηη=Rhxij.(83)Cyij1Δxξ+Cyij2Δxη+Cyij3Δyξ+Cyij4Δyη+Cyij5Δxξξ+Cyij6Δxξη+Cyij7Δxηη+Cyij8Δyξξ+Cyij9Δyξη+Cyij10Δyηη=Rhyij.where(84a)Cxij1=P1+2Xηηxξ−Xξηxη+gR1+B111xs+B112xt.(84b)Cxij2=P2+2Xξξxη−Xξηxξ+gR2+B121xs+B122xt.(84c)Cxij3=2Xηηyξ−Xξηyη+B211xs+B212xt.(84d)Cxij4=2Xξξyη−Xξηyξ+B221xs+B222xt.(84e)Cxij5=g22+A1111xs+A1112xt.(84f)Cxij6=−2g12+A1121xs+A1122xt.(84g)Cxij7=g11+A1221xs+A1222xt.(84h)Cxij8=A2111xs+A2112xt.(84i)Cxij9=A2121xs+A2122xt.(84j)Cxij10=A2221xs+A2222xt,and(85a)Cyij1=2Yηηxξ−Yξηxη+B111ys+B112yt.(85b)Cyij2=2Yξξxη−Yξηxξ+B121ys+B122yt.(85c)Cyij3=P1+2Yηηyξ−Yξηyη+gR1+B211ys+B212yt.(85d)Cyij4=P2+2Yξξyη−Yξηyξ+gR2+B221ys+B222yt.(85e)Cyij5=A1111ys+A1112yt.(85f)Cyij6=A1121ys+A1122yt.(85g)Cyij7=A1221ys+A1222yt.(85h)Cyij8=g22+A2111ys+A2112yt.(85i)Cyij9=−2g12+A2121ys+A2122yt.(85j)Cyij10=g11+A2221ys+A2222yt,and where(86a)Rhxij=−Gx−gR1xξ−gR2xη.(86b)Rhyij=−Gy−gR1yξ−gR2yη.In equations (84a)–(84j) to (86a)–(86b), the curvature control functions can be deactivated by setting all the Rp,Almnp and Blkp 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 x and y coordinates from iteration n to n+1, we have(87a)Δxξ=Δxi+1,j−Δxi−1,j2Δξ.(87b)Δyξ=Δyi+1,j−Δyi−1,j2Δξ.(87c)Δxη=Δxi,j+1−Δxi,j−12Δη.(87d)Δyη=Δyi,j+1−Δyi,j−12Δη.(87e)Δxξξ=Δxi+1,j−2Δxi,j+Δxi−1,jΔξ2.(87f)Δyξξ=Δyi+1,j−2Δyi,j+Δyi−1,jΔξ2.(87g)Δxηη=Δxi,j+1−2Δxi,j+Δxi,j−1Δη2.(87h)Δyηη=Δyi,j+1−2Δyi,j+Δyi,j−1Δη2.(87i)Δxξη=Δxi+1,j+1−Δxi+1,j−1−Δxi−1,j+1+Δxi−1,j−14ΔξΔη.(87j)Δyξη=Δyi+1,j+1−Δyi+1,j−1−Δyi−1,j+1+Δyi−1,j−14ΔξΔη.In the mesh generation process, the curvilinear coordinates,ξ and η, are the independent variables, and for convenience, they are set equal to the mesh indices i and j, respectively. This implies that ξ and η vary from 1 to imax and jmax, respectively, and thus Δξ=Δη=1 everywhere in the mesh. Therefore, substituting the finite-difference expressions from equations (87a)–(87j) into equations (82) and (83), yields the following:(88)−14Cxij6Δxi+1,j−1−12Cxij2−Cxij7Δxi,j−1+14Cxij6Δxi−1,j−1+12Cxij1+Cxij5Δxi+1,j−2Cxij5+Cxij7Δxi,j−12Cxij1−Cxij5Δxi−1,j+14Cxij6Δxi+1,j+1+12Cxij2+Cxij7Δxi,j+1−14Cxij6Δxi−1,j+1−14Cxij9Δyi+1,j−1−12Cxij4−Cxij10Δyi,j−1+14Cxij9Δyi−1,j−1+12Cxij3+Cxij8Δyi+1,j−2Cxij8+Cxij10Δyi,j−12Cxij3−Cxij8Δyi−1,j+14Cxij9Δyi+1,j+1+12Cxij4+Cxij10Δyi,j+1−14Cxij9Δyi−1,j+1=Rhxij,and(89)−14Cyij6Δxi+1,j−1−12Cyij2−Cyij7Δxi,j−1+14Cyij6Δxi−1,j−1+12Cyij1+Cyij5Δxi+1,j−2Cyij5+Cyij7Δxi,j−12Cyij1−Cyij5Δxi−1,j+14Cyij6Δxi+1,j+1+12Cyij2+Cyij7Δxi,j+1−14Cyij6Δxi−1,j+1−14Cyij9Δyi+1,j−1−12Cyij4−Cyij10Δyi,j−1+14Cyij9Δyi−1,j−1+12Cyij3+Cyij8Δyi+1,j−2Cyij8+Cyij10Δyi,j−12Cyij3−Cyij8Δyi−1,j+14Cyij9Δyi+1,j+1+12Cyij4+Cyij10Δyi,j+1−14Cyij9Δyi−1,j+1=Rhyij.For the ease of presentation, we introduce the parametersJxijk and Jyijk to represent the coefficients in equations (88) and (89) as follows:(90)Jxij1=14Cxij6;Jxij2=−12Cxij2−Cxij7;Jxij3=12Cxij1+Cxij5Jxij4=−2Cxij5+Cxij7;Jxij5=−12Cxij1−Cxij5;Jxij6=12Cxij2+Cxij7Jxij7=14Cxij9;Jxij8=−12Cxij4−Cxij10;Jxij9=12Cxij3+Cxij8Jxij10=−2Cxij8+Cxij10;Jxij11=−12Cxij3−Cxij8;Jxij12=12Cxij4+Cxij10,and(91)Jyij1=14Cyij6;Jyij2=−12Cyij2−Cyij7;Jyij3=12Cyij1+Cyij5Jyij4=−2Cyij5+Cyij7;Jyij5=−12Cyij1−Cyij5;Jyij6=12Cyij2+Cyij7Jyij7=14Cyij9;Jyij8=−12Cyij4−Cyij10;Jyij9=12Cyij3+Cyij8Jyij10=−2Cyij8+Cyij10;Jyij11=−12Cyij3−Cyij8;Jyij12=12Cyij4+Cyij10.Thus, inserting equations (90) and (91) into equations (88) and (89), respectively, we obtain the following(92)Jxij5Δxi−1,j−Jxij1Δxi+1,j−1+Jxij2Δxi,j−1+Jxij1Δxi−1,j−1+Jxij3Δxi+1,j+Jxij4Δxi,j+Jxij1Δxi+1,j+1+Jxij6Δxi,j+1−Jxij1Δxi−1,j+1−Jxij7Δyi+1,j−1+Jxij8Δyi,j−1+Jxij7Δyi−1,j−1+Jxij9Δyi+1,j+Jxij10Δyi,j+Jxij11Δyi−1,j+Jxij7Δyi+1,j+1+Jxij12Δyi,j+1−Jxij7Δyi−1,j+1=Rhxij,and(93)Jyij5Δxi−1,j−Jyij1Δxi+1,j−1+Jyij2Δxi,j−1+Jyij1Δxi−1,j−1+Jyij3Δxi+1,j+Jyij4Δxi,j+Jyij1Δxi+1,j+1+Jyij6Δxi,j+1−Jyij1Δxi−1,j+1−Jyij7Δyi+1,j−1+Jyij8Δyi,j−1+Jyij7Δyi−1,j−1+Jyij9Δyi+1,j+Jyij10Δyi,j+Jyij11Δyi−1,j+Jyij7Δyi+1,j+1+Jyij12Δyi,j+1−Jyij7Δyi−1,j+1=Rhyij.Equations (92) and (93) are the complete discretized forms of the grid generation equations that must be solved at each grid point of the mesh. As can be seen from these equations, the second-order finite difference expressions result in a 9-point stencil surrounding each grid point (i, j), as illustrated in Figure 3 for a 4×4 grid.Figure 3
4×4 grid with a 9-point stencil.Equations (92) and (93) yield a linear system of equations involving a nonsymmetric Jacobian matrix, with entries given by the Jxijk and Jyijk terms. This matrix is populated by sweeping the mesh one row at a time along the i direction, and placing the coefficients multiplying the x and y coordinates of every point in adjacent rows. This ordering of the mesh nodes yields the matrix with the smallest bandwidth of nonzero entries about the diagonal, as illustrated by the matrix sparsity pattern shown in Figure 4. As shown, the matrix consists of three diagonal bands, each consisting of a block tridiagonal band of 2 × 2 blocks. The distance of the upper and lower bands from the diagonal is only 2(jmax - 5) columns, where jmax is the upper limit of the smaller of the two mesh indices, i,j. As shown in Figure 4(a), the total bandwidth of the nonzero elements is very small relative to the overall size of the matrix.Figure 4
Matrix sparsity pattern. (a) Sparsity pattern. (b) Details of upper-left corner.
(a)(b)To solve this system, a number of solvers and preconditioners, drawn from the Portable Extensible Toolkit for Scientific Computation (PETSc) library [28], are employed. The solution vector obtained from the PETSc solution gives the changes in the coordinates, Δxi,j and Δyi,j, at each grid point. The updated values of the coordinates at iteration n+1 are then obtained as follows:(94)xi,jn+1=xi,jn+Δxi,jn+1,yi,jn+1=yi,jn+Δyi,jn+1.


### 3.2. Evaluation of Solvers and Preconditioners

The PETSc library is a collection of data structures and routines for the solution of scientific and engineering problems formulated by partial differential equations. It consists of a number of libraries which contain modules for creating both sequential and parallel vectors and matrices as well as modules dealing with linear and nonlinear solvers. It is structured such that each module can be invoked by a set of calling sequences facilitating the solution of PDE’s [29]. Once the interface to the PETSc library is established, all solvers and preconditioners are accessible. This approach enables the evaluation of various solvers and preconditioners for a particular problem to determine the best combination for tackling the problem at hand.


#### 3.2.1. Selected Solvers and Preconditioners

In the applications presented Section4, virtually all methods available in the PETSc library were evaluated, and the set of preconditioners and solvers that turned out to be the most effective were identified as follows:Preconditioners:ILU: Incomplete lower-upper (LU) factorization.SOR: Successive over-relaxation variant of Gauss-Seidel.EISENSTAT: SSOR (symmetric successive over-relaxation, symmetric Gauss-Seidel) incorporating Eisenstat’s “trick” to reduce computation.Solvers:LU: Direct solver based on lower-upper factorization.BiCGSTAB: Stabilized version of Biconjugate Gradient method.GMRES: Generalized Minimal Residual method, optionally with accelerated restart.PIPELINED GCR: Pipelined Generalized Conjugate Residual method.


### 3.3. Implementation in Fortran

The numerical solution of the grid generation equations with curvature functions was implemented in a fortran code. The simplified algorithm describing the implementation of the solution method is summarized in Algorithm1 below.Algorithm 1: Algorithm for the implementation in Fortran.
(1)
Read the initial algebraic mesh(2)
Initialize various PETSc variables(3)
Forj=2⟶jmax−1do(4)
Fori=2⟶imax−1do(5)
Set the exact no. of nonzeros in each row of the matrix depending on various conditions(6)
End(7)
End(8)
Set up the PETSc matrix and vector(9)
Foriter=1⟶niterdo(10)
Forj=2⟶jmax−1do(11)
Fori=2⟶imax−1do(12)
(1) Calculate derivatives ofx and y w.r.t. to ξ and η (equations (87a)–(87j))(13)
(2) Evaluate the metric tensor components (equations (14)–(16))(14)
(3) Calculate the control functions (stretching and orthogonality) (equations (19) and (26a)-(26b))(15)
(4) Evaluate the curvature control terms (equations (10) and (28a)-(28b))(16)
(5) Evaluate the RHS vector components forx and y (equations (86a)-(86b))(17)
(6) Evaluate the coefficients of the linearized equations (equations (92)-(93))(18)
(7) Populate the matrix entries in the nonzero locations(19)
End(20)
End(21)
(1) Assemble the PETSc vector and matrix(22)
(2) Solve the linear system using any of the KSP methods(23)
(3) Extract the solution vector(24)
(4) Update the values ofx and y(25)
End(26)
(1) Delete the PETSc objects(27)
(2) Write the output file


## 3.1. Setup of the Matrix Equation

Writing the finite-difference expressions in (34) in terms of the change in x and y coordinates from iteration n to n+1, we have(87a)Δxξ=Δxi+1,j−Δxi−1,j2Δξ.(87b)Δyξ=Δyi+1,j−Δyi−1,j2Δξ.(87c)Δxη=Δxi,j+1−Δxi,j−12Δη.(87d)Δyη=Δyi,j+1−Δyi,j−12Δη.(87e)Δxξξ=Δxi+1,j−2Δxi,j+Δxi−1,jΔξ2.(87f)Δyξξ=Δyi+1,j−2Δyi,j+Δyi−1,jΔξ2.(87g)Δxηη=Δxi,j+1−2Δxi,j+Δxi,j−1Δη2.(87h)Δyηη=Δyi,j+1−2Δyi,j+Δyi,j−1Δη2.(87i)Δxξη=Δxi+1,j+1−Δxi+1,j−1−Δxi−1,j+1+Δxi−1,j−14ΔξΔη.(87j)Δyξη=Δyi+1,j+1−Δyi+1,j−1−Δyi−1,j+1+Δyi−1,j−14ΔξΔη.In the mesh generation process, the curvilinear coordinates,ξ and η, are the independent variables, and for convenience, they are set equal to the mesh indices i and j, respectively. This implies that ξ and η vary from 1 to imax and jmax, respectively, and thus Δξ=Δη=1 everywhere in the mesh. Therefore, substituting the finite-difference expressions from equations (87a)–(87j) into equations (82) and (83), yields the following:(88)−14Cxij6Δxi+1,j−1−12Cxij2−Cxij7Δxi,j−1+14Cxij6Δxi−1,j−1+12Cxij1+Cxij5Δxi+1,j−2Cxij5+Cxij7Δxi,j−12Cxij1−Cxij5Δxi−1,j+14Cxij6Δxi+1,j+1+12Cxij2+Cxij7Δxi,j+1−14Cxij6Δxi−1,j+1−14Cxij9Δyi+1,j−1−12Cxij4−Cxij10Δyi,j−1+14Cxij9Δyi−1,j−1+12Cxij3+Cxij8Δyi+1,j−2Cxij8+Cxij10Δyi,j−12Cxij3−Cxij8Δyi−1,j+14Cxij9Δyi+1,j+1+12Cxij4+Cxij10Δyi,j+1−14Cxij9Δyi−1,j+1=Rhxij,and(89)−14Cyij6Δxi+1,j−1−12Cyij2−Cyij7Δxi,j−1+14Cyij6Δxi−1,j−1+12Cyij1+Cyij5Δxi+1,j−2Cyij5+Cyij7Δxi,j−12Cyij1−Cyij5Δxi−1,j+14Cyij6Δxi+1,j+1+12Cyij2+Cyij7Δxi,j+1−14Cyij6Δxi−1,j+1−14Cyij9Δyi+1,j−1−12Cyij4−Cyij10Δyi,j−1+14Cyij9Δyi−1,j−1+12Cyij3+Cyij8Δyi+1,j−2Cyij8+Cyij10Δyi,j−12Cyij3−Cyij8Δyi−1,j+14Cyij9Δyi+1,j+1+12Cyij4+Cyij10Δyi,j+1−14Cyij9Δyi−1,j+1=Rhyij.For the ease of presentation, we introduce the parametersJxijk and Jyijk to represent the coefficients in equations (88) and (89) as follows:(90)Jxij1=14Cxij6;Jxij2=−12Cxij2−Cxij7;Jxij3=12Cxij1+Cxij5Jxij4=−2Cxij5+Cxij7;Jxij5=−12Cxij1−Cxij5;Jxij6=12Cxij2+Cxij7Jxij7=14Cxij9;Jxij8=−12Cxij4−Cxij10;Jxij9=12Cxij3+Cxij8Jxij10=−2Cxij8+Cxij10;Jxij11=−12Cxij3−Cxij8;Jxij12=12Cxij4+Cxij10,and(91)Jyij1=14Cyij6;Jyij2=−12Cyij2−Cyij7;Jyij3=12Cyij1+Cyij5Jyij4=−2Cyij5+Cyij7;Jyij5=−12Cyij1−Cyij5;Jyij6=12Cyij2+Cyij7Jyij7=14Cyij9;Jyij8=−12Cyij4−Cyij10;Jyij9=12Cyij3+Cyij8Jyij10=−2Cyij8+Cyij10;Jyij11=−12Cyij3−Cyij8;Jyij12=12Cyij4+Cyij10.Thus, inserting equations (90) and (91) into equations (88) and (89), respectively, we obtain the following(92)Jxij5Δxi−1,j−Jxij1Δxi+1,j−1+Jxij2Δxi,j−1+Jxij1Δxi−1,j−1+Jxij3Δxi+1,j+Jxij4Δxi,j+Jxij1Δxi+1,j+1+Jxij6Δxi,j+1−Jxij1Δxi−1,j+1−Jxij7Δyi+1,j−1+Jxij8Δyi,j−1+Jxij7Δyi−1,j−1+Jxij9Δyi+1,j+Jxij10Δyi,j+Jxij11Δyi−1,j+Jxij7Δyi+1,j+1+Jxij12Δyi,j+1−Jxij7Δyi−1,j+1=Rhxij,and(93)Jyij5Δxi−1,j−Jyij1Δxi+1,j−1+Jyij2Δxi,j−1+Jyij1Δxi−1,j−1+Jyij3Δxi+1,j+Jyij4Δxi,j+Jyij1Δxi+1,j+1+Jyij6Δxi,j+1−Jyij1Δxi−1,j+1−Jyij7Δyi+1,j−1+Jyij8Δyi,j−1+Jyij7Δyi−1,j−1+Jyij9Δyi+1,j+Jyij10Δyi,j+Jyij11Δyi−1,j+Jyij7Δyi+1,j+1+Jyij12Δyi,j+1−Jyij7Δyi−1,j+1=Rhyij.Equations (92) and (93) are the complete discretized forms of the grid generation equations that must be solved at each grid point of the mesh. As can be seen from these equations, the second-order finite difference expressions result in a 9-point stencil surrounding each grid point (i, j), as illustrated in Figure 3 for a 4×4 grid.Figure 3
4×4 grid with a 9-point stencil.Equations (92) and (93) yield a linear system of equations involving a nonsymmetric Jacobian matrix, with entries given by the Jxijk and Jyijk terms. This matrix is populated by sweeping the mesh one row at a time along the i direction, and placing the coefficients multiplying the x and y coordinates of every point in adjacent rows. This ordering of the mesh nodes yields the matrix with the smallest bandwidth of nonzero entries about the diagonal, as illustrated by the matrix sparsity pattern shown in Figure 4. As shown, the matrix consists of three diagonal bands, each consisting of a block tridiagonal band of 2 × 2 blocks. The distance of the upper and lower bands from the diagonal is only 2(jmax - 5) columns, where jmax is the upper limit of the smaller of the two mesh indices, i,j. As shown in Figure 4(a), the total bandwidth of the nonzero elements is very small relative to the overall size of the matrix.Figure 4
Matrix sparsity pattern. (a) Sparsity pattern. (b) Details of upper-left corner.
(a)(b)To solve this system, a number of solvers and preconditioners, drawn from the Portable Extensible Toolkit for Scientific Computation (PETSc) library [28], are employed. The solution vector obtained from the PETSc solution gives the changes in the coordinates, Δxi,j and Δyi,j, at each grid point. The updated values of the coordinates at iteration n+1 are then obtained as follows:(94)xi,jn+1=xi,jn+Δxi,jn+1,yi,jn+1=yi,jn+Δyi,jn+1.


## 3.2. Evaluation of Solvers and Preconditioners

The PETSc library is a collection of data structures and routines for the solution of scientific and engineering problems formulated by partial differential equations. It consists of a number of libraries which contain modules for creating both sequential and parallel vectors and matrices as well as modules dealing with linear and nonlinear solvers. It is structured such that each module can be invoked by a set of calling sequences facilitating the solution of PDE’s [29]. Once the interface to the PETSc library is established, all solvers and preconditioners are accessible. This approach enables the evaluation of various solvers and preconditioners for a particular problem to determine the best combination for tackling the problem at hand.


### 3.2.1. Selected Solvers and Preconditioners

In the applications presented Section4, virtually all methods available in the PETSc library were evaluated, and the set of preconditioners and solvers that turned out to be the most effective were identified as follows:Preconditioners:ILU: Incomplete lower-upper (LU) factorization.SOR: Successive over-relaxation variant of Gauss-Seidel.EISENSTAT: SSOR (symmetric successive over-relaxation, symmetric Gauss-Seidel) incorporating Eisenstat’s “trick” to reduce computation.Solvers:LU: Direct solver based on lower-upper factorization.BiCGSTAB: Stabilized version of Biconjugate Gradient method.GMRES: Generalized Minimal Residual method, optionally with accelerated restart.PIPELINED GCR: Pipelined Generalized Conjugate Residual method.


## 3.2.1. Selected Solvers and Preconditioners

In the applications presented Section4, virtually all methods available in the PETSc library were evaluated, and the set of preconditioners and solvers that turned out to be the most effective were identified as follows:Preconditioners:ILU: Incomplete lower-upper (LU) factorization.SOR: Successive over-relaxation variant of Gauss-Seidel.EISENSTAT: SSOR (symmetric successive over-relaxation, symmetric Gauss-Seidel) incorporating Eisenstat’s “trick” to reduce computation.Solvers:LU: Direct solver based on lower-upper factorization.BiCGSTAB: Stabilized version of Biconjugate Gradient method.GMRES: Generalized Minimal Residual method, optionally with accelerated restart.PIPELINED GCR: Pipelined Generalized Conjugate Residual method.


## 3.3. Implementation in Fortran

The numerical solution of the grid generation equations with curvature functions was implemented in a fortran code. The simplified algorithm describing the implementation of the solution method is summarized in Algorithm1 below.Algorithm 1: Algorithm for the implementation in Fortran.
(1)
Read the initial algebraic mesh(2)
Initialize various PETSc variables(3)
Forj=2⟶jmax−1do(4)
Fori=2⟶imax−1do(5)
Set the exact no. of nonzeros in each row of the matrix depending on various conditions(6)
End(7)
End(8)
Set up the PETSc matrix and vector(9)
Foriter=1⟶niterdo(10)
Forj=2⟶jmax−1do(11)
Fori=2⟶imax−1do(12)
(1) Calculate derivatives ofx and y w.r.t. to ξ and η (equations (87a)–(87j))(13)
(2) Evaluate the metric tensor components (equations (14)–(16))(14)
(3) Calculate the control functions (stretching and orthogonality) (equations (19) and (26a)-(26b))(15)
(4) Evaluate the curvature control terms (equations (10) and (28a)-(28b))(16)
(5) Evaluate the RHS vector components forx and y (equations (86a)-(86b))(17)
(6) Evaluate the coefficients of the linearized equations (equations (92)-(93))(18)
(7) Populate the matrix entries in the nonzero locations(19)
End(20)
End(21)
(1) Assemble the PETSc vector and matrix(22)
(2) Solve the linear system using any of the KSP methods(23)
(3) Extract the solution vector(24)
(4) Update the values ofx and y(25)
End(26)
(1) Delete the PETSc objects(27)
(2) Write the output file


## 4. Applications

The solution process developed herein was applied to the generation of two C-meshes around a two-dimensional airfoil. The first mesh is suitable for the solution of the Euler equations, while the second is a more refined mesh suitable for the Navier–Stokes equations, in which the height of the first layer of cells on the airfoil boundary is set to10−6×c, where c is the airfoil chord. The Euler mesh consists of 321 × 49 (15,729) points, yielding a matrix with dimensions 29986×29986, while the Navier–Stokes mesh contains 321 × 97 (31,137) points, resulting in a matrix with dimensions 60610×60610. These matrices are very sparse, with the number of nonzero entries being equal to 530,980 for the Euler mesh, and 1,081,060 for the Navier-Stokes mesh, such that only approximately 0.06% and 0.03% of the matrix entries are nonzero, respectively. Figure 5 shows the initial algebraic mesh and domain boundaries. The far field regions of the Euler and Navier–Stokes meshes are virtually identical, as the meshes only differ in the region close to the airfoil surface. Figures 6(a) and 6(b) show the two meshes in the vicinity of the airfoil boundary. The initial meshes are generated using transfinite interpolation from the boundaries.Figure 5
Far-field of initial algebraic Euler and Navier–Stokes meshes.Figure 6
Closeup of initial algebraic meshes. (a) Euler mesh. (b) Navier–Stokes mesh.
(a)(b)Both the Euler and Navier–Stokes meshes were generated with and without curvature control functions, in order to compare the properties of the meshes as well as the convergence characteristics of the solution process. In addition, for the meshes generated without curvature functions, optimizations were done with and without orthogonality functions. The reason for this is that meshes generated with curvature functions are inherently more orthogonal than meshes generated without them, and for many applications, including the ones presented here, orthogonality functions are not required. Therefore, in the following figures the effects of turning on and off both the curvature functions and the orthogonality functions are presented and discussed.The first optimizations were performed on the Euler meshes and the results are shown in Figures7 to 10. Figure 7 shows the far field of the Euler meshes optimized without curvature functions, where 7(a) and 7(b) show the meshes optimized without and with orthogonality functions, respectively. Closeup views of these meshes are provided in Figure 8. It is clear from these figures that orthogonality functions are required to obtain a high-quality mesh in this case, as is virtually always the case when using a Laplace-based operator. To highlight the impact of the curvature functions on orthogonality, Figures 9 and 10 show the Euler mesh optimized with curvature functions, but without orthogonality functions. These figures can be compared to Figures 7(a) and 8(a), respectively, to see the impact of the curvature functions on the mesh. As can be seen in these figures, the mesh optimized with curvature functions is fully orthogonal on the airfoil boundary, without the need for dedicated orthogonality functions, in contrast to the mesh optimized without curvature functions. Also, evident in Figure 9 is that the mesh optimized with curvature functions retains the global curvature of the C-mesh topology, as well as the mesh clustering prescribed via the boundary point distributions. In contrast, the mesh optimized without curvature functions concentrates the mesh points near the airfoil leading edge. This tendency of Laplace-based operators can also be seen when comparing Figures 8(a) to 10. In the latter figure, it can be seen that for the mesh optimized with curvature functions the height of the first layer of cells around the airfoil is consistent with the point distribution prescribed on the boundaries (as shown in Figure 6(a)), whereas for the mesh optimized without curvature functions this spacing is significantly altered.Figure 7
Far-field of Euler meshes optimized without curvature functions. (a) Optimized without orthogonality functions. (b) Optimized with orthogonality functions.
(a)(b)Figure 8
Closeup of Euler meshes optimized without curvature functions. (a) Optimized without orthogonality functions. (b) Optimized with orthogonality functions.
(a)(b)Figure 9
Far-field of Euler mesh optimized with curvature control functions but without orthogonality functions, showing the impact of curvature functions when compared to Figure7(a).Figure 10
Closeup of Euler mesh optimized with curvature control functions but without orthogonality functions, showing the impact of curvature functions when compared to Figure8(a).Similar conclusions can be drawn for the optimizations of the Navier–Stokes mesh, as shown in Figures11 to 14. Figure 11 shows the far field of the Navier-–Stokes meshes optimized without curvature functions, with and without orthogonality functions. Closeup views of these meshes are provided in Figure 12. Figures 13 and 14 show the Navier–Stokes mesh optimized with curvature functions, but without orthogonality functions. These figures can be compared to Figures 11(a) and 12(a), respectively, to see the impact of the curvature functions on the mesh. Again, the mesh optimized with curvature functions is fully orthogonal on the airfoil boundary and retains the global curvature of the C-mesh in both the near field and far field regions.Figure 11
Far-field of the Navier–Stokes meshes optimized without curvature functions. (a) Optimized without orthogonality functions. (b) Optimized with orthogonality functions.
(a)(b)Figure 12
Closeup of the Navier–Stokes meshes optimized without curvature functions. (a) Optimized without orthogonality functions. (b) Optimized with orthogonality functions.
(a)(b)Figure 13
Far-field of the Navier–Stokes mesh optimized with curvature control functions but without orthogonality functions, showing the impact of curvature functions when compared to Figure11(a).Figure 14
Closeup of the Navier–Stokes mesh optimized with curvature control functions but without orthogonality functions, showing the impact of curvature functions when compared to Figure12(a).


### 4.1. Convergence of the Solution Process

The purpose of the sample solutions presented herein is to validate the Newton linearization of the equations and the setup of the system matrix. The results presented do not constitute an exhaustive evaluation of the various solvers and preconditioners that can be used to solve the system, but rather provide a range of results illustrating the efficiency of the solution process for two-dimensional applications. Most available combinations of preconditioners and solvers in the PETSc library were tested for the cases presented here; however, only the more effective combinations are included in the figures. In all cases, all internal settings of the various PETSc preconditioners and solvers were set to their default values, except for the relative tolerance for the stopping criterion for solver subiterations, which was set to10−4. For more information on the detailed PETSc settings for the various solvers, the reader if referred to the PETSc manual [29].The convergence curves for the generation of the Euler mesh with curvature control functions are shown in Figures15 and 16. Figure 15 shows the convergence as a function of iteration, while Figure 16 shows the convergence as a function of time. As shown in Figure 15, no reduction in the maximum residual was achieved in the first 5 iterations, due to the fact that the curvature functions were activated progressively in these early iterations, reaching their full value at iteration 5. In all cases presented here, a gradual activation of the curvature functions was required to smooth out any spurious curvatures in the initial algebraic mesh and ensure convergence. As shown in the convergence plots, once the curvature functions are fully activated, the solution converges rapidly.Figure 15
Convergence for the Euler mesh vs. iteration, with curvature functions.Figure 16
Convergence for the Euler mesh vs. time, with curvature functions.As shown in Figure16, the most efficient solver for the generation of the Euler mesh with curvature functions is the direct “LU” solver. This solver achieved machine accuracy in only 18 iterations, with an elapsed time of 3 seconds, running on a single Intel i7-7700K, 4.2 GHz core.The convergence curves for the generation of the Navier–Stokes mesh with curvature functions are shown in Figures17 and 18. Figure 17 shows the convergence as a function of iteration, and Figure 18 shows the convergence as a function of time. Here again, the curvature functions were activated progressively in the first 5 iterations. As shown in Figure 18, the most efficient solver in this case is again the “LU” solver, achieving machine accuracy in 26 iterations, with an elapsed time of 16 seconds. Also interesting in this figure is the performance of the “EISENSTAT” preconditioner with the “PIPELINED GCR” solver; although this combination did not achieve machine accuracy, it achieved what can be considered to be sufficient “engineering” accuracy, and did so in significantly less elapsed time than the other preconditioner-solver combinations.Figure 17
Convergence of the Navier–Stokes mesh vs. iteration, with curvature functions.Figure 18
Convergence of the Navier–Stokes mesh vs. time, with curvature functions.Based on these results, the most effective overall solver for the solution of the grid generation equations with curvature control functions is the “LU” solver, as it performed well for both the Euler and Navier–Stokes meshes. Hence, using this solver, the meshes were also generated without curvature control functions to illustrate the impact of the latter on convergence rates. Figures19 and 20 show the comparison of the convergence curves for both the Euler and Navier–Stokes meshes with and without curvature functions. Figure 19 shows the convergence as a function of iteration, and Figure 20 shows the convergence as a function of time. As shown, the inclusion of curvature functions reduces the convergence rate and increases the elapsed time to achieve machine accuracy. However, for practical applications, where a reduction of 10−5 in the residual is typically sufficient, the addition of curvature functions does not significantly increase the elapsed time.Figure 19
Convergence rates vs. iteration, with and without curvature functions.Figure 20
Convergence rates vs. time, with and without curvature functions.It is perhaps not surprising that the “LU” solver outperforms the other solvers in most cases presented here, given that it is a direct solver. However, it is important to note that this will not necessarily be the case for large three-dimensional applications, for which the cost of a direct solver may be too high. What is also observed in the applications presented herein is that the curvature functions degrade the conditioning of the system matrix, as evidenced by an increase in the condition number. For both the Euler and Navier–Stokes meshes generated without curvature functions, the condition number of the matrices does not exceedO104 throughout the iterative solution process. In contrast, for Euler meshes generated with curvature control functions, the condition number of the matrices increases to O105, and for Navier–Stokes meshes the conditions numbers reach O106 in some cases. This indicates that the very high cell aspect ratios in the boundary layer of the Navier–Stokes meshes have a more detrimental impact on the conditioning of the matrix when curvature functions are included. This is likely the reason that the “SOR + BiCGSTAB”, “SOR + GMRES”, and “EISENSTAT + PIPELINED GCR” combinations did not achieve machine accuracy for the Navier–Stokes meshes with curvature functions, as shown in Figures 17, 18 and 20.Figure21 shows the order of convergence achieved for the Euler and Navier–Stokes meshes, with and without curvatures functions. The order of convergence, p, is calculated as follows:(95)p=limn⟶∞logRmaxn+1/RmaxnlogRmaxn/Rmaxn−1,where Rmaxn is the maximum residual anywhere in the mesh at iteration n.Figure 21
Order of convergence for the Euler and Navier–Stokes meshes, with and without curvature functions.As shown, the numerical solutions without curvature functions achieve quadraticp=2 convergence for both the Euler and Navier–Stokes meshes. This indicates that the very high cell aspect-ratios of O105 in the boundary layer of the Navier–Stokes mesh do not significantly impact the convergence of the solution. This can also be observed in Figure 19, where the converge rates for the Euler and Navier–Stokes meshes generated without curvature functions are virtually identical. In contrast, the solutions for the meshes generated with curvature functions achieve orders of convergence of only 1.3 to 1.8, with the Euler mesh peaking at 1.77 and the Navier–Stokes mesh peaking at 1.36.It should be noted that the convergence curves shown in Figures15 to 21 pertain to meshes generated without orthogonality functions. The reason for this is that at the time of this writing, the implementation of orthogonality functions was not fully integrated with the newly implemented solution process, and thus the convergence rates for optimizations with orthogonality functions could not be compared to those obtained without orthogonality functions. Furthermore, since the optimizations with curvature functions did not require orthogonality functions, the latter were not included in the comparisons of convergence rates in order to isolate the impact of the curvature functions on convergence. Thus, the differences in convergence rates shown in Figures 19 and 20 are due solely to the addition of curvature functions. The proper integration of the orthogonality functions with the newly implemented PETSc-based solution process will be part of future work.


## 4.1. Convergence of the Solution Process

The purpose of the sample solutions presented herein is to validate the Newton linearization of the equations and the setup of the system matrix. The results presented do not constitute an exhaustive evaluation of the various solvers and preconditioners that can be used to solve the system, but rather provide a range of results illustrating the efficiency of the solution process for two-dimensional applications. Most available combinations of preconditioners and solvers in the PETSc library were tested for the cases presented here; however, only the more effective combinations are included in the figures. In all cases, all internal settings of the various PETSc preconditioners and solvers were set to their default values, except for the relative tolerance for the stopping criterion for solver subiterations, which was set to10−4. For more information on the detailed PETSc settings for the various solvers, the reader if referred to the PETSc manual [29].The convergence curves for the generation of the Euler mesh with curvature control functions are shown in Figures15 and 16. Figure 15 shows the convergence as a function of iteration, while Figure 16 shows the convergence as a function of time. As shown in Figure 15, no reduction in the maximum residual was achieved in the first 5 iterations, due to the fact that the curvature functions were activated progressively in these early iterations, reaching their full value at iteration 5. In all cases presented here, a gradual activation of the curvature functions was required to smooth out any spurious curvatures in the initial algebraic mesh and ensure convergence. As shown in the convergence plots, once the curvature functions are fully activated, the solution converges rapidly.Figure 15
Convergence for the Euler mesh vs. iteration, with curvature functions.Figure 16
Convergence for the Euler mesh vs. time, with curvature functions.As shown in Figure16, the most efficient solver for the generation of the Euler mesh with curvature functions is the direct “LU” solver. This solver achieved machine accuracy in only 18 iterations, with an elapsed time of 3 seconds, running on a single Intel i7-7700K, 4.2 GHz core.The convergence curves for the generation of the Navier–Stokes mesh with curvature functions are shown in Figures17 and 18. Figure 17 shows the convergence as a function of iteration, and Figure 18 shows the convergence as a function of time. Here again, the curvature functions were activated progressively in the first 5 iterations. As shown in Figure 18, the most efficient solver in this case is again the “LU” solver, achieving machine accuracy in 26 iterations, with an elapsed time of 16 seconds. Also interesting in this figure is the performance of the “EISENSTAT” preconditioner with the “PIPELINED GCR” solver; although this combination did not achieve machine accuracy, it achieved what can be considered to be sufficient “engineering” accuracy, and did so in significantly less elapsed time than the other preconditioner-solver combinations.Figure 17
Convergence of the Navier–Stokes mesh vs. iteration, with curvature functions.Figure 18
Convergence of the Navier–Stokes mesh vs. time, with curvature functions.Based on these results, the most effective overall solver for the solution of the grid generation equations with curvature control functions is the “LU” solver, as it performed well for both the Euler and Navier–Stokes meshes. Hence, using this solver, the meshes were also generated without curvature control functions to illustrate the impact of the latter on convergence rates. Figures19 and 20 show the comparison of the convergence curves for both the Euler and Navier–Stokes meshes with and without curvature functions. Figure 19 shows the convergence as a function of iteration, and Figure 20 shows the convergence as a function of time. As shown, the inclusion of curvature functions reduces the convergence rate and increases the elapsed time to achieve machine accuracy. However, for practical applications, where a reduction of 10−5 in the residual is typically sufficient, the addition of curvature functions does not significantly increase the elapsed time.Figure 19
Convergence rates vs. iteration, with and without curvature functions.Figure 20
Convergence rates vs. time, with and without curvature functions.It is perhaps not surprising that the “LU” solver outperforms the other solvers in most cases presented here, given that it is a direct solver. However, it is important to note that this will not necessarily be the case for large three-dimensional applications, for which the cost of a direct solver may be too high. What is also observed in the applications presented herein is that the curvature functions degrade the conditioning of the system matrix, as evidenced by an increase in the condition number. For both the Euler and Navier–Stokes meshes generated without curvature functions, the condition number of the matrices does not exceedO104 throughout the iterative solution process. In contrast, for Euler meshes generated with curvature control functions, the condition number of the matrices increases to O105, and for Navier–Stokes meshes the conditions numbers reach O106 in some cases. This indicates that the very high cell aspect ratios in the boundary layer of the Navier–Stokes meshes have a more detrimental impact on the conditioning of the matrix when curvature functions are included. This is likely the reason that the “SOR + BiCGSTAB”, “SOR + GMRES”, and “EISENSTAT + PIPELINED GCR” combinations did not achieve machine accuracy for the Navier–Stokes meshes with curvature functions, as shown in Figures 17, 18 and 20.Figure21 shows the order of convergence achieved for the Euler and Navier–Stokes meshes, with and without curvatures functions. The order of convergence, p, is calculated as follows:(95)p=limn⟶∞logRmaxn+1/RmaxnlogRmaxn/Rmaxn−1,where Rmaxn is the maximum residual anywhere in the mesh at iteration n.Figure 21
Order of convergence for the Euler and Navier–Stokes meshes, with and without curvature functions.As shown, the numerical solutions without curvature functions achieve quadraticp=2 convergence for both the Euler and Navier–Stokes meshes. This indicates that the very high cell aspect-ratios of O105 in the boundary layer of the Navier–Stokes mesh do not significantly impact the convergence of the solution. This can also be observed in Figure 19, where the converge rates for the Euler and Navier–Stokes meshes generated without curvature functions are virtually identical. In contrast, the solutions for the meshes generated with curvature functions achieve orders of convergence of only 1.3 to 1.8, with the Euler mesh peaking at 1.77 and the Navier–Stokes mesh peaking at 1.36.It should be noted that the convergence curves shown in Figures15 to 21 pertain to meshes generated without orthogonality functions. The reason for this is that at the time of this writing, the implementation of orthogonality functions was not fully integrated with the newly implemented solution process, and thus the convergence rates for optimizations with orthogonality functions could not be compared to those obtained without orthogonality functions. Furthermore, since the optimizations with curvature functions did not require orthogonality functions, the latter were not included in the comparisons of convergence rates in order to isolate the impact of the curvature functions on convergence. Thus, the differences in convergence rates shown in Figures 19 and 20 are due solely to the addition of curvature functions. The proper integration of the orthogonality functions with the newly implemented PETSc-based solution process will be part of future work.


## 5. Conclusions and Recommendations

In the work presented herein, the grid generation equations with and without curvature control functions were discretized using second-order finite differences and linearized using Newton’s method. To solve the resulting linear system, a number of solvers and preconditioners available in the PETSc library were implemented and evaluated. The efficiency of the solution process was demonstrated with applications to two-dimensional airfoil meshes, and the impact of the curvature functions on the overall convergence was investigated and quantified.In future work, an investigation will be undertaken to shed light on the underlying characteristics impacting the conditioning of the system to further enhance the solution process. An efficient methodology will also be developed to properly integrate the orthogonality functions with the newly-implemented solution process.Ultimately the methodology presented in this work must be generalized to three-dimensional meshes. The challenges arising from the latter, in particular the substantial computing resources required to generate large meshes, will need to be addressed. To this end, a full parallelization of the solution process must also be implemented, and other solution techniques not included herein, such as multigrid methods, will need to be evaluated.



---

*Source: 7266801-2022-07-08.xml*