First order hyperbolic approach for Anisotropic Diffusion equation

TL;DR

This paper introduces a high-order finite difference solver based on the first-order hyperbolic system method for anisotropic diffusion problems, achieving uniform accuracy regardless of anisotropy degree.

Contribution

The paper develops a fifth-order hyperbolic method that simplifies constructing accurate schemes independent of anisotropy levels, with extensions to magnetized electron simulations.

Findings

Achieves high accuracy with fifth-order schemes

Demonstrates uniform accuracy regardless of anisotropy degree

Successfully applied to magnetized electron test case

Abstract

In this paper, we present a high order finite difference solver for anisotropic diffusion problems based on the first-order hyperbolic system method. In particular, we demonstrate that the construction of a uniformly accurate fifth-order scheme that is independent of the degree of anisotropy is made straightforward by the hyperbolic method with an optimal length scale. We demonstrate that the gradients are computed simultaneously to the same order of accuracy as that of the solution variable by using weight compact finite difference schemes. Furthermore, the approach is extended to improve further the simulation of the magnetized electrons test case previously discussed in Refs.[J. Comput. Phys., 284 (2015) 59-69 and 374 (2018) 1120-1151]. Numerical results indicate that these schemes are capable of delivering high accuracy and the proposed approach is expected to allow the hyperbolic method to be successfully applied to a wide variety of linear and nonlinear problems with anisotropic diffusion.

Tables

Table 1. Table 1: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for ν o ​ p ​ t subscript 𝜈 𝑜 𝑝 𝑡 \nu_{opt} and ν = 1 𝜈 1 \nu=1 by 3rd order explicit, 5th order explicit and compact schemes for anisotropic diffusion problem given in Example 4.1 .

	${\nu }_{opt}$						$\nu =1$
Number	Upwind-3E		Upwind-5E		Upwind-5C		Upwind-5C
of points	Error	Order	Error	Order	Error	Order	Error	Order
$8^2$	4.95E-05		2.81E-06		3.36E-06		3.82E-03
${16}^2$	2.60E-06	4.25	1.12E-07	4.65	1.20E-07	4.80	1.07E-04	5.16
${32}^2$	3.59E-07	2.86	6.15E-09	4.19	6.41E-09	4.23	2.30E-06	5.54
${64}^2$	5.12E-08	2.81	2.27E-10	4.76	2.40E-10	4.74	5.64E-08	5.35
${128}^2$	6.84E-09	2.90	7.57E-12	4.91	8.10E-12	4.89	1.54E-09	5.20

Table 2. Table 2: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for ν o ​ p ​ t subscript 𝜈 𝑜 𝑝 𝑡 \nu_{opt} by 3rd order explicit, 5th order explicit and compact schemes for anisotropic diffusion problem given in Example 4.2a .

Number	Upwind-3E		Upwind-5C		Upwind-5E
of points	error	order	error	order	error	order
${16}^2$	3.61E-03	-	2.75E-04	-	1.00E-03	-
${32}^2$	5.13E-04	2.82	7.15E-06	5.26	3.77E-05	4.73
${64}^2$	6.71E-05	2.93	2.09E-07	5.09	1.22E-06	4.95
${128}^2$	8.44E-06	2.99	6.39E-09	5.04	3.81E-08	5.00
${256}^2$	1.05E-06	3.00	1.98E-10	5.01	1.19E-09	5.01

Table 3. Table 3: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for D | | = 10 9 D_{||}=10^{9} , Example 4.2b .

Number	Upwind-3E		Upwind-5E		Upwind-5C
of points	error	order	error	order	error	order
$8^2$	5.81E-02		3.03E-02		3.18E-02
${16}^2$	6.90E-03	2.13	6.03E-03	2.33	5.16E-03	2.62
${32}^2$	1.78E-03	1.95	1.34E-03	2.17	1.24E-03	2.06
${64}^2$	3.04E-04	2.55	2.10E-04	2.68	1.92E-04	2.69
${128}^2$	4.11E-05	2.89	2.64E-05	2.99	2.06E-05	3.22

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for D | | = 10 9 D_{||}=10^{9} , Example 4.3b .

Number	Upwind-3E		Upwind-5C		Upwind-5E
of points	error	order	error	order	error	order
${16}^2$	4.00E-04		1.21E-04		1.80E-04
${32}^2$	5.39E-05	2.89	1.52E-05	3.00	1.86E-05	3.27
${64}^2$	6.94E-06	2.96	1.89E-06	3.01	2.00E-06	3.22
${128}^2$	8.78E-07	2.98	2.12E-07	3.15	2.24E-07	3.15
${256}^2$	1.10E-07	2.99	2.49E-08	3.09	2.75E-08	3.03

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for primary variable T by 3rd order explicit, 5th order explicit and compact schemes for anisotropic diffusion problem given in Example 4.4a .

Number	Upwind-3E		Upwind-5C		Upwind-5E
of points	error	order	error	order	error	order
$D_{\left|\right|}={10}^0$
${16}^2$	1.69E-04	-	4.30E-06	-	4.73E-06	-
${32}^2$	2.29E-05	2.89	1.50E-07	4.84	1.68E-07	4.81
${64}^2$	2.97E-06	2.95	4.73E-09	4.98	5.47E-09	4.94
${128}^2$	3.77E-07	2.98	1.47E-10	5.01	1.73E-10	4.98
${256}^2$	4.74E-08	2.99	4.56E-12	5.01	5.42E-12	4.99
$D_{\left|\right|}={10}^9$
${16}^2$	8.94E-05		3.32E-06		3.54E-06
${32}^2$	1.23E-05	2.87	1.27E-07	4.71	1.35E-07	4.71
${64}^2$	1.62E-06	2.92	4.16E-09	4.93	4.54E-09	4.90
${128}^2$	2.08E-07	2.96	1.31E-10	4.99	1.45E-10	4.97
${256}^2$	2.64E-08	2.98	4.10E-12	5.00	4.60E-12	4.98

Table 6. Table 6: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for D | | = 10 9 D_{||}=10^{9} , Example 4.4b .

Number	Upwind-3E		Upwind-5C		Upwind-5E
of points	error	order	error	order	error	order
${16}^2$	3.07E-03	-	5.99E-04	-	8.87E-04	-
${32}^2$	3.80E-04	3.01	1.18E-05	5.67	2.83E-05	4.97
${64}^2$	4.73E-05	3.00	2.47E-07	5.58	9.08E-07	4.96
${128}^2$	5.86E-06	3.01	5.89E-09	5.39	2.79E-08	5.02
${256}^2$	7.28E-07	3.01	1.65E-10	5.16	8.69E-10	5.00

Equations

$$\frac{\partial T}{\partial t}=-\nabla .\text{q}+\text{S},\text{q}=-(\text{D}.\nabla T),$$

$$\displaystyle\quad\mathbf{D}=\left[\begin{array}[]{ccc}D_{xx}&D_{xy}\\ \\ D_{yx}&D_{yy}\\ \end{array}\right]=\left[\begin{matrix}D_{||}\cos^{2}{\left(\beta\right)}+D_{\perp}\sin^{2}{\left(\beta\right)}&\frac{1}{2}\left(D_{||}-D_{\perp}\right)\sin{\left(2\beta\right)}\\ \\ \frac{1}{2}\left(D_{||}-D_{\perp}\right)\sin{\left(2\beta\right)}&D_{||}\sin^{2}{\left(\beta\right)}+D_{\perp}\cos^{2}{\left(\beta\right)}\end{matrix}\right],$$

$$0$$

$$=\frac{\partial }{\partial x}\left(D_{xx}\frac{\partial T}{\partial x}+D_{xy}\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial y}\left(D_{yx}\frac{\partial T}{\partial x}+D_{yy}\frac{\partial T}{\partial y}\right)+S.$$

$$\frac{\partial T}{\partial \tau }-\frac{\partial }{\partial x}\left(D_{xx}g+D_{xy}h\right)-\frac{\partial }{\partial y}\left(D_{xy}g+D_{yy}h\right)$$

$$\frac{\partial g}{\partial \tau }-\frac{\partial T}{\partial x}$$

$$\frac{\partial h}{\partial \tau }-\frac{\partial T}{\partial y}$$

$${\mathbf{P}}^{-1}\frac{\partial Q}{\partial \tau }+\frac{\partial {\mathbf{E}}_x}{\partial x}+\frac{\partial {\mathbf{E}}_y}{\partial y}=\mathbf{S},$$

$$\displaystyle\mathbf{Q}=\left[\begin{array}[]{ccc}T\\ g\\ h\end{array}\right],\mathbf{E}_{x}=\left[\begin{array}[]{ccc}-D_{xx}g-D_{xy}h\\ -T\\ 0\end{array}\right],$$

$${\mathbf{P}}^{-1}=\left[\begin{array}{ccc}1& 0& 0\\ 0& T_r& 0\\ 0& 0& T_r\end{array}\right],$$

$${\mathbf{PA}}_n=\mathbf{P}\frac{\partial ({\mathbf{E}}_x{n}_x+{\mathbf{E}}_y{n}_y)}{\partial Q}=\left[\begin{array}{ccc}0& -D_{xx}{n}_x-D_{xy}{n}_y& -D_{xy}{n}_x-D_{yy}{n}_y\\ -\frac{n_x}{T_r}& 0& 0\\ -\frac{n_y}{T_r}& 0& 0\end{array}\right],$$

$${\lambda }_1=0,{\lambda }_{2,3}=\pm \sqrt{\frac{1}{T_r}\left(D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2\right)},$$

$${\mathbf{R}}_n=\left[\begin{array}{ccc}\frac{D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2}{n_y\sqrt{\frac{1}{T_r}\left(D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2\right)}}& 0& -\frac{D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2}{n_y\sqrt{\frac{1}{T_r}\left(D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2\right)}}\\ \frac{n_x}{n_y}& -\frac{D_{xy}{n}_x+D_{yy}{n}_y}{D_{xx}{n}_x+D_{xy}{n}_y}& \frac{n_x}{n_y}\\ 1& 1& 1\end{array}\right].$$

$${\mathbf{\Lambda }}_{\mathbf{n}}=\sqrt{\frac{1}{T_r}\left(D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2\right)}\left[\begin{array}{ccc}-1& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right]$$

$$|{\mathbf{PA}}_n|={\mathbf{R}}_n|{\mathbf{\Lambda }}_n|{\mathbf{R}}_n^{-1}=\lambda \left[\begin{array}{ccc}1& 0& 0\\ 0& \frac{n_x\left(D_{xx}{n}_x+D_{xy}{n}_y\right)}{D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2}& \frac{n_x\left(D_{xy}{n}_x+D_{yy}{n}_y\right)}{D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2}\\ 0& \frac{n_y\left(D_{xx}{n}_x+D_{xy}{n}_y\right)}{D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2}& \frac{n_y\left(D_{xy}{n}_x+D_{yy}{n}_y\right)}{D_{xx}{n}_x^2+2D_{xy}{n}_x{n}_y+D_{yy}{n}_y^2}\end{array}\right],$$

$${\mathbf{Q}}^{\beta }={\mathbf{Q}}_0\exp (i({\beta }_{x}x+{\beta }_{y}y)/h),$$

$$\frac{d{\mathbf{Q}}_0}{d\tau }=\mathbf{M}{Q}_0,$$

$$\mathbf{M}=\left[\begin{matrix}0&\frac{i(D_{xx}\beta_{x}+{\color[rgb]{0,0,0}D_{xy}}\beta_{y})}{h}&\frac{i(D_{xy}\beta_{x}+{\color[rgb]{0,0,0}D_{yy}}\beta_{y})}{h}\\ \\ \frac{i\beta_{x}}{hT_{r}}&-1/T_{r}&0\\ \\ \frac{i\beta_{y}}{hT_{r}}&0&-1/T_{r}\end{matrix}\right].$$

$${\lambda }_{1,2}$$

$${\lambda }_3$$

$$h^2-4T_r\left(D_{xx}{\beta }_x^2+2D_{xy}{\beta }_x{\beta }_y+D_{yy}{\beta }_y^2\right)<0,$$

$$T_r>\frac{h^2}{4\left(D_{xx}{\beta }_x^2+2D_{xy}{\beta }_x{\beta }_y+D_{yy}{\beta }_y^2\right)}.$$

$$T_r=\frac{1}{(2\pi {)}^2\left(D_{xx}+2D_{xy}+D_{yy}\right)}.$$

$$T_r=\frac{L_r^2}{{\nu }_{opt}},\text{where}{\nu }_{opt}=D_{xx}+2D_{xy}+D_{yy},$$

$$L_r=\frac{2h}{\pi /\tilde{N}(\pi /\tilde{N}+4)},\text{where}\tilde{N}=\frac{1}{h\sqrt{(\frac{1}{L_x^2}+\frac{1}{L_y^2})}},$$

$$x_{j+\frac{1}{2}}=x_a+j\Delta x,y_{k+\frac{1}{2}}=y_a+k\Delta y$$

$$\Delta x=\frac{x_b-x_a}{N_x},\Delta y=\frac{y_b-y_a}{N_y}$$

$${\frac{\partial Q}{\partial t}|}_{j,k}+{\frac{\partial E_x}{\partial x}|}_{j,k}+{\frac{\partial E_y}{\partial y}|}_{j,k}={S|}_{j,k}$$

$$\displaystyle\rm{4^{th}\ order:};\Bigl{|}\frac{\partial E_{x}}{\partial x}\Bigr{|}_{j,k}$$

$$\displaystyle\rm{6^{th}\ order:};\Bigl{|}\frac{\partial E_{x}}{\partial x}\Bigr{|}_{j,k}$$

$$+\frac{3}{640}\left(E_{j+\frac{5}{2},k}-E_{j-\frac{5}{2},k}\right)],$$

$$5^{\mathrm{th}}\mathrm{order}\mathrm{compact}\{\begin{array}{l}\frac{1}{2}{Q}_{j-\frac{1}{2}}^L+Q_{j+\frac{1}{2}}^L+\frac{1}{10}{Q}_{j+\frac{3}{2}}^L=\frac{1}{10}{Q}_{j-1}+Q_j+\frac{1}{2}{Q}_{j+1}\\ \frac{1}{10}{Q}_{j-\frac{1}{2}}^R+Q_{j+\frac{1}{2}}^R+\frac{1}{2}{Q}_{j+\frac{3}{2}}^R=\frac{1}{2}{Q}_j+Q_{j+1}+\frac{1}{10}{Q}_{j+2}.\end{array}$$

$$3^{\mathrm{rd}}\mathrm{order}$$

Full text

1 Introduction

Anisotropic diffusion occurs in many physical applications in which the rate of diffusion in a certain direction can be orders of magnitude higher than the other. Thermal conductivity in fusion plasmas [1], image processing [2], biological process, and medical imaging [3] are some of the examples. Diffusion tensors can be extremely anisotropic in magnetized and high-temperature plasma process which poses a challenging problem for computational simulations. Due to extreme anisotropy, the diffusion phenomenon is effectively aligned with the magnetic field lines. Such alignment may lead to the parallel diffusion coefficient being orders of magnitude, up to $10^{9}$ in fusion plasmas, larger than perpendicular diffusion coefficient.

These extreme anisotropies can have strict requirements on the numerical methods used to model the anisotropic diffusion and magnetohydrodynamics equations since any misalignment of the grid can lead to significant numerical diffusion in the perpendicular direction to the magnetic field line. One approach is to solve the equations on a computational grid that is aligned with the applied magnetic field which can automatically take care of the directionality of the diffusion coefficients. This approach, known as the magnetic field-aligned mesh (MFAM), has been successfully used in modeling plasma propulsion devices and fusion plasmas [4, 5]. This approach can have problems in the case of crossing field lines where local non-alignment unavoidable. Also, one has to reconstruct the mesh for the situations involving time-varying field lines, or magnetic field induce flow, and it can be a cumbersome task. It is beneficial to develop numerical methods that are applicable even on non-aligned grids with minimal numerical diffusion.

Gunter et al. [6] has developed “Symmetric” and “asymmetric” scheme which are simple and easy to implement with low perpendicular numerical pollution. Sovinec et al. [7] have used a higher-order finite element method in the direction of the larger diffusion coefficient to reduce the numerical diffusion. On the other hand, Degond et al. have developed Asymptotic preserving schemes (AP) in a series of papers [8, 9] for the strongly anisotropic diffusion equation in which the equation itself is split into two parts, a limit problem for infinite anisotropy and the original singular perturbation problem. Mentrelli and Negulescu [10] further extended the AP schemes to the generalized nonlinear cases, and Chacon et al. [11] applied a more generic semi-Lagrangian approach to unsteady anisotropic diffusion. Van Es et al. [12] has developed a second-order accurate aligned finite difference method for anisotropic diffusion problems. Another important aspect of the anisotropic diffusion equation is the positivity of the temperature or space potential. A numerical scheme should be robust against the development of nonphysical negative temperatures during the simulation and must satisfy positivity and monotonicity. For example, negative temperatures, in Hall thruster modeling can lead to decreased electron currents and negative Joule power density near the cathode region. Positivity has been attained by nonlinear schemes proposed in the literature, imposed via limiters by Kuzmin et al.[13] and Sharma and Hammett [14], but limiters can lead to low-order spatial accuracy and are also limited to moderate anisotropies of the order $\approx 10^{3}$.

Recently an upwind formulation for diffusion equation has been introduced by Nishikawa [15] based on the residual-distribution (RD) method and later with a finite-volume method [16]. The mathematical strategy of this approach is to split the second order partial differential equation into a set of first-order differential equations by adding new variables and pseudo-time advancement terms such that the diffusion equation is reformulated as a hyperbolic system. This radical approach has been shown to offer several advantages over conventional methods, such as accelerated convergence for steady state solution and higher order of accuracy for both primary and gradient variables, as demonstrated for diffusion [16], the incompressible/compressible Navier-Stokes equations [17, 18], third-order dispersion equations [19], an incompressible magnetohydrodynamics model [20], an elliptic distance-function model [21], and so on. The original approach of Nishikawa has been further extended to a constant diffusion tensor by Lou et al. [22, 23] discretized with the reconstructed discontinuous Galerkin scheme (rDG). The hyperbolic approach has also been implemented for an anisotropic diffusion equation, based on high-order finite-volume schemes in Ref.[24] (see Example 4.1.5 in the reference). Even though fifth order accuracy was demonstrated for the test cases considered, it has been known that the method yields significantly large errors for increasing degree of anisotropy. To address this issue, we will extend the analysis in Ref.[25] and derive an optimal length scale that defines a relaxation time in a hyperbolic formulation of anisotropic diffusion equations. As we will demonstrate, the optimal length scale renders the method independent of the degree of anisotropy, thus achieving the same level of errors for a wide range of anisotropy. The method is developed based on a local-preconditioning formulation and extended to variable-coefficient and nonlinear anisotropic diffusion equations.

The rest of the paper is organized as follows. Hyperbolic approach and derivation of the preconditioned formulation for the anisotropic diffusion are described in Section 2. The cell-centered explicit and compact upwind finite difference upwind schemes and the implementation of boundary conditions are presented in Section 3. Several test cases in one and two-dimensional problems for anisotropic diffusion equation are presented in Section 4. These numerical experiments validate our numerical scheme and corroborate the high-order accuracy and implementation of boundary conditions. Finally, Section 5 summarizes our conclusions.

2 Hyperbolic formulation for anisotropic diffusion equation

In this section, the hyperbolic approach for anisotropic diffusion equation is briefly explained. Anisotropic thermal diffusion equation can be described as

[TABLE]

where ${T}$, $\mathbf{q}$, $\mathbf{S}$ and $\mathbf{D}$ represent temperature, heat flux, source terms and the diffusion tensor respectively. For a two-dimensional problem the diffusion tensor is given by

[TABLE]

where $\mathbf{\beta}$ is the angle between the grid and the magnetic field lines and $D_{\perp}$ and $D_{||}$ are the perpendicular and parallel diffusion coefficients, shown in Fig. 1, respectively. The diffusion tensor $\mathbf{D}$ is assumed to be positive definite: ${\bf x}^{t}\mathbf{D}{\bf x}>0$ for an arbitrary non-zero vector ${\bf x}$.

By dropping the time derivative and substituting q in the equation, we get the following steady state anisotropic diffusion equation

[TABLE]

Note that the physical time derivative can be retained for unsteady simulations. This paper focuses on steady problems and the developed steady solver can be employed to solve unsteady residual equations for unsteady simulations. See Refs.[26, 27, 28] for the hyperbolic method applied to unsteady problems.

By introducing new variables $g$ and $h$, where $g=\frac{\partial T}{\partial x}$ and $h=\frac{\partial T}{\partial y}$, and the system is made hyperbolic by adding pseudo-time derivatives with respect to all variables as follows:

[TABLE]

This formulation is proposed by Lou et al. [22, 23]. The Eq.(3) can be written as a preconditioned system as follows,

[TABLE]

where

[TABLE]

and the preconditioning matrix is given by,

[TABLE]

where $T_{r}$ is the relaxation time which will be explained later. For an arbitrary $T_{r}$, this system is equivalent to the anisotropic diffusion equation (2) in the steady state or as soon as the pseudo time terms are dropped. The preconditioned hyperbolic formulation was first introduced in Ref.[29] to extend the hyperbolic method to nonlinear equations and simplify the analysis of the eigenstructure of the Navier-Stokes equations by eliminating the need to differentiate the viscosity given as a function of a solution. It is also deliberately constructed as a preconditioned conservative form, so that any discretization method designed for a conservation law can be directly applied [29]. For the anisotropic diffusion equation, a similar formulation requires us to define four extra variables corresponding to the four diffusive flux components: $D_{xx}g$, $D_{xy}h$, $D_{xy}g$ and $D_{yy}h$. To minimize the number of equations, we introduce, instead, the derivatives of the temperature as extra variables. As a result, the diffusion coefficients remain in the flux, and thus they need to be differentiated in deriving the flux Jacobian. However, as pointed out in Ref.[30], to construct a numerical dissipation matrix, the flux Jacobian does not have to be exactly derived, and the coefficients can be frozen and thus not needed to be differentiated. This type of formulation is very useful for developing high-order methods [31, 32]. Note that the preconditioned formulation is necessary, for whatever variables are chosen, to apply the hyperbolic method to variable-coefficient and nonlinear equations because $T_{r}$ is not a constant and thus cannot be included in the flux vectors [30]. To construct an upwind flux, consider the preconditioned Jacobian of the flux projected along an arbitrary vector, $n=(n_{x},n_{y})$:

[TABLE]

where the diffusion coefficients have been assumed to be constant for the sake of derivation (but they can be functions of space or solutions in the final form). This system is hyperbolic since $\mathbf{PA}_{n}$ has real eigenvalues, which are given by

[TABLE]

and linearly independent right-eigenvectors:

[TABLE]

The absolute Jacobian $|\mathbf{PA}_{n}|$, which is required for current upwind flux approach, is constructed by right-eigenvector matrix $\mathbf{R}_{n}$ and the diagonal eigenvalue-matrix, $\mathbf{\Lambda_{n}}$,

[TABLE]

as follows:

[TABLE]

where $\lambda$=${\sqrt{\frac{1}{T_{r}}\left(D_{xx}n_{x}^{2}+2D_{xy}n_{x}n_{y}+D_{yy}n_{y}^{2}\right)}}$.

2.1 Fourier analysis for optimum relaxation time and length scale

Due to the significant difference between the eigenvalues, $10^{9}$ and 1, the condition number (ratio of the maximum eigenvalue to the minimum eigenvalue) of the system increases and degrades the convergence performance and also the accuracy of the numerical result. The optimal choice of relaxation time improves the condition number of the problem and thereby improves the accuracy such that the results are independent of the degree of anisotropy. In this subsection, we derive optimal length scale and relaxation time such that the error remains the same regardless of the degree of anisotropy. The importance of length scale and relaxation time for the isotropic diffusion are discussed in Ref. [33]. Optimal values for $T_{r}$ is derived by following the procedure briefly described in [15, 33] for a two-dimensional hyperbolic diffusion system. In short, we substitute a Fourier mode into the first-order version of the finite-difference or finite-volume scheme in the semi-discrete form (on a regular grid) to derive a pseudo-time evolution of the amplitude of the Fourier mode, and derive the optimal values by requiring that the eigenvalues of the evolution matrix (the RHS of the ODE for the Fourier mode amplitude) become complex, so that the Fourier mode will propagate rather than purely get damped. If the Fourier mode propagates, the scheme has a property of removing errors by propagation in addition to error damping, which is much faster than pure damping. The slow convergence of relaxation schemes for elliptic problems comes from the pure damping, and the hyperbolic method improves the convergence by introducing error propagation into the scheme.

Consider the preconditioned form of the Eq.(4). To derive an optimal $T_{r}$, we ignore the physical time derivative $\frac{\partial Q}{\partial t}$ and focus on the pseudo-time system. Our goal is to define $T_{r}$ such that all possible Fourier modes will propagate by the system in the pseudo time. Define the Fourier mode as

[TABLE]

on a uniform grid of spacing h, where $i$ = $\sqrt{-}1$ and Q_{0}=({\color[rgb]{0,0,0}T_{0}},g_{0},h_{0}) is the vector of amplitudes. Then, substitute it into the first-order system to obtain

[TABLE]

where

[TABLE]

The eigenvalues of $\mathbf{M}$ are given by

[TABLE]

The third eigenvalue corresponds to the inconsistency-damping mode associated with the source terms [15]. We focus on the remaining pair and determine $T_{r}$ to make them complex conjugate, so that the Fourier mode propagates. That is, we require z $<$ 0 for all possible frequencies:

[TABLE]

which, since $\mathbf{D}$ is positive definite, leads to

[TABLE]

The lower bound reaches the maximum for the smoothest discrete mode: $\beta_{x}$ = $\beta_{y}$ = $\pi$h, and therefore the inequality is satisfied for all possible frequencies if we set

[TABLE]

This implies that we define $T_{r}$ as

[TABLE]

where ${L}_{r}=1/(2\pi)$ is a relaxation length. In the hyperbolic approach proposed by Lou et al. [22], based on reconstructed Discontinuous Galerkin scheme (rDG), $\nu_{opt}$ is taken as 1. It will be shown here that such a hyperbolic diffusion scheme can be high order accurate, but it is not independent of the diffusion coefficient and loses accuracy with an increasing degree of anisotropy. It has been shown that the relaxation length associated with the hyperbolic diffusion system is a dimensional quantity and must be properly scaled by a reference length for a given problem to guarantee dimensional-invariance of numerical solutions [25]. For each reference length chosen, an optimal value of the corresponding non-dimensionalized relaxation length needs to be found to achieve optimal performance of a hyperbolic solver, called an optimal reference length. The optimal value is given by $\frac{1}{2\pi}$ for a unit square domain, but it is not necessarily optimal for a general domain. A practical approach for obtaining a reference length is derived in Ref. [25]. For a rectangular domain, (x, y) $\in$ $[0,L_{x}]\times[0,L_{y}]$, the following equation can be used,

[TABLE]

where $h$ is grid spacing for a given uniform grid but has a very minor impact on $L_{r}$ when $h\ll 1$. This approach is considered in this current paper and is found important for magnetized electrons test case.

3 Review of numerical schemes and boundary conditions

In this section the procedures of cell centered upwind finite difference schemes [34, 35] and boundary conditions [24] are briefly explained. The domain is discretized uniformly into a Cartesian grid of size $N_{x}\times N_{y}$, where $N_{x}$ and $N_{y}$ are the numbers of cells in $x$-and $y$-coordinate directions, and the domain is defined as $I_{j,k}=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right]\times\left[y_{k-\frac{1}{2}},y_{k+\frac{1}{2}}\right]$ for $1\leq j\leq N_{x}$, $1\leq k\leq N_{y}$, where

[TABLE]

and

[TABLE]

The Eq.(4) in semi-discrete form by finite difference scheme can be written as

[TABLE]

where \frac{\partial E_{x}}{\partial x}\big{|}_{j,k}, \frac{\partial E_{y}}{\partial y}\big{|}_{j,k}, and $S|_{j,k}$ are the approximations of the flux derivatives in the $x$ and $y$ directions, and source term at cell centers $(x_{j},y_{k})$, shown in Fig. 2, respectively.

3.1 Overview of upwind finite difference schemes

For the evaluation of the flux derivatives, Lele has developed a family of compact schemes in [36]. Through a detailed Fourier analysis, it has been shown that compact schemes on a cell-centered mesh, staggered, provide better resolution in comparison with the node-centered, collocated, mesh. The cell-centered method involves interpolation of variables to cell interfaces, which can be carried out either by a centered stencil or an upwind biased stencil. Nagarajan et al. [37] and Boersma [38] used cell-centered compact schemes with a centered stencil to simulate compressible flows, and their numerical tests indicate that their methods are quite robust. Unfortunately, central interpolations are known to cause oscillations in the presence of discontinuities and are not suitable for the current hyperbolic approach. On the other hand, Deng and Zhang [34] developed a class of weighted nonlinear compact schemes (WCNS) based on the idea of Weighted Essentially Non-oscillatory approach (WENO). In this approach, the interpolation of variables is carried out through a nonlinear combination of several lower-order stencils, which are upwind biased, similar to the methodology of the WENO scheme [39]. These schemes capture the discontinuities without spurious oscillations and recover formal order of accuracy in smooth regions. Deng et al. [40] has also developed a family of hybrid cell-interface and cell-node compact linear schemes. In this paper, both linear and nonlinear weighted compact schemes are used for spatial discretization. It should be noted that the finite-difference schemes considered here are more efficient and cheaper for the nonlinear problems than the high-order finite volume methods (FVM) which require additional flux computations at quadrature points [41]. To achieve higher order accuracy, third or more, in FVM, the integrals of the fluxes are discretized utilizing a high-order Gaussian quadrature with suitable Gaussian integration points over the faces of the control volume. Due to these difficulties, the higher order finite volume approach is at least 3 to 12 times expensive compared to the finite difference methods in multidimensional problems. Nonomura et al. have discussed the advantages of these schemes over high-order finite volume schemes for multi-component flows in detail(readers can refer [42] for further information). The WCNS consists of three parts: high-order flux difference for flux derivatives, high-order interpolation for cell-interface values, and upwind flux evaluation at the cell interface.

In the first part, a high-order flux difference, cell-interface-to-node, scheme is adapted to calculate the flux derivatives, which are

[TABLE]

where $E_{j+\frac{1}{2},k}$ are fluxes approximated at cell interfaces by linear or nonlinear interpolations. The derivatives in the $y$ direction are computed similarly. The flux derivatives can also be computed compact difference schemes, but Deng and Zhang [34] report that the resolution of the explicit difference scheme is almost same as that of the compact difference scheme and also the explicit differencing is computationally more efficient than the compact scheme. Nonomura et al.[35] proposed a family of cell interface-and-node-to-node differencing schemes for the computation of the first derivatives which are robust than Eq.(24) but in our computations, we found no difference between the results.

In the second part, high-order linear or nonlinear interpolations are used for the computation of $E_{j+\frac{1}{2}}$ in Eq.(24) which requires knowledge of the solutions $Q^{R}$ and $Q^{L}$, shown in Fig. 2. The solutions of $Q^{R}$ and $Q^{L}$ are obtained by upwind interpolation. The ${Q}_{j+\frac{1}{2}}^{(L)}$ is biased to left and similarly ${Q}_{j+\frac{1}{2}}^{(R)}$ is biased to the right for correct upwinding. The interpolation can be either explicit or implicit in space, also known as compact schemes. Details of the implementation of nonlinear interpolation schemes (WCNS) [34, 43] are given in Appendix A. Compact linear interpolation polynomials, denoted as U-5C, for the left and right interface are as follows,

[TABLE]

Polynomials of accuracy up to $9^{th}$ order are considered for the computations. Explicit third and fifth order linear interpolation formulas, denoted as U-3E and U-5E, are given below

[TABLE]

Finally, in the third part, the upwind flux will be computed by an appropriate Riemann solver. One of the advantages of WCNS is that the interpolation procedure is independent of the numerical flux calculations, unlike the WENO interpolation. The numerical flux can be evaluated by both flux vector splitting and flux difference-splitting schemes (see [34] for details). We evaluate the numerical flux by the following upwind flux:

[TABLE]

where $|\mathbf{A}|$ is the flux Jacobian, given in Eq.(11). In the local preconditioning method [44], the numerical flux $\mathbf{E}_{j+\frac{1}{2},k}$ is constructed based on the preconditioned Jacobian. $\mathbf{P}^{-1}$ is multiplied in the dissipation term in the above expression to cancel the effect of $\mathbf{P}$ that will be multiplied in the final residual which is a standard procedure in the construction of upwind fluxes for a preconditioned PDE. The complete procedure can be summarized as follows :

Interpolate the conservative variables from cell-node to cell interface, Eq.(25) or (26). 2. 2.

compute the upwind flux evaluation at the cell interface, Eq.(27), and 3. 3.

evaluate the derivative by cell-interface-to-node differencing, Eq.(24). Fourth order differencing formula is used for third order interpolation, and sixth order differencing formula is used for fifth order explicit and implicit interpolation formulae.

3.2 Time discretization

After discretizing the spatial derivative, we form a set of ordinary differential equations: e.g., at a cell $(j,k)$,

[TABLE]

where

[TABLE]

and integrate it in the pseudo time by the following third order TVD Runge-Kutta method, to drive the solution to the steady state:

[TABLE]

where $\Delta\tau$ is the pseudo time step given by the following equation

[TABLE]

Time integration is performed with a CFL $=$ 0.2 for all the schemes until the residual is reduced by ten orders of magnitude in the L1 norm.

3.3 Boundary conditions

In this section, the implementation of boundary conditions is described. In the hyperbolic approach the Neumann boundary conditions are also implemented as Dirichlet boundary condition through the gradient variables, and therefore we describe only Dirichlet boundary conditions, more details are given in Ref. [24]

**Primary variable $T$

**The physical domain is extended through ghost cells and the Dirichlet boundary condition is employed at the cell interface, say $T_{\frac{1}{2}}$. The $3^{rd}$ and $5^{th}$ order accurate polynomials for the left and right boundaries are given by Eqs. (32) and (33).

[TABLE]

[TABLE] 2. 2.

**Gradient variables $g$ and $h$

**In isotropic diffusion problems boundary conditions are required only for one of the gradient variables, i.e. $g$ in x-direction and $h$ in y-direction, but in anisotropic diffusion boundary values should be specified for both the gradient variables, $g$ and $h$. The $3^{rd}$ and $5^{th}$ order cell-center to cell-center extrapolation formula are given by Eq.(34) and 35,

[TABLE]

[TABLE]

4 Numerical tests

In this section, various test cases for anisotropic diffusion equation are solved by the upwind schemes and the boundary conditions described earlier. All the simulations are carried until the residuals are dropped below $10^{-10}$ in $L_{1}$ norm by the TVD-RK time integration. Uniform cartesian meshes are used for all the problems. For problems with constant diffusion tensor, $D_{\perp}$ is taken as 1.

4.1 Grid aligned test case

Example 4.1. The following anisotropic diffusion equation identified in Ref. [13], with a spatial domain $(x,y)\in[0,1]\times[0,1]$, has been considered.

[TABLE]

where the diffusion tensor and source term are

[TABLE]

where $\gamma$ is varied from 0 to 9 and the exact solution for this test case is given by,

[TABLE]

which is also imposed as Dirichlet boundary conditions. In this test case, the diffusion coefficients $D_{xy}$ and $D_{yx}$ are zero. Therefore this test case may be considered as a grid-aligned situation where angle $\beta=0^{\circ}$, shown in Fig. 1. Grid refinement study is carried out for all the upwind schemes for two different relaxation times $\nu=1$ and $\nu_{opt}$. It can be seen from Fig. 3 and Table 1 that the design order of accuracy is obtained for both the relaxation times for the primary variable, $T$, for all the schemes, but the $L_{2}$ error obtained by $\nu_{opt}$ is much lower than $\nu=1$. Fig. 3 shows that the error remains the same for increasing value of $D_{||}$, which is varied from $10^{0}$ to $10^{9}$, for $\nu_{opt}$ whereas it increases for $\nu=1$. It is also important to point out that the convergence speed of the computation also improves significantly and is independent of the anisotropy by including the diffusion coefficient in $T_{r}$. It will be shown in further numerical tests that the choice of $\nu_{opt}$ can make the scheme independent of anisotropy and also the angle of misalignment. The solution contours obtained by U-5E on a 32 $\times$ 32 grid and computed values at the geometric center for various schemes are shown in Figs. 4 and 4 respectively.

4.2 Constant angle of misalignment

In the following test cases the applicability of the approach to full diffusion tensor, with a constant angle of misalignment $\beta$ shown in Eq.(1), is analyzed. For all problems in this subsection, we set $D_{\perp}=1$ and $D_{||}=10^{\gamma}$, where $\gamma$ varies again 0 to 9.

Example 4.2a. Second, we consider a simple steady diffusion problem given in Ref. [12] with the following exact solution, which simulates a temperature peak,

[TABLE]

where $\omega=10$ and the angle of misalignment $\beta$ is set to a constant value of $30^{\circ}$. The source term is calculated by plugging Eq.(39) into Eq. (1). The simulations are conducted with grid refinements from 8$\times$8 to 128$\times$128. Figure 5 shows the results obtained with $\nu_{opt}$ that the design order accuracy is obtained for all the schemes and fifth order compact schemes are significantly more accurate in comparison with explicit schemes. As expected, the relaxation time with $\nu_{opt}$ makes the schemes independent of the anisotropy. Figures 5 and 5 show the solution contours and order of accuracy for an anisotropy of $10^{9}$ at a constant $\beta$ of $30^{\circ}$ respectively. All the schemes can capture the temperature peak without any numerical diffusion.

It can be easily confirmed from the Figs. 6(a), 6(b), and 6(c) that the schemes are independent of both angle and degree of anisotropy even for an extreme an anisotropy of $10^{9}$ as the angle of misalignment, $\beta$, is varied from ${0^{\circ}}$ to ${90^{\circ}}$. It can be seen from Fig. 6 that the results obtained by $\nu_{opt}$ are independent of anisotropy whereas the error increases with anisotropy for $\nu=1$.

Example 4.2b. In this test case, the following temperature distribution from Ref. [45] is considered to test the effect of steep temperature gradients,

[TABLE]

where $D_{||}=1,10^{3},10^{6}$, and $10^{9}$ and $\omega$ is taken as 6, and the angle of misalignment is set to a constant value of $45^{\circ}$ for the anisotropy of $10^{9}$. Finite-volume schemes proposed in Ref. [45] were unable to capture the sharp gradients in temperature and were not independent of anisotropy. Only third-order boundary conditions are considered for all the schemes due to steep gradients close to the boundary. It can be seen from Fig. 7 that the hyperbolic scheme can also handle situations with severe anisotropy as well as sharp gradients without losing accuracy. Table 3 show that the third order of accuracy is obtained for the primary variable, $T$, for all the schemes.

4.3 Variable diffusion tensor

In the earlier test cases the angle of misalignment is constant, and thereby the diffusion tensor is constant in space. It is found from the earlier test cases the $T_{r}$ with $\nu_{opt}$ is required to make the schemes independent of the diffusion tensor. In the following test cases, we further test the approach where diffusion coefficients are space varying which can also be thought of as varying angle of misalignment. The physical domain is [0, 1] $\times$ [0, 1] for both the test cases. The preconditioned formulation employed in this study is necessary for variable diffusion tensor problems [30], where $\nu_{opt}$ is a function of space. Again, we set $D_{\perp}=1$ for all problems considered below.

Example 4.3a For the first test case, the exact solution is given by $T(x,y)=1-\tanh\left(\frac{(x-0.5)^{2}+(y-0.5)^{2}}{0.01}\right)$ and the diffusion tensor given by Lou et al. [22, 23] is modified by including extreme anisotropy through $D_{||}$,

[TABLE]

where $\gamma=9$. Again, the grid refinement study is carried out for all three upwind schemes, and the design order of accuracy is achieved for all them, as shown in Fig. 8(b). Once again, the compact scheme is more accurate than explicit interpolation. Solution contours computed on a grid size of 96$\times$96 is shown in Fig. 8(a).

Example 4.3b In the second test case, the exact solution is

[TABLE]

where a = 2 and b =5. The diffusion is tensor is varied by the angle $\beta$ as follows,

[TABLE]

where $D_{||}=10^{9}$ and $\mathbf{\beta}=\arctan\left(x+y\right)$. Again, grid refinement study is carried out for all three upwind schemes and the design order of accuracy is achieved for all them, as shown in Fig. 9. All the schemes are 3rd order accurate because of the boundary conditions as shown in Table 4. Solution contours computed on a grid size of 64$\times$64 is shown in Fig. 9.

4.4 Nonlinear examples

Until now we have considered only linear problems where the diffusion tensor is either constant or varying in space. In real-world applications of MHD or Hall thruster simulations, the diffusion coefficient can be a function of the temperature. In this section, we demonstrate the applicability of the approach for nonlinear anisotropic diffusion problems. Note that the preconditioned formulation needs to be employed for nonlinear problems [30], where $\nu_{opt}$ depends on the solution value.

Example 4.4a First test case is a steady state nonlinear diffusion problem in a unit square, i.e.,

[TABLE]

where the diffusion coefficient, the exact solution, and the source term given by

[TABLE]

The problem reduces to isotropic diffusion, which is considered in Ref. [25], for $D_{||}=1$. Table 5 shows the $L_{2}$ errors for $D_{||}$ = $10^{0}$ and $10^{9}$. It can be seen that the scheme are independent of the degree of anisotropy and design order of accuracy is obtained for all them.

Example 4.4b In the final test case we consider the following full diffusion tensor,

[TABLE]

and $g$ and $h$ are temperature gradients $\frac{\partial T}{\partial x}$ and $\frac{\partial T}{\partial y}$ respectively. Such temperature dependence is proposed by Lafleur et al. [46] for magnetized plasmas in Hall thrusters. The exact solution for this problem is the same as that of the Example 4.2a. Table 6 shows the $L_{2}$ errors and order of convergence, and we can observe that the design order of accuracy is obtained for all the numerical schemes. There is no qualitative difference between the nonlinear problem and the results obtained in Example 4.2a.

4.5 Magnetized electron test case:

Final numerical problem is the magnetized electron test case, which has anisotropic diffusion, considered in [47, 24]. The hyperbolized mass and momentum equations, readers can refer the aforementioned references for further details, are given by

[TABLE]

where $n_{\rm e}$, \vec{u}_{\rm e}{\color[rgb]{0,0,0}=(u_{x},u_{y})}, $\phi$, $\nu_{\rm{ion}}$, and $T_{\rm e}$, are the electron number density, electron velocity, space potential, ionization collision frequency, and electron temperature respectively. The electron mobility tensor $\left[\mu\right]$, similar to $[D]$ in Eq.(1), can be expressed as follows,

[TABLE]

For simplified analysis, the following values are assumed for all the test calculations,

[TABLE]

Then, the system (46) can be written in the form (4) with

[TABLE]

where $(\hat{u}_{x},\hat{u}_{y})=\left[\mu\right]^{-1}\vec{u}_{\rm e}$, and

[TABLE]

Note that the relaxation time $T_{r}$ has been introduced, which allows us to incorporate the dimensionally-consistent optimal length scale, i.e., Equation (19) with Equation (20), without changing steady solutions that we seek. In contrast, the formulation in the previous studies [24, 47] corresponds to the following:

[TABLE]

and thus no consideration was given on the optimal length scale nor the dimensional consistency. Even though the results obtained with the previous formulation are non-oscillatory [24], further studies indicated that the results are not independent of the scale. As we demonstrate below, the optimal definition of $T_{r}$ is very important for achieving anisotropy- and scale-independent solutions. Note also that the above formulation is in the preconditioned conservative form with fluxes $\mathbf{E}_{x}$ and $\mathbf{E}_{y}$, and therefore can be discretized straightforwardly by any method designed for conservative systems, including the method described in Section 3.1. The previous formulation, on the other hand, is in the quasi-linear form with coefficient matrices, and therefore requires the evaluation of the coefficient matrices at each cell center.

It should be noted that the above formulation is slightly different from Equation (4) because the diffusion tensor is not included in the first equation: it gives

[TABLE]

where $\mathbf{P}$ is defined by Equation (50). Fortunately, it leads to the same set of eigenvalues (14) for the matrix ${\bf M}$ in Equation (13) when the Fourier mode is substituted (although the matrix itself is different). Therefore, the relaxation time as derived in Equation (19) is optimal for this system also. The dissipation term is, however, slightly different from Equation (11) because, although the eigenvalues are the same,

[TABLE]

the right-eigenvectors are different:

[TABLE]

which results in the following dissipation matrix:

[TABLE]

We consider the test case calculated by Kawashima et al. [47] with uniformly angled magnetic lines of force at $\theta=45^{\circ}$ from the vertical, shown in Fig. 12. Uniform magnetic confinement ratio, ${\mu}_{||}/{\mu}_{\perp}$ of 1000 is considered throughout the domain. Dirichlet boundary conditions for the non-dimensional space potential are defined at the left and right side boundaries and zero-flux conditions, which are also Dirichlet, are used for the top and bottom boundaries, shown in Eq.(56).

[TABLE]

The simulations are carried out by WCNS-Z scheme, as it is shown in Ref. [24] that linear schemes can cause spurious oscillations. The numerical results computed on a 96 $\times$ 96 grid for a mobility ratio of $10^{3}$. Once again, WCNS-Z scheme along WENO extrapolation (see [24, 48, 49] for further details) are found suitable to reduce the spurious oscillations, as seen in Figs. 13. The numerical results are invariant of the dimensions, and also the positivity of the space potential is also preserved. The domain size of $[0,200]\times[0,100]$ is considered in Refs.[47, 24, 50]. Our scheme yields the same solution independent of the domain size.

5 Conclusions

The important findings of the paper are summarized as follows

Optimal $T_{r}$ and $L_{r}$, preconditioning matrix, are derived through Fourier analysis, and the numerical scheme found to be independent of the degree of anisotropy as well as the angle of misalignment. The preconditioned formulation is suitable for solving variable-coefficient and nonlinear anisotropic diffusion equation. 2. 2.

Consistent third and fifth order accuracies are obtained for anisotropic diffusion by the hyperbolic approach using finite difference schemes. Very high-order compact schemes are found to improve the accuracy for most of the test cases significantly. 3. 3.

Results for magnetized electron test case, with the updated preconditioning matrix along with the optimal length scale formula, are oscillation free for wide-ranging dimensions. It also reiterates the conclusions given in [24] that WENO interpolation approach along with WENO extrapolation for boundary conditions can result in oscillation free results. Dimensionally consistent formulation given in [25] is also required for scale-invariant results.

The extension to curvilinear coordinates using WCNS[51] and field aligned unstructured meshes by discontinuous galerkin methods[23, 52] is currently ongoing and the results will be discussed elsewhere.

Acknowledgements

A.S. is supported by Monbukagakusho fellowship during this work.

Appendix A

In this appendix, the fifth order weighted nonlinear compact scheme is explained. For simplicity, the interpolation polynomials to the left side of the cell interface at $x_{j+\frac{1}{2}}$ are only presented here. The WCNS numerical flux is constructed as follows:

[TABLE]

where

[TABLE]

and $\omega_{k}$ are the nonlinear weights. These nonlinear weights are computed as follows:

[TABLE]

and

[TABLE]

[TABLE]

Here $\epsilon=10^{-6}$ is a small constant to prevent division by zero, and $m$ is equal to 2. In this paper, we considered an improved version of WCNS proposed by Borges et al. in [53]. The improved nonlinear weights are as follows:

[TABLE]

[TABLE]

where the smoothness indicators $\beta_{k}$’s are the same as those given in Equations (61), $\epsilon=10^{-40}$, and $\tau$ is the smoothness indicator of the large stencil. The variable $p$ is used to tune the dispersive and dissipative properties of the scheme. It is reported by Borges et al. [53] that the scheme becomes more dissipative when $p$ is increased. In this paper, denoted as WCNS-Z, $p$= 2 is employed for the magnetized test case for optimal fifth order accuracy.

6 References