Abstract
A series of compact implicit schemes of fourth and sixth orders are developed for solving differential equations involved in geodynamics simulations. Three illustrative examples are described to demonstrate that high-order convergence rates are achieved while good efficiency in terms of fewer grid points is maintained. This study shows that high-order compact implicit difference methods provide high flexibility and good convergence in solving some special differential equations on nonuniform grids.
1. Introduction
As one of the earliest and most commonly used numerical methods for computing discrete solutions of differential equations, finite-difference methods are based on the Taylor series expansion on a set of grids, most commonly, a set of uniformly spaced grids. The advantages of finite-difference methods lie in two aspects. They are simple in formulation and implementation, and they are highly scalable. However, one of the disadvantages is that large stencils are usually required when high-order accuracy is stipulated. In addition, the nominal order of accuracy may not be maintained, and often it decreases when nonuniform grids are used in discretization.
Much development has been achieved in finite-difference methods. A number of papers related to various finite-difference algorithms have been published over the past decades. Kreiss [1] was among the first who pioneered the research on compact implicit methods. Hirsh [2] developed and applied a fourth-order compact scheme to Burgers' equation in fluid mechanics. Adam [3] developed fourth-order compact schemes for parabolic equations on uniform grids. Hoffman [4] studied the truncation errors of the centered finite-difference method on both uniform and nonuniform grids and discussed the accuracy of these schemes after applying grid transformations. Dennis and Hudson [5] studied a fourth-order compact scheme to approximate Navier-Stokes-type operators. Rai and Moin [6] presented several finite-difference solutions for an incompressible turbulence channel flow. They also compared their solutions of high-order upwind schemes with solutions obtained from spectral methods. Their work demonstrated that finite-difference methods can be a good approach for the direct numerical simulation of turbulence. Lele [7] demonstrated the spectral-like resolution of the classical Padé schemes and the flexibility of compact finite-difference methods. He developed compact difference methods for mid-point interpolation and filtering with high-order formal accuracy. Yu et al. [8] solved an unsteady 2D Euler equation with a sixth-order compact difference scheme in space and a fourth-order Runge-Kutta scheme in time, which provide a good example of high-order (both spatial and temporal) finite-difference methods. Mahesh [9] presented high-order compact schemes with both the first and second order derivatives evaluated simultaneously. Gamet et al. [10] developed a fourth-order compact scheme for approximating the first-order derivative with a full inclusion of metrics directly on a nonuniform mesh with variable grid spacing. Using this method, a direct numerical simulation of a compressible turbulent channel flow was conducted. Dai and Nassar [11, 12] developed high-order compact difference schemes for high-order derivatives in time and space in solving a 3D heat transport problem at the microscale. Sun and Zhang [13] proposed a new high-order finite-difference discretizaion strategy based on the Richardson extrapolation technique and on an operator interpolation scheme for solving one- and two-dimensional convection-diffusion equations.
There are several numerical models developed over the past 15 years to describe the origin of the geomagnetism, that is, the observed geomagnetic field is generated and maintained by (strongly turbulent) convective flow in the Earth's liquid outer core (geodynamo theory). In these models, Navier-Stokes-type partial differential equations (PDEs) for rapidly rotating, magnetohydrodynamic (MHD) systems are solved via various numerical schemes. The physical systems are in general three-dimensional. In the earliest approach, pseudospectral algorithms are used to solve the PDES in the spatial domain [14]. While this approach is efficient for shared-memory computing systems, it is not very scalable for distributive computing systems. Later, finite-difference algorithms are then employed on the grid points along the radius (of the spherical coordinate system), while spectral algorithms are used on different spherical surfaces that are defined by the radial grid points [15]. More recently, finite-elements algorithms are used to solve the dynamo problems in the spatial domain [16].
However, the finite-difference algorithms in numerical geodynamo simulation need to be considered carefully due to the boundary layers developed in the fluid core. Because of the rapid rotation, the fluid viscosity in the reference frame corotating with the fluid is very small. In general, the viscous effect in the nondimensionalized equations is defined by the Ekman number , which can be very small (e.g., in the Earth's core, ). Consequently, boundary layers of the thickness develop in the system (e.g., Greenspan [17], Kuang and Bloxham [15]). This type of boundary layers, often called the Ekman layer in rotating fluid mechanics, is common in all dynamo simulation. To resolve the boundary layers, often second-order and fourth-order finite-difference schemes are used with varying grid sizes for computational efficiency and convergence of numerical solutions [15]. Could higher order finite-difference algorithms be used to improve numerical accuracy and computational efficiency?
The objective of this paper is to develop sixth-order compact implicit schemes for solving parabolic equations in geodynamo simulations as in [15, 18]. These compact implicit schemes are constructed of three piecewise nonuniform grids to achieve flexibility and efficiency in discretization. The purpose is to advance compact schemes to handle complex geometries with high-order accuracy. This will help to increase the convergence rate in geodynamo models, which are among the most computationally challenging simulations to date. The specific objective of this research is to acquire compact difference schemes of sixth-order accuracy on piecewise nonuniform grids for the differential equations involved in the geodynamo model [15, 18] in the radial direction.
In this paper, we intend to develop sixth-order compact implicit schemes with both the Jacobian transformation (JTr) and the full inclusion of metrics (FIM) on nonuniform grids. We also describe a new method which combines these two approaches. These methods demonstrate high-order accuracy and fast convergence rates in the numerical solutions. In addition, we wish to show the ability of resolving a thin boundary layer with fewer grid points than the traditional finite-difference schemes with the same or a higher order of accuracy. To be efficient in spatial discretization, zeros of Chebyshev polynomials are adopted as the grid points for computation [19], although other options for nonuniform grids were also explored.
The outline of this paper is as follows. We first show that compact implicit schemes combined with nonuniform grids are able to fully resolve the thin boundary layers arising in geo-dynamo simulations with fourth- and sixth-order accuracy. The differential equation here describes a boundary layer profile and is essentially a two-point boundary value problem. The computational interval is discretized via a nonuniform mesh so as to improve the efficiency in spatial resolution. Secondly, we demonstrate the reliability of these schemes subject to different boundary layers, that is, different boundary conditions. A similar order of accuracy and rate of convergence are obtained with Neumann and Robin boundary conditions. Finally, a series of parabolic equations from a geodynamo model [15, 18] is studied using the compact implicit schemes developed earlier in this paper. A combination of the full inclusion of metrics on a nonuniform grid and the Jacobian transformation is proposed and tested. The convergence rate and accuracy order are studied and compared. The full inclusion of metrics on a nonuniform grid is a good alternative approach to solving certain kinds of differential equations when a Jacobian transformation is not applicable. Based on numerical experiments in this paper, it is found that for different parameters in the governing equation, the best convergence rate is achieved with slightly different compact implicit schemes, that is, different linear combinations of neighboring values. This research demonstrates the flexibility and potential of high-order finite-difference methods on nonuniform grids.
2. Mathematical Problems from Geodynamo Simulation
In Kuang and Bloxham's geodynamo simulation model [15, 18], spherical coordinates are used and is the dimensionless radius. The toroidal and poloidal components of the magnetic field and velocity (in the Earth's fluid outer core) are expanded in terms of spherical harmonics along the spherical surface. The coefficients of the spherical harmonic expansion for the magnetic and velocity field become functions of time and the radial distance from the center of the Earth. These coefficients, denoted by , satisfy the following generic second order parabolic partial differential equation: where is the degree of the spherical harmonic functions, is the source term, and is the nondimensional distance measured from the center of the Earth (scaled by the mean radius at the core-mantle boundary). Along the radial direction, based on the structural characteristics and the electromagnetic properties, the domain is divided into three major subdomains: the inner core (), the outer core (), and the layer . The innermost radius is sufficiently small such that asymptotic conditions can be employed. The dimensionless radii and are at the inner core boundary and at the top of the layer, respectively. A variety of boundary conditions, such as Dirichlet, Neumann, and Robin type, are specified at the boundaries between the subdomains. These boundaries contain large gradients of the state variables and serve as the bridges between the unknown state variables in different subdomains along the radial direction. These large gradients are the characteristics of the boundary layers. Therefore, an accurate understanding of these gradients and also the profiles of the boundary layers is vital to geodynamic simulations.
Because a boundary layer is usually one or two orders of magnitude smaller than the typical length scale of the geometry under investigation, the velocity gradients in shear stress or the temperature gradients in heat flux inside a boundary layer are one or two orders of magnitude larger. Consequently the profile of the boundary layers is of great importance in fluid dynamics, heat transfer, and especially computational geodynamo simulations where large values of velocity, and temperature gradients occur at the inner-outer core boundary (dimensionless radius about 0.3) and the core-mantle boundary (dimensionless radius 1). Therefore, numerical solutions of boundary layer problems and reliable predictions of the thermal gradient and shear stress in boundary layers are vital in predicting thermal load, hydraulic pressure drop, velocity and temperature distributions. In addition, an efficient and reliable resolution of a boundary layer is central to the entire simulation since a boundary layer contains the most important information and also it propogates the information into the interior domain. Usually a boundary layer is the obstacle for high-order accuracy and numerical stability in a time-dependent problem [12, 18].
2.1. Compact Difference Methods on Nonuniform Grids
Nonuniform grids are commonly used in modeling when large gradients are expected, local details are desired, or complex geometries are encountered. Previous research work by Hirt and Ramshaw [20] and Roache [21] demonstrated that the finite-difference approximation on a nonuniform grid has lower order of formal accuracy than on a uniform grid.
Generally speaking, two different approaches are commonly used to implement compact finite-difference methods on nonuniform grids. The first is called the Jacobian Transformation (JTr), which involves mapping a nonuniform grid to a uniform grid by means of a Jacobian transformation (described in Section 2.2) and calculating derivatives on the uniform grid to maintain formal accuracy. finite-difference approximations are therefore constructed on the equidistant grid in the transformed space. For this reason, the overall accuracy depends strongly on the transformation, according to Hoffman [4]. One advantage of JTr is that schemes for approximating high-order derivatives are relatively easy to derive. Another advantage is that JTr promises formal accuracy provided that the transformation is continuous and smooth. The disadvantages of JTr come in two aspects. First, sometimes the transformed equations are even more complicated to solve than the original equations. Second, some transformations can give rise to loss of information at certain locations [10] and result in some computational errors.
The other approach is called Fully Integrated Metrics or Fully Included Metrics (FIMs). With the FIM, discretization of a differential equation and evaluation of derivatives are performed directly on a nonuniform grid based on the Taylor series expansion (see Lele [7]). One advantage of the FIM method is that it is straightforward and does not need to transform the equation from one space to another. Another advantage is that it is a reliable alternative when the JTr does not work well. However, a major disadvantage of the FIM is that the approximations implemented on a nonuniform grid do not yield the nominal order of accuracy that can be achieved on an equidistant grid. Therefore the accuracy of the FIM is not guaranteed. Another disadvantage is that high-order schemes on a nonuniform grid are difficult to derive, especially when a large stencil is involved.
2.2. Benchmark Problem 1
Equation (2.1) describes the time evolution of the coefficients of the spherical harmonic expansions for both the magnetic and velocity fields. To mimic the time-stationary solutions of (2.1) on a nonuniform grid, we need to consider a steady state problem. In addition, to study the profile of a thin boundary layer and to resolve the large gradient inside it, local high resolution is necessary. To demonstrate the ability to fully resolve a boundary layer, a two-point boundary value problem as in what follows is studied, subject to the Dirichlet boundary conditions: at both and . This differential equation describes the profile of a velocity boundary layer. The exact solution of this problem is available for comparison and error analysis. In this problem, is the distance from the origin measured in a nonuniform scale. To avoid potential confusion in notation, is used instead of for the radius at the inner-outer core boundary and is the controlling parameter for the thickness of the boundary layer. The smaller the , the thinner the boundary layer. The range of in this study is from to .
To provide enough resolution and to keep the linear system of manageable size as in the geodynamo simulation [15, 22], the zeros of Chebyshev polynomials were ultimately chosen as the grid points to discretize the interval between the inner-outer core boundary () and the core-mantle boundary (), although other options of nonuniform grids were also attempted. Using the JTr to solve this problem, a nonlinear and smooth transformation was created to map the nonuniform grid into a uniform grid on which high-order compact implicit schemes were developed. This mapping from the nonuniform grid to a uniform grid is given as Via this mapping, is an indirect function of . This is why we kept the partial differentiation symbols in (2.2). The Jacobian of this transformation is A straightforward way of solving (2.2) is to express and as the derivatives with respect to via the Jacobian: Then substitute (2.5), (2.6), and (2.3) into (2.2) to obtain where and are functions of only. Next we use implicit Padé approximations for and to solve for (, ), or alternatively, (, ).
However, in geodynamo models, is needed in other computations. Instead of solving (, ), we need to solve (, ) in the following way. First we let , as in (2.5), and then we express as the first-order derivative of with respect to : Substituting (2.8) into (2.2), the original second-order equation becomes the first-order equation: We rearrange (2.5) into Therefore two first-order differential equations (2.9) and (2.10) form a linear system in terms of the state variables (, ), where . Both (2.9) and (2.10) can be approximated with the same Padé schemes. The following fourth- and sixth-order compact implicit schemes were designed to solve the aforementioned linear system.
For the interior grids, a fourth-order scheme consisting of a 3-point stencil is where the derivative is with respect to at the grid point . An alternative choice is However, within the boundary layer, is almost 0, but is large, therefore it is constructive to use more information from . For this reason, (2.11) is slightly better.
Similarly, for the interior grids, a sixth-order scheme with a 5-point stencil is
The advantage of this scheme is that it gives rise to a diagonally dominant linear system, which is theoretically nonsingular. As to the indexing of the state variables, and were arranged alternatily in the sequence , to reduce the bandwidth of the matrix.
It is a common practice to use schemes of one order lower in accuracy for the boundary points than for the interior points. However, in order to resolve the thin boundary layer with high-order accuracy, we maintain the same truncation order. At the two-end points, the following fourth-order compact schemes were used to approximate the first-order derivative:
Correspondingly, the sixth-order schemes for the end points, , and were derived. As the spatial resolution increases, values in adjacent rows near the end points become close. This causes the matrix to be close to ill-conditioned and therefore it is difficult to invert. To avoid catastrophic cancellation in partial pivoting, different schemes were used for adjacent points near the end points. At point 1, (2.15) was used to approximate the derivative in (2.9): At point 2, (2.16) was used for the derivatives in (2.9) and (2.10): At point , (2.17) was used for the derivatives in (2.9) and (2.10): At point , (2.18) was used for the derivative in (2.9):
Figure 1(a) shows three sets of solutions: the exact solution and the approximate solutions from the sixth- and fourth-order compact schemes. No visible difference is observed at this scale. On the left is a thin boundary layer with a velocity gradient about two orders of magnitude larger than the solution itself. On the right is a laminar boundary layer with a mild velocity gradient. Figure 1(b) presents the solution errors (exact solution-approximate solution) in the sixth- and fourth-order compact schemes. Even though there is a large velocity gradient inside the thin boundary layer at the left end, both schemes are able to capture the characteristics of the boundary layer, with the error in the sixth-order scheme reduced to about of that in the fourth-order approximation.
(a)
(b)
Figure 2 shows the convergence rates of the fourth-order compact scheme for and . The convergence rate was calculated from the ratio of the logarithm of the norm of an error (either in the solution or the solution's derivative) and the logarithm of the number of grid points, using the data points after an initial transient region in the plot. On the left is the logarithm of the norm of the solution error versus the logarithm of the number of grid points; on the right is the logarithm of the norm of the error in the derivative of the solution versus the logarithm of the number of grid points. Beyond a threshold of , the convergence rates for the solution are at and at . The convergence rates for are and , respectively.
(a)
(b)
Figure 3 presents the convergence rate of the sixth-order compact scheme for and . The left plot is the logarithm of the norm of the error in solution versus the logarithm of the number of grid points. The right plot is the logarithm of the norm of the error in the derivative of the solution versus the logarithm of the number of grid points. Beyond a threshold of , the convergence rates for the solution are at and at . The convergence rates for are and , respectively.
(a)
(b)
Error analysis indicates that at low-to-medium resolution, that is, about 50 to 250 grid points, as in Figure 1(b), the errors for both and are mainly uniformly distributed over all grids. However, at higher resolution with the number of grid points over , the leading error shifts into the thin boundary layer at the left end. To enhance the local resolution for the derivative in the thin boundary layer, more points within the stencil should be included. For example, for a sixth-order scheme with a 5-point stencil, to enhance the resolution for the derivative of , it would be helpful to use up to 5 points in approximating , as opposed to (2.15) and (2.16).
2.3. Benchmark Problem 2
The general form of a solution is determined by the governing differential equation. However, specific types of boundary conditions also influence the solution of a given problem. To explore the numerical solution of a similar problem as in the previous section, but with different boundary conditions, we modify (2.2) slightly: subject to a Dirichlet boundary condition at and a Robin boundary condition at . The exact solution of this problem is also available for comparison and error analysis. The Robin (also called natural) boundary condition at is numerically weaker than a Dirichlet condition because it imposes less restriction on the solution itself. Because of this change in the boundary condition, the profile of the solution changes accordingly. We use the JTr method, keep the same discretization and the compact difference schemes as in Section 2.2.
Figure 4(a) shows the exact solution and the sixth-order approximate solution of the problem in (2.19) with Dirichlet and Robin boundary conditions. Due to the weak boundary condition at the right end point, the solution at this point is a finite-value instead of being zero. The boundary condition is somewhat like a slip boundary condition, which causes the solution to depend on at the right end point. Figure 4(b) presents the error of the approximate solution from the sixth-order compact scheme. The leading error is still in the thin boundary layer near the left end point. Except for the left boundary layer, the error is basically uniformly distributed over the interval.
(a)
(b)
Figure 5 presents the convergence rate of the sixth-order compact scheme for and in terms of the norm of the error for (left panel) and (right panel). After a transient region, the convergence rates for the solution are at and at . The convergence rates for are and , respectively.
(a)
(b)
These two examples demonstrate that the proposed compact implicit schemes are able to fully capture the characteristics of these boundary layer problems with the projected order of accuracies under strong- and weak-boundary conditions.
3. Sample Problems from Geodynamo Simulation
Next we study a set of sample problems that arose in geodynamo simulation. This section is arranged as follows. First, the problem is set up mathematically and its exact solution is derived. Then the numerical solutions of these geodynamo sample problems are presented using two different methods—Jacobian transformation and fully integrated metrics. After that, a combined method is proposed and tested for its performance in solving the same problems. Numerical results showing the sixth-order accuracy from both methods are presented and comparisons are made afterward.
3.1. Problem Setup
After demonstrating the accuracy and flexibility of high-order compact schemes on nonuniform grids in the previous benchmark problems, we turn our attention back to the geodynamo differential equation (2.1). The Crank-Nicolson method was used in the temporal discretization to maintain second-order accuracy in time while the forcing term was treated explicitly. After applying some spatial discretization to be discussed, (2.1) becomes a set of algebraic equations with the parameter involved, as shown in the following linear system: where and are matrices resulting from the temporal and spatial discretization with the Crank-Nicolson scheme and compact implicit schemes; is the corresponding forcing vector and represents the unknown state variables. The solution at the time level can be computed recursively from the previous time level as where is the amplification matrix: Based on von Neumann stability analysis, the spectral radius of is no more than one due to the Crank-Nicolson method. Therefore an unconditional stability in time integrations is guaranteed. As long as the data and are explicitly given in the initial condition, we are able to acquire the values for the unknown state variables at any time via time marching with a temporal error on the order of the square of the time step (). So far we have addressed the stability issue for (2.1). If we have the property of consistency for this problem, based on the Lax Equivalence theorem, we will be able to find convergent solutions.
Now we discuss the consistency, that is, we address the question how to form the matrices and , including the boundary conditions. To this end, we eliminate time and consider an equivalent steady-state (time-independent) problem for (2.1): in which the forcing term is defined by where , , and . The aforementioned three intervals in (3.5) correspond to the inner core (IC), the outer core (OC), and the layer (). is the degree of the spherical harmonic expansion. Because the magnetic monopole does not exist in nature, the magnetic flux starts from the dipole component and ranges from 2 to 33 in [15, 22].
The boundary conditions for (3.4) are two Robin conditions:
Equations (3.4), (3.5), and (3.6) define a series of problems of the energy spectrum in the spectral space in terms of the parameter . Exact solutions can be obtained analytically so as to perform error analysis and detailed comparison. The details of the derivation of the exact solutions of the series of problems is described in the appendix.
3.2. Sixth-Order Approximations With JTr Method
As in Section 2.2, the zeros of Chebyshev polynomials are used to discretize the inner core, the outer core, and the layer. A general transformation between the uniform grid and the nonuniform grid is given as The associated Jacobian of this transformation is Since is of special interest in geodynamo models [15, 22], we would like to keep this quantity as a state variable. Based on the Chain Rule, after introducing , we split (3.4) into two first-order differential equations as in what follows in terms of and : Both equations in (3.9) can be approximated with the same scheme at a guaranteed order of accuracy in the transformed space.
The boundary conditions do not need to be transformed since we solve for and here, and and appear naturally in the boundary conditions. The state variables and were arranged alternately in the sequence , in the state variable index to reduce the bandwidth of the matrix as in Section 2.2.
In terms of the approximation schemes, one approach, similar to Section 2.2, is described in what follows. For the interior points, we used (2.13) to approximate the first derivative. This numerical scheme leads to a symmetric matrix since the stencil is centered. For the two end points near , we used and for the two end points near , we used Results from this approach show stable and consistent sixth-order convergence as anticipated.
Considering and are arranged alternately, an asymmetric approximation could also perform well if used properly. In the second approach, (3.10) with indices shifted backward by one was used for interior points. For the two end points near , we used For the two end points near , we used
The second approach turns out to be slightly better than the first one in decreasing the norms of the solution error and the derivative error, called the overall error in this paper, by a factor of about 2 to 4. Therefore, the numerical results based on the JTr are obtained by the second approach. Since there are three intervals from to , the number of grid points in each interval affects the overall error. Therefore, the simulation result with the smallest overall error is obtained when the errors in these three intervals are approximately the same order. In fact, there are many different combinations of the way that grid points are distributed in these three intervals at a given total number of grid points. The subsequent results are given at the smallest overall error. The numerical results indicate that at different total numbers of grid points, the minimum overall error is achieved at slightly different grid point distributions among these three intervals.
3.3. Sixth-Order Approximations with Combined Method
We have demonstrated the fourth and sixth-order accuarcy and consistency in convergence for using the JTr in two benchmark problems. A proper combination of the JTr and FIM is expected to give more flexibility and maintain the same order of accuracy. This combined method is believed to be especially suitable for the unique problems in geodynamo modeling with three nonuniform discretization regions and the radius being close to zero on the left end point in (3.9). We demonstrate that this combined method is capable of generating reliable numerical solutions for the geodynamo sample problems as in (3.4).
Technically, this combined method integrates the JTr and the FIM, then it applies to different grid points. For the same geodynamo sample problem, the JTr was used at most of the grid points except the two adjacent points near the end . This will provide good formal accuracy because of the high-order compact schemes developed in the previous sections. For these two adjacent points near the end , where is close to zero and is very large, the FIM was used to reduce the local truncation error, since the JTr may introduce a relatively large error in case of low resolution. It is unnecessary to use the FIM for other points including the end point on the other side of the inner core, the end points of both the outer core, and the layer. Our study showed that at these points, the FIM does not noticeably help to reduce the discretization error.
Next we give the details of the combined method for all points in three intervals. For the end point in the inner core closest to and the end points in the layer, the combined method uses (3.10) and (3.11). For the interior points in both the inner core and the layer, it uses (2.13). For the end points in the outer core, it uses (3.12) and (3.13). For the interior points in the outer core, it uses (3.10).
For the two grid points nearest to , we develop the following numerical schemes based on the FIM. First, by introducing , (3.4) can be reduced into two first order differential equations: To approximate the first-order derivative near , a 3-point stencil difference approximation based on the FIM is We define the radii at the grid points , , and as , , and , respectively. Then we expand every term at the point with the Taylor series expansion and combine terms of the same order. By matching terms of the same order on each side, we obtain
Since the grid spacing is much smaller near the end point than in the middle, that is, is smaller than , the truncation error will actually be smaller than the corresponding error in a uniform grid at the same number of grid points. Therefore a fourth-order scheme will yield slightly higher than fourth-order norminal accuracy. We chose a fourth-order scheme with a 3-point stencil such that the first five equations in (3.17) are satisfied exactly. The leading order of the truncation error is The first five equations in (3.17) form a closed linear system for five unknowns out of the six unknown coefficients in (3.15), with one degree of freedom. For point 1, we choose so that (3.15) becomes After we define , , and , the coefficients are For point 2, we choose , and (3.15) becomes where the coefficients are
We obtained numerical results from purely JTr for all grid points in three sections, and also results from the combined method. Although the two methods are not completely different, numerical results show dramatic differences.
Figure 6 presents a series of numerical solutions of the sample geodynamo problem (3.4) with the JTr method. The left plot is the solution ; the right plot is the derivative of the solution . It is observed from the left plot that the dipole component () is dominant; it contains over 90% of the total energy of . Therefore the magnetic field of the Earth behaves like a giant magnet with two opposite polarities although it is actually a superposition of magnetic multipoles. As becomes larger, the magnitude of decreases rapidly. For example, at , reduces to less than 10% of the value at . When , this magnitude drops to less than 1%. There is a saddle point at the inner-outer core boundary near where the derivative decreases on the left side and increases on the other for all the values. Therefore is not smooth and is discontinuous at this point. This is also why we solve the first derivative directly instead of the second derivative in the formulation.
(a)
(b)
Figure 7 shows the local error distribution across the entire region for from the results with the JTr and the combined method. The state variable is on the left and the derivative of the state variable is on the right. The resolution was set at 100 grid points. The magnitude of the error for both and are about 7 orders of magnitude less than the magnitude of the variables themselves. Therefore there is no visual difference between the numerical solution and the exact solution, and the solution plots were omitted. However, in the left plot, the leading error for is located around the inflection point near . Errors at the left boundary within the inner core and at the right boundary in the layer are also slightly larger than the neighboring points. In general, the JTr and the combined method produced comparable orders of accuracy. In the right plot, the leading error for is at the left boundary in the inner core. The JTr produces a slightly higher error than the combined method.
(a)
(b)
Figure 8 demonstrates the local error distribution for from the results obtained with the JTr and the combined method. The error in solution variable (a) is on the left and the error in the derivative of the solution variable (b) is on the right. The resolution was set at 190 grid points. The leading errors for both and have shifted to the location around the inflection point. Additionally, the leading errors increased slightly with increased resolution. This is because as increases by one unit, the coefficient of , , which is called the stiffness of (3.4), almost doubles its value. As is much less than one in most grid points, this will make the stiffness increase very fast. Therefore (3.4) becomes very stiff, and the discretization matrix equation (3.3) becomes close to ill-conditioned, and even close to singular as the grid points increase over 1000. However at current CPU capacities, it is still difficult to use 1000 grid points in the radial direction in the MoSST model [15, 18]. Both plots in Figure 8 show that the results from the combined method have larger leading error than the result from the JTr around the inflection point with approximately 0.34. This is because the FIM schemes developed in this paper are nominally of fourth-order accuracy, which is two orders lower than the sixth-order accuracy schemes based on the JTr used for the rest of the points. It is very difficult to develop a sixth-order compact scheme based on the FIM since for the full stencil of 5 points, 10 unknown coefficients will be involved in a general notation and 7 equations will be solved with 3 degrees of freedom.
(a)
(b)
As becomes large, both methods become less efficient although they are still able to provide accurate numerical solutions. The convergence properties of both methods deteriorate, especially the combined method since it uses fourth-order schemes on the left end. Figure 9 demonstrates the convergence properties of both methods at . Since the value of is small, the stiffness of (3.4) is very small and both the JTr and the combined method provide excellent convergence rates of above , as listed in Table 1. As the spatial resolution becomes higher than 250 grid points, the coefficients in the approximation schemes at adjacent points based on the combined method tend to be close in value. Therefore the matrix generated from the combined method tends to be close to ill-conditioned earlier than that from the JTr.
(a)
(b)
As becomes even larger, the stiffness of (3.4) becomes more than the square of . This large stiffness makes (3.4) hard to solve and the discretized linear system close to singular. Therefore, the convergence rate drops noticeably, especially for the combined method. For example, when , the combined method has a convergence rate 6.05 for and 6.01 for . But when , the convergence rates for the combined method become 4.34 and 4.54, respectively, as shown in Figure 10. As the number of grid points increases beyond 250, the error in the combined method starts to grow as shown in Figure 10. This is because the schemes developed from the FIM with the nominal fourth-order accuracy are becoming incompatible with the sixth-order schemes developed from the JTr. On the other hand, the results from purely JTr demonstrate to be robust and maintain the expected order of convergence.
(a)
(b)
4. Conclusions
The benchmark differential equations involved in this paper describe parametrically adjustable boundary layer profiles. These problems are essentially similar to a two-point boundary value problem. In addition, a series of parabolic differential equations that arise from the spectral discretization of a geodynamo model is studied with the compact implicit schemes developed in this paper.
The sixth-order compact implicit schemes with both Jacobian transformation and full inclusion of metrics on nonuniform grids for better efficiency are developed. Additionally, a new method which integrates these two approaches is developed and tested. These methods demonstrate high-order accuracy, fast convergence rate, and robust adaptibility to Dirichlet, Neumann, and Robin boundary conditions. In addition, it is demonstrated in the numerical results that high-order compact methods with the JTr and FIM techniques are capable of resolving a thin boundary layer with relatively fewer grids than the traditional finite-difference schemes with the same order of accuracy. Convergence rate and accuracy order are studied and compared in detail.
The reason that the result from the JTr is slightly better than that from the FIM is because as approaches , where , the Jacobian () is close to zero as in (2.4). Therefore, is one magnitude smaller than although both are much larger than one at . Consequently, the stiffness in (3.9) is about one order less than that in (3.14). However, the combined method which uses the FIM directly on a nonuniform grid is always a good backup and alternative approach in solving certain kind of differential equations when the JTr is not applicable or results in loss of information.
Based on these numerical experiments, it is found that for different numbers of total grid points in the inner core, the outer core, and the layer, the best convergence rate is obtained with slightly different combinations of the numbers of grid points in all intervals. This also demonstrates the high flexibility and potential of the sixth-order compact finite-difference schemes developed in this paper. These methods provide some reference for improving the ongoing development of the geodynamo model.
Appendix
Exact Solutions
The exact solutions of the series of problems, described by (3.4), (3.5), and (3.6) in geodynamo modeling [15, 22], and the derivatives, have the following form: where is the particular solution of the differential equation: After solving the aforementioned equation, the particular solution has the form: The coefficients , , , and in (A.1) can be obtained from the connectivity relations: Equation (A.4) is a system for these four coefficients. The explicit expressions for these coefficients are given as where , , , and can be computed from (A.3) as
Acknowledgments
This work was supported by the NSF Mathematical Geophysics Program under the grant EAR-0327875. Our appreciation also goes to our colleagues including Dr. Weiyuan Jiang and Dr. Bernd Schroeder for providing computing support, constructive suggestions, and proofreading.