Received 15 February 2015; accepted 26 March 2015; published 30 March 2015

ABSTRACT

A new numerical approach, called the “subdomain Chebyshev spectral method” is presented for calculation of the spatial derivatives in a curved coordinate system, which may be employed for numerical solutions of partial differential equations defined in a 2D or 3D geological model. The new approach refers to a “strong version” against the “weak version” of the subspace spectral method based on the variational principle or Galerkin’s weighting scheme. We incorporate local non- linear transformations and global spline interpolations in a curved coordinate system and make the discrete grid exactly matches geometry of the model so that it is achieved to convert the global domain into subdomains and apply Chebyshev points to locally sampling physical quantities and globally computing the spatial derivatives. This new approach not only remains exponential convergence of the standard spectral method in subdomains, but also yields a sparse assembled matrix when applied for the global domain simulations. We conducted 2D and 3D synthetic experiments and compared accuracies of the numerical differentiations with traditional finite difference approaches. The results show that as the points of differentiation vector are larger than five, the subdomain Chebyshev spectral method significantly improve the accuracies of the finite differ- ence approaches.

Keywords:

Numerical Differentiation, Chebyshev Spectral Method, Curved Coordinate System, Arbitrary Topography

1. Introduction

With development of computer technology, computer modeling that imitates the underground earth interior becomes now a powerful tool to geological, environmental and geophysical problems, e.g. numerical simulations of tectonic dynamics [1] [2] , underground water flow or geothermal system [3] -[5] , seismic wave propagation [6] -[9] and geo-electromagnetic fields [10] [11] . Mathematically, computer modeling seeks approximate solutions of partial differential equations (PDE) that govern geological evolution or geophysical fields in the Earth, subject to a Dirichlet or Neumann boundary conditions. The simplest method of computer modeling is finite- difference approach to the partial derivatives of PDE, e.g. various high-order finite-difference formulae [12] and staggered grid schemes [13] . Due to simple principle and easy implementation, they become now widely used for many geological and geophysical problems [14] [15] . However, these methods often employ regular grids, and consequently suffer from complex geometries of geological models, e.g. various free-surface topography, such as hill, trench and folded rocks, which are often encountered in practice. In these cases, the finite-difference approach has to refine the free-surface topography and subsurface interfaces by stepwise curves or stepwise blocks [16] , which cannot exactly match the geometries of geological models. To overcome the disadvantage, many researchers prefer the so-called “weak” solutions of PDE, i.e. finite-element method [17] or subspace spec- tral element method [18] [19] , which is based on the variational principle or Galerkin’s weighting scheme. Alternatively, one has to define a set of coordinate transformations and adapt the finite-difference approach to a curved coordinate system so as to match the free-surface topography and subsurface interfaces of geological models [20] -[22] .

Spectral method, e.g. Chebyshev spectral method is superior to the traditional finite-difference approach in numerical differentiations because of exponential convergence [23] . Due to its high accuracy, easy computer programming, and applicability to complex model geometries, the spectral method becomes one of the most popular accurate solvers of PDE [23] -[25] . However, the most disadvantageous aspect of the spectral methods is expensive consumption of computer memory and CPU time, because it employs global collocation points in calculations of the spatial derivatives, resulting in a fully filling-in differentiation matrix [23] . Therefore, it suffers the disadvantage in computer modeling for a large 3D geological or geophysical model [24] [25] .

In this paper, we present a subdomain Chebyshev spectral method to calculate spatial derivatives in a curved coordinate system, which is similar to subspace spectral element method (weak form) [18] that divides the global domain into subdomains (elements) exactly matching the geometry of the model, and then Chebyshev points are applied to discretization of physical quantities and calculation of the spatial derivatives in each subdomain. The Chebyshev collocation points for the spatial differentiations only involve the local physical samples instead of the global values; therefore, the method saves huge computer memory for a 3D computer modeling. Such discretization results in a sparse assembled matrix like the “weak” version of subspace spectral element method [18] [19] and Gaussian quadrature grid approach [9] . We conduct 2D and 3D synthetic experiments and compare the results with the finite-difference approaches. The results show that the new method is superior to the traditional finite-difference method in numerical differentiations and applicable to complex geological models having arbitrary free-surface topography and multiple subsurface interfaces.

2. Strong Solution

In general, computer modeling is to solve the following PDE [26] :

, (1)

subject to some Dirichlet or Neumann boundary conditions. Here is often a second-order linear differential operator depending on model vector m and partial derivatives and (i, j = 1, 2, 3). The vectors and represent location and physical field in a 2D or a 3D model domain Ω. The right-hand side vector is considered as a source of the physical field F defined in Ω. The vectors F and m are different physical quantities according to geological or geophysical problems, e.g. in tectonic dynamics, F may consist of pressure field p, deviatoric stress and velocity field v, and m comprises the density and viscosity: of the Earth [1] [2] ; in seismic wave modeling, F is the ground displacement vector u, or composed of speed v and stress tensor, and m includes density and elastic moduli: [8] [9] ; in geo-electromagnetics, F becomes either electric field E or magnetic field H, and m involves electric permittivity, magnetic permeability and conductivity: [7] [10] [11] . All these quantities are functions of the spatial coordinates.

In order to match free-surface topography, e.g. hill, trench and ridge, one may employ the global coordinate transformations [20] -[22] :

(2)

mapping the computational coordinates into the physical coordinates. In general, we assume F be a component of the vector F, and according to the chain law, we have the first and second derivatives:

, (3)

. (4)

Here, summation convention of the repeated integer subscripts k and l has been applied. Accordingly, the discrete versions of Equations (3) and (4) can be written in the following forms

(5)

(6)

where

(7)

(8)

Here, (or) and are components of the Jacobian and Hessian matrices respectively, both can be analytically calculated from the given Equation (2); and are the first and second differentiation vectors and have the following generalised forms:

(9)

Here, and stand for components of the vectors and; and are the number of points in the -direction and -directions, e.g. in Cartesian system, if Equation (2) is, then the Jacobian and Hessian matrices become and respectively, so Equations (7) and (8) are

reduced to the simplest forms: and, which become the standard finite-difference or

spectral differentiation vectors [23] . According to Equation (9), the discrete forms of the field quantity F are ob-

tained, i.e., and, (, ,

), corresponding to and, respectively.

Substituting Equations (5) and (6) into (1), and then applying the PDE to every point results in 3N linear equations:

, (10)

which have 3N unknowns of and gives a linear system. Solving the linear system subject to Dirichlet or Neumann boundary conditions, one obtains the so-called “strong” numerical solutions: due to the same smoothness of as F defined in Equation (1). If and are determined by the neighbour points of, e.g. using a fourth-order finite-difference approach, the assembled matrix of Equation (10) becomes sparse and saves much computer memory. Then, the linear system can be efficiently solved by iterative or parallel solvers for a large 3D geological model [27] [28] . Therefore, the crucial step of geological or geophysical modeling is to find accurate and efficient differentiation vectors and in terms of Equations (7) and (8).

It should be noticed that the differentiation vectors and given by Equations (7) and (8), are appli-

cable only if the Jacobian and Hessian matrices are calculable in Ω. This means that Equ-

ation (2) must be smooth (differentiable) in the domain. In the next section, we present a set of the local coordinate transformations which are applicable for arbitrary free-surface topography and make the numerical differentiation vectors and calculable.

3. Coordinate Transformation

In order to exactly match free-surface topography and subsurface interfaces of the Earth, we divide the model domain into non-overlapping subdomains, e.g. in 3D case, we have , where and are the free-surface topography and subsurface interfaces if they are present. If there is not any physical interface in the model domain, are pure mathematical boundaries of the subdomains. In each subdomain, we write Equation (2) into the following forms:

(11)

which maps the local coordinates into the global coordinates. Here and represent the size and central point of the subdomain, respectively, e.g. in the z-direction, we have

(12)

and both are functions of x and y. Therefore, Equation (11) gives non-zero components of the Jacobian matrix:

(13)

and non-zero components of the Hessian matrix:

(14)

Equations (13) and (14) indicate that the Jacobian and Hessian matrices depend on the derivatives, , and. According to Equation (12), these derivatives require the subdomain

boundaries to be differentiable in the global domain. For this requirement, we apply cubic-

spline interpolation functions to the boundaries [29] [30] , e.g. in a 3D case,

, (15)

where are coefficients defined in the subdomain and determined from the known samples

at the boundaries of the subdomains. Spline interpolation theory has shown that arbitrary scattered or regularly-distributed free-surface topography or subsurface interface can be expressed by Equation (15), and these expressions are differentiable everywhere in the model domain. Therefore, Equation (15) guarantees the partial derivatives, , and for computing the Jacobian and Hessian matrices involved in Equations (13) and (14).

4. Subdomain Chebyshev Spectral Method

In the transformed subdomain, we apply Chebyshev points [23]

(16)

to three coordinates, where, and are numbers of the points in the three directions. Geo- metrically, these points are visualized as the projections on of equispaced points on the upper half of the unit circle. Figure 1 gives the two subdomains of the 1D case. One can easily extend it to 2D and 3D cases. Substituting the Chebyshev points into equation (11), one obtains the global coordinates that de- monstrate a one-to-one mapping between the two coordinate systems. All these Chebyshev points form a grid that discretizes the model parameters m and physical field vector F in Equation (1). This grid may be called “subdomain Chebyshev differentiation grid”, which is similar to the “Gaussian quadrature grid” we used to sample the physical quantities in the elements [8] [9] . Basing on the subdomain Chebyshev differentiation grid, we construct the Lagrange interpolation polynomials with the Chebyshev points in the subdomain, so that we have the first differentiation vectors

(17)

and second differentiation vectors

(18)

where are cardinal functions, and are their first and second derivatives. These differentiation vectors have exponential convergences of numerical differentiations in the subdomains and avoid the Runge phenomenon¾the errors increase exponentially in the equispaced case [23] .

Substitutions of Equations (13) and (17) for (7), we obtain general forms of the first differentiation vectors with respect to the global coordinates:

(19)

Similarly, substituting Equations (13), (14) and (19) into (8), we obtain the second differentiation vectors with respect to the global coordinates. These substitutions are more or less tedious and involve computations of the second derivatives, and for and. We omitted them here for saving space. However, the following simple method may be used to calculate these second differentiation vectors, e.g. according to Equation (5), we have following identity:

(20)

and obtain

(21)

Here, we used and to stand for their components so as to clearly indicate the

relationship between and. Calculating, gives the same result as Equation (21). This shows that the second differentiation vector given by Equation (21) satisfies Clairaut’s theorem, and can be obtained by multiplication of two first differentiation vectors and. The advantage of Equation (21) is avoidance of computing the second derivatives of coordinate transformations, so that computation of the second differentiation vectors becomes much simpler than using Equation (6).

It should be noticed that Equations (19) and (21) are valid for the points inside the subdomains. At a boundary point of the subdomains, they may yield different values of the partial derivatives due to multiple subdomains’ sharing. For the “strong” solution, the partial derivatives at the boundary points must be unique, for which we pick up half of Chebyshev points from the connecting subdomains (see points in the “boundary domain” in Figure 1) and construct the Lagrange interpolation polynomials given by Equation (17) and (18). Consequently, Equations (19) and (21) become well defined (uniqueness) for the partial derivatives at the boundary points. Apparently, the differentiation vectors at the boundary points differ from the ones at the points inside the subdomains due to the different distributions of Chebyshev points in the boundary domains and the ordinary subdomains (see Figure 1). This difference does not break down Equations (19) and (21), and on the contrary, it makes them applicable to any point in the whole domain. We call these differentiation schemes the “subdomain Chebyshev spectral method”, because it applies the Chebyshev points to the subdomain rather than the globe domain presented in standard spectral method [23] . Particularly, if we replace the Chebyshev points with equispaced points in the subdomains, Equations (19) and (21) become the central finite-difference approaches. It shows that the subdomain Chebyshev spectral method may regress to the finite-difference approaches.

5. Numerical Experiments

To demonstrate the capability of the new approach, we designed a 2D and a 3D geological model (see Figure 2).

The former is a valley and the latter represents a ridge, both have two interfaces under the ground. With the given model geometries, we used subdomain Chebyshev points to discretise the model domain and obtained the subdomain Chebyshev differentiation grids. Figure 2 shows the grids of the two models, from which one can see 10 ´ 8 and 10 ´ 10 ´ 10 bent subdomains of the 2D and 3D geological models, respectively. The grids exactly match free-surface topography and subsurface interfaces.

In order to examine the accuracy of the subdomain Chebyshev spectral method, we employed the following testing functions as the physical quantities defined in the two models:

(22)

(23)

Here, L_x, L_y, and L_z are lengths of the model domain in three coordinate directions. From the two functions, one can analytically calculate the first and second partial derivatives, which are used to check the accuracies of the subdomain Chebyshev spectral method. Figure 3 gives five partial derivatives (left column) and their absolutely relative errors (right column) of the 2D model, computed by the presented method. In these experiments, we applied seven Chebyshev points in the horizontal and vertical directions in each subdomain, and 11 points for the 1D spline interpolation functions matching the free-surface topography and subsurface interfaces. The left column shows that and are odd functions of x;, and are even functions of x; and equals, all of which are predicted from exact derivatives of Equation (22). The right column demonstrates that accuracies of the derivatives with respect to the vertical direction, i.e. and are much higher than the derivatives with respect to horizontal direction, i.e., and. This is because of the linear transformation from to z and non-linear relationship between and x in Equation (11), which defines the changes of free-surface topography and subsurface interfaces in the horizontal direction. However, one can find that all absolutely relative errors are less than 0.3%. These results prove that the subdomain Chebyshev spectral method can provide accurate spatial derivatives in a 2D geological model.

In order to compare the subdomain Chebyshev spectral method with the traditional finite-difference approach, we repeated the above numerical experiments with finite-difference method and compared the average and maximum absolutely relative errors of the two methods. Meanwhile, five algorithms have been implemented¾ three finite-difference approaches (FDM0, FDM1, FDM2) and two subdomain Chebyshev spectral methods (SCSM1, SCSM2). Here, the number “0” denotes the algorithm applying Lagrange interpolation functions to the free-surface topography and subsurface interfaces instead of spline interpolation functions. The integers “1” and “2” stand for the algorithms using Equation (8) and (21), respectively, for computing the second differentiation vectors, i.e., and. As mentioned above, the difference between the two algorithms is inclusion or exclusion of the high-order derivatives of the free-surface topography and subsurface interfaces, i.e., and.

Figure 4 shows the average absolutely relative errors of five partial derivatives computed by the subdomainChebyshev spectral methods and finite-difference approaches. From these results, four facts are observed: firstly,

replacement of the global spline interpolations with local Lagrange interpolations, to match free-surface topography and subsurface interfaces, will lead numerical differentiation to fail in convergence (see FEM0-curves in -diagram, -diagram and -diagram of Figure 4). This is because Lagrange interpolations of the topography and subsurface interfaces cannot guarantee differentiability of the coordinate transformations defined by Equation (11) [29] , which are required in the computation of the differentiation vectors and given by Equations (7) and (8); secondly, the convergence curves of the derivatives and computed by the two methods are slightly different (see -diagram and -diagram in Figure 4). This is because of linear transformation from to (see Equation (11)); thirdly, the two algorithms of FDM2 and SCSM2 performs

better convergences than the algorithms of FDM1 and SCSM1 in computations of the second derivatives, because the former does not require high-order derivatives and actually reduces the order of smooth- ness in the coordinate transformations; finally, increasing the points in the subdomains (larger than five), the subdomain Chebyshev spectral methods significantly improve accuracies of the finite-difference approaches, and show achievement of almost one-order higher accuracies than the finite-difference methods.

In the 3D case (see Figure 2), we once again applied the subdomain Chebyshev spectral method to compute nine partial derivatives and their absolutely relative errors. Figure 5 shows an example in which seven points were also employed for the subdomain Chebyshev differentiation vector, and 10 ´ 10 points were applied to 2D cubic spline interpolations defining the free-surface topography and subsurface interfaces. In this figure, we give two panels: (a) and (b), presenting the computed derivatives and absolutely relative errors, respectively; both have three sections at y = −20, 0, +20 km showing the patterns of the derivatives and distributions of the absolutely relative errors in the model domain. From these results, one can see that patterns of the numerical results are completely consistent with the analytic derivatives, e.g. Figure 5(a) show that is an even function of x and y and that is an odd function of x and y; meanwhile, equals. Figure 5(b) shows that most of the errors appear at the boundaries of the subdomains and places of the rough free-surface and subsurface interfaces; the maximum absolutely relative errors are less than 0.7% for the second derivatives. These performances show that the subdomain Chebyshev spectral method is an accurate tool for spatial derivatives in such a 3D geological model.

Figure 6 gives overall comparisons of the subdomain Chebyshev spectral method with the finite-difference approaches in the 3D case. The average absolutely relative errors of nine spatial derivatives are plotted against points of the differentiation vectors. These results exhibit similar performances to the 2D case, e.g. the conver- gence curves of and have slight difference between two methods due to the linear transformation

from to z. Other curves demonstrate that, as the points of differentiation vectors are larger than 5, the subdomain Chebyshev spectral method significantly improve accuracies of the finite-difference approaches. Therefore, through the 3D experiments, we once again show that the subdomain Chebyshev spectral method is superior to the finite-difference approximation in numerical differentiation with curved coordinate system.

6. Conclusions

We have demonstrated a subdomain Chebyshev spectral method to numerical differentiations in a curved coordinate system, which may be employed to seek a “strong” solution of PED defined in ageological or a geophysical model. The new method employs a set of non-linear local coordinate transformations, global spline interpolations and subdomain Chebyshev points to make the discretization grid automatically match free-surface topography and subsurface interfaces, and guarantees the numerical solution converging to the “strong” solution of partial differential equations. The new method possesses advantages from the traditional finite-difference approaches, i.e. simple discretization procedure and sparse assembled matrix, but it has stronger capability to deal with complex geometries of geological models and better accuracies in numerical differentiations. Methodologically speaking, the subdomain Chebyshev spectral method may replace the finite-difference approaches in geological or geophysical computer modelling and enables us to easily hand free-surface topography and significantly improve accuracy of modelling results.

2D and 3D synthetic experiments demonstrate the accuracies and capabilities of the subdomain Chebyshev spectral method. As points of the differentiation vectors are large than 5, the subdomain Chebyshev spectral method can reach one-order higher accuracy in numerical differentiations than the traditional finite-difference approaches in curved coordinates systems. Two algorithms of the subdomain Chebyshev spectral method (SCSM1, SCSM2) are both valid for numerical differentiations in PDE, but the algorithm SCSM2 is more efficient in computations than SCSM1, due to exclusion of computing the high-order derivatives of the coordinate trans- formations.

Acknowledgements

This work was supported by a Discovery Project (DP1093110) of the Australia Research Council. The authors thank Mr. Craig Patten for his assistance in using high-performance computing facility at eResearch SA in Australia.

Cite this paper

Bing Zhou, Graham Heinson, Aixa Rivera-Rios,, (2015) Subdomain Chebyshev Spectral Method for 2D and 3D Numerical Differentiations in a Curved Coordinate System. Journal of Applied Mathematics and Physics, 03, 358-370. doi: 10.4236/jamp.2015.33047

References