High-Order Compact Implicit Difference Methods For Parabolic Equations in Geodynamo Simulation

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.


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

Mathematical Problems from Geodynamo Simulation
In Kuang and Bloxham's geodynamo simulation model 15, 18 , spherical coordinates are used and r is the dimensionless radius.The toroidal and poloidal components of the magnetic field B and velocity V 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 L is the degree of the spherical harmonic functions, f is the source term, and r 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 r 0 ≤ r < r i , the outer core r i ≤ r ≤ 1 , and the D layer 1 < r ≤ r d .The innermost radius r 0 is sufficiently small such that asymptotic conditions can be employed.The dimensionless radii r i and r d are at the inner core boundary and at the top of the D 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 .

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.

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: φ 0 at both r r 0 and r 1.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, r is the distance from the origin measured in a nonuniform scale.To avoid potential confusion in notation, r 0 is used instead of r i 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 0.5 to 0.002.
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 r r 0 and the core-mantle boundary r 1 , 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 r ∈ r 0 , 1 to a uniform grid x ∈ 0, π is given as

Advances in Mathematical Physics
Via this mapping, φ is an indirect function of x.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 ∂φ/∂r and ∂ 2 φ/∂r 2 as the derivatives with respect to x via the Jacobian: ∂φ ∂r Then substitute 2.5 , 2.6 , and 2.3 into 2.2 to obtain where ξ and η are functions of x only.Next we use implicit Padé approximations for ∂ 2 φ/∂x 2 and ∂φ/∂x to solve for φ, ∂φ/∂x , or alternatively, φ, ∂ 2 φ/∂x 2 .However, in geodynamo models, ∂φ/∂r is needed in other computations.Instead of solving φ, ∂φ/∂x , we need to solve φ, ∂φ/∂r in the following way.First we let u ∂φ/∂r, as in 2.5 , and then we express ∂ 2 φ/∂r 2 as the first-order derivative of u with respect to x: Substituting 2.8 into 2.2 , the original second-order equation becomes the first-order equation: We rearrange 2.5 into

2.10
Therefore two first-order differential equations 2.9 and 2.10 form a linear system in terms of the state variables φ, u , where u ∂φ/∂r.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 φ i is with respect to x at the grid point x i .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, u and φ were arranged alternatily in the sequence u 1 , φ 1 , u 2 , φ 2 , . .., 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:

2.14
Correspondingly, the sixth-order schemes for the end points, 1, 2, N, and N 1 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 N, 2.17 was used for the derivatives in 2.9 and 2.10 : At point N 1, 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 1/6 of that in the fourth-order approximation.
Figure 2 shows the convergence rates of the fourth-order compact scheme for 0.01 and 0.002.The convergence rate was calculated from the ratio of the logarithm of the L 2 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 L 2 norm of the solution error versus the logarithm of the number of grid points; on the right is the logarithm of the L 2 norm of the error in the derivative of the solution versus the logarithm of the number of grid points.Beyond a threshold of N 40, the convergence rates for the solution φ are 4.1 at 0.01 and 3.9 at 0.002.The convergence rates for ∂φ/∂r are 3.7 and 3.5, respectively.
Figure 3 presents the convergence rate of the sixth-order compact scheme for 0.01 and 0.002.The left plot is the logarithm of the L 2 norm of the error in solution versus the logarithm of the number of grid points.The right plot is the logarithm of the L 2 norm of the error in the derivative of the solution versus the logarithm of the number of grid points.Beyond a threshold of N 40, the convergence rates for the solution φ are 6.5 at 0.01 and 6.0 at 0.002.The convergence rates for ∂φ/∂r are 5.8 and 6.0, respectively.

Advances in Mathematical Physics
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 ∂φ/∂r are mainly uniformly distributed over all grids.However, at higher resolution with the number of grid points over 1000, 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 ∂φ/∂r, as opposed to 2.15 and 2.16 .

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 φ 0 at r r 0 and a Robin boundary condition ∂φ/∂r φ 0 at r 1.The exact solution of this problem is also available for comparison and error analysis.The Robin also called natural boundary condition at r 1 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.
Figure 5 presents the convergence rate of the sixth-order compact scheme for 0.01 and 0.002 in terms of the L 2 norm of the error for φ left panel and ∂φ/∂r right panel .After a transient region, the convergence rates for the solution φ are 6.8 at 0.01 and 7.2 at 0.002.The convergence rates for ∂φ/∂r are 5.6 and 6.0, respectively.
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.

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.

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 L involved, as shown in the following linear system: where A and B are matrices resulting from the temporal and spatial discretization with the Crank-Nicolson scheme and compact implicit schemes; F is the corresponding forcing vector and − → φ represents the unknown state variables.The solution at the time level n 1 can be computed recursively from the previous time level n as where G is the amplification matrix:

3.3
Based on von Neumann stability analysis, the spectral radius of G 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 − → φ 0 and F 0 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 Δt 2 .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 A and B, 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 f r is defined by

3.5
where r 0 0.034286, r i 0.34286, and r d 1.0057143.The aforementioned three intervals in 3.5 correspond to the inner core IC , the outer core OC , and the D layer D .L 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 L ranges from 2 to 33 in 15, 22 .
The boundary conditions for 3.4 are two Robin conditions: L 1 r φ 0, at r r 0 ; and ∂φ ∂r L r φ 0, at r r d .

3.6
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 L. 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.

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 D layer.A general transformation between the uniform grid x ∈ 0, π and the nonuniform grid r ∈ a 1 , a 2 is given as The associated Jacobian of this transformation is

3.9
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 u ∂φ/∂r here, and φ and u appear naturally in the boundary conditions.The state variables u and φ were arranged alternately in the sequence u 1 , φ 1 , u 2 , φ 2 , . .., 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 a 1 , we used Advances in Mathematical Physics and for the two end points near a 2 , we used Results from this approach show stable and consistent sixth-order convergence as anticipated.
Considering u 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 a 1 , we used For the two end points near a 2 , we used The second approach turns out to be slightly better than the first one in decreasing the L 2 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 r 0 to r d , 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.

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 r 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 r r 0 .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 r r 0 , where r 0 is close to zero and 1/r 2 0 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 Advances in Mathematical Physics 15 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 D 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 r i and the end points in the D layer, the combined method uses 3.10 and 3.11 .For the interior points in both the inner core and the D 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 r 0 , we develop the following numerical schemes based on the FIM.First, by introducing ∂φ/∂r u, 3.4 can be reduced into two first order differential equations: To approximate the first-order derivative near r r 0 , a 3-point stencil difference approximation based on the FIM is 3.15 We define the radii at the grid points i, i 1, and i 2 as r i , r i 1 , and r i 2 , respectively.Then we expand every term at the point i 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, r i 1 − r i is smaller than r i 2 − r i 1 , the truncation error will actually be smaller than the corresponding error in a uniform grid at the same number of grid points.Therefore a fourthorder 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

Advances in Mathematical Physics
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 η 0 so that 3.15 becomes

3.19
For point 2, we choose ξ 0, and 3.15 becomes where the coefficients are

3.21
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 ∂φ/∂r.It is observed from the left plot that the dipole component L 2 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 L becomes larger, the magnitude of φ decreases rapidly.For example, at L 10, φ reduces to less than 10% of the value at L 2. When L 28, this magnitude drops to less than 1%.There is a saddle point at the inner-outer core boundary near r 0.34 where the derivative decreases on the left side and increases on the other for all the L values.Therefore ∂φ/∂r is not smooth and ∂ 2 φ/∂r 2 is discontinuous at this point.This is also why we solve the first derivative directly instead of the second derivative in the formulation.
Figure 7 shows the local error distribution across the entire region for L 2 from the results with the JTr and the combined method.The state variable φ is on the left and the derivative of the state variable ∂φ/∂r is on the right.The resolution was set at 100 grid points.The magnitude of the error for both φ and ∂φ/∂r 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 r 0.34.Errors at the left boundary within the inner core and at the right boundary in the D 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 ∂φ/∂r is at the left boundary in the inner core.The JTr produces a slightly higher error than the combined method.
Figure 8 demonstrates the local error distribution for L 28 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 ∂φ/∂r b is on the right.The resolution was set at 190 grid points.The leading errors for both φ and ∂φ/∂r have shifted to the location around the inflection point.Additionally, the leading errors increased slightly with increased resolution.This is because as L increases by one unit, the coefficient of φ, L L 1 /r 2 , which is called the stiffness of 3.4 , almost doubles its value.As r 2 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 r 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.
As L 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 L 2. Since the value of L is small, the stiffness of 3.4 is very small and both the JTr and the combined method provide excellent convergence rates of above 6, 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.
As L becomes even larger, the stiffness of 3.4 becomes more than the square of L. 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 L 25, the combined method has a convergence rate 6.05 for φ and 6.01 for ∂φ/∂r.But when L 28, 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   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 D 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 finitedifference 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:

A.4
Equation A.4 is a 4 × 4 system for these four coefficients.The explicit expressions for these coefficients are given as A.6

Figure 1 :
Figure 1: a Solutions of the benchmark problem 1. b Solution errors of the sixth-and fourth-order compact schemes.The total grid points were N 100 and the parameter was 0.01.

Figure 2 :Figure 3 :
Figure 2: Comparisons of convergence rates of the fourth-order compact schemes for problem 1 at different values of : solution error φ a ; derivative error ∂φ/∂r b .

Figure 4 :
Figure 4: a Solutions of the benchmark problem 2. b The solution error of the sixth-order compact scheme.The resolution was set at N 100 and the parameter was 0.01.

Figure 5 :
Figure 5: Comparisons of convergence rates of the sixth-order compact schemes for the problem 2 at different values of : solution error φ a ; derivative error ∂φ/∂r b .

Figure 6 :
Figure 6: Comparisons of the solution φ a and the derivative of the solution ∂φ/∂r b for the sample geodynamo problem.Results were obtained from the sixth-order compact schemes with the Jacobian transformation at different values of L.

Figure 7 :
Figure 7: Comparisons of the error in solution φ a and the derivative of the solution ∂φ/∂r bfor the sample geodynamo problem.Results were obtained separately from the sixth-order compact schemes with the Jacobian transformation and the combined method at L 2. The total number of grid points in three subdomains is 100.

Figure 8 :
Figure 8: Comparisons of the error in solution φ a and in the derivative of the solution ∂φ/∂r b for the sample geodynamo problem.Results were obtained separately from the sixth-order compact schemes with the Jacobian transformation at L 28 on a nonuniform grid with 190 points.

Figure 9 :
Figure 9: Comparisons of convergence rates of the sixth-order compact schemes for the sample geodynamo problem.Results were obtained separately from the Jacobian transformation and the combined method at L 2: solution error φ a ; derivative error ∂φ/∂r b .

Table 1 :
A summary of the convergence rates for the state variable φ and its derivative ∂φ/∂r at different values of L.
, φ p r i , φ p 1 , and φ p 1 can be computed from A.3 as