First order hyperbolic approach for Anisotropic Diffusion equation
Amareshwara Sainadh Chamarthi Hiroaki Nishikawa Kimiya Komurasaki
physics.comp-ph



























AA_address]Amareshwara Sainadh Chamarthifn1 [cor1]Corresponding author. 
 
E-mail address: s.chamarthi@al.t.u-tokyo.ac.jp (Amareshwara Sainadh  Ch.).
BB_address]Hiroaki Nishikawafn1
AA_address]Kimiya Komurasaki
[AA_address]Department of Aeronautics and Astronautics, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-8656, Japan
[BB_address]National Institute of Aerospace, 100 Exploration Way, Hampton, VA 23666, USA

[fn1]These authors contributed equally to the work.

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.



Hyperbolic system, Anisotropic Diffusion, Compact finite-difference, Nonlinear Diffusion










    


§ 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 <cit.>, image processing <cit.>, biological process, and medical imaging <cit.> 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 <cit.>. 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. <cit.> has developed “Symmetric” and “asymmetric” scheme which are simple and easy to implement with low perpendicular numerical pollution. Sovinec et al. <cit.> 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 <cit.> 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 <cit.> further extended the AP schemes to the generalized nonlinear cases, and Chacon et al. <cit.> applied a more generic semi-Lagrangian approach to unsteady anisotropic diffusion. Van Es et al. <cit.> 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.<cit.> and Sharma and Hammett <cit.>, but limiters can lead to low-order spatial accuracy and are also limited to moderate anisotropies of the order ≈ 10^3. 

Recently an upwind formulation for diffusion equation has been introduced by Nishikawa <cit.> based on the residual-distribution (RD) method and later with a finite-volume method <cit.>. 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 <cit.>, the incompressible/compressible Navier-Stokes equations <cit.>, third-order dispersion equations <cit.>, an incompressible magnetohydrodynamics model <cit.>,  
an elliptic distance-function model <cit.>, and so on. The original approach of Nishikawa has been further extended to a constant diffusion tensor by Lou et al. <cit.> 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.<cit.> (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.<cit.> 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.



§ 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

    ∂ T/∂ t=-∇ . q+ S,   q=-(D.∇ T),

where T, 𝐪, 𝐒 and 𝐃 represent temperature, heat flux, source terms and the diffusion tensor respectively. For a two-dimensional problem the diffusion tensor is given by

    𝐃 = [ [ D_xx D_xy;     ; D_yx D_yy;      ]]  = [D_||cos^2 (β ) + D_⊥sin^2 (β )   1/2(D_|| - D_⊥) sin (2 β )
    
    1/2(D_|| -D_⊥) sin (2 β )    D_||sin^2 (β ) + D_⊥cos^2 (β )],

where β is the angle between the grid and the magnetic field lines and D_⊥ and D_|| are the perpendicular and parallel diffusion coefficients, shown in Fig. <ref>, respectively. The diffusion tensor 𝐃 is assumed to be positive definite: x^t 𝐃 x > 0 for an arbitrary non-zero vector x.















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

    0   =∇ .(D∇ T) + S,
       =∂/∂ x(D_xx∂ T/∂ x+D_xy∂ T/∂ y)+∂/∂ y(D_yx∂ T/∂ x+D_yy∂ T/∂ y) + S.

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.<cit.> for the hyperbolic method applied to unsteady problems.

























By introducing new variables g and h, where g=∂ T/∂ x and h=∂ T/∂ y, and the system is made hyperbolic by adding pseudo-time derivatives with respect to all variables as follows:

    ∂ T/∂τ-   ∂/∂ x(D_xx g+D_xy h) - ∂/∂ y(D_xy g+D_yy h)    = S ,
    ∂ g/∂τ - ∂ T/∂ x   =  -g ,
    ∂ h/∂τ - ∂ T/∂ y   =  -h.


This formulation is proposed by Lou et al. <cit.>. The Eq.(<ref>) can be written as a preconditioned system as follows,
 
    𝐏^-1∂ 𝐐/∂τ
                      +∂𝐄_x/∂x
                      +∂𝐄_y/∂y=𝐒,
 
   where


    𝐐= [ [ T; g; h ]], 
    𝐄_x = [ [ -D_xx g -D_xy h;              -T;               0 ]], 
       𝐄_y = [ [ - D_xy g - D_yy h;                 0;                -T ]],
    𝐒 = [ [ -S; -g; -h ]],
 
    and the preconditioning matrix is given by,
 
    𝐏^-1 = [1     0     0
    0     T_r     0 
     0     0     T_r ],

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 (<ref>) in the steady state or as soon as the pseudo time
terms are dropped. The preconditioned hyperbolic formulation was first introduced in Ref.<cit.>  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 <cit.>.
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.<cit.>, 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 <cit.>. 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 <cit.>. To construct an upwind flux, 
consider the preconditioned Jacobian of the flux projected along an arbitrary vector, n=(n_x, n_y):

    𝐏𝐀_n
    =𝐏∂(𝐄_x n_x+𝐄_y n_y)/∂ Q 
    = [0     - D_xx n_x - D_xy n_y    - D_xy n_x - D_yy n_y
    
    - n_x/T_r    0     0
    - n_y/T_r    0     0],

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 𝐏𝐀_n has real eigenvalues, which are given by

    λ_1 = 0,   λ_2,3 = ±√(1/T_r(D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2)),

 and linearly independent right-eigenvectors:

    𝐑_n= [D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2/n_y√(1/T_r(D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2))    0     - D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2/n_y√(1/T_r(D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2))
    n_x/n_y    - D_xy n_x + D_yy n_y/D_xx n_x + D_xy n_y   n_x/n_y
    1     1     1].

The absolute Jacobian |𝐏𝐀_n|, which is required for current upwind flux approach, is constructed by right-eigenvector matrix 𝐑_n and the diagonal eigenvalue-matrix, Λ_𝐧,

    Λ_𝐧 = √(1/T_r(D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2))[-1     0     0
    0     0     0
    0     0     1]

as follows:

    |𝐏𝐀_n|=𝐑_n|Λ_n| 𝐑^-1_n= λ[1     0     0
    0    n_x(D_xx n_x + D_xy n_y)/D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2   n_x(D_xy n_x + D_yy n_y)/D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2
    0    n_y(D_xx n_x + D_xy n_y)/D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2   n_y(D_xy n_x + D_yy n_y)/D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2],

where λ=√(1/T_r(D_xx n_x^2 + 2 D_xy n_x n_y + D_yy n_y^2)).



 §.§ 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. <cit.>. Optimal values for T_r is derived by following the procedure briefly described in <cit.> 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.(<ref>). To derive an optimal T_r, we ignore the physical time derivative ∂ Q/∂ 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


    𝐐^β = 𝐐_0 exp (i(β_x x +β_y y)/h),

on a uniform grid of spacing h, where i = √(-)1 and Q_0 = (T_0, g_0, h_0) is the vector of amplitudes. Then, substitute it into the first-order system to obtain

    d𝐐_0/dτ = 𝐌Q_0,

where 
 
    𝐌 = [ 0    i(D_xxβ_x+D_xyβ_y)/h   i(D_xyβ_x+D_yyβ_y)/h
    
    iβ_x/hT_r    -1/T_r     0 
    
    iβ_y/hT_r    0     -1/T_r ].

The eigenvalues of 𝐌 are given by



    λ_1,2   =   -  h ±√(z)/2 h T_r ,    z = h^2-4T_r(D_xxβ_x^2 + 2 D_xyβ_xβ_y + D_yyβ_y^2), 
    λ_3   =  - 1/T_r.

The third eigenvalue corresponds to the inconsistency-damping mode associated with the source terms <cit.>. 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:


    h^2-4T_r(D_xxβ_x^2 + 2 D_xyβ_xβ_y + D_yyβ_y^2) < 0,

which, since 𝐃 is positive definite, leads to


    T_r > h^2/4 (D_xxβ_x^2 + 2 D_xyβ_xβ_y + D_yyβ_y^2).

The lower bound reaches the maximum for the smoothest discrete mode: β_x = β_y = πh, and therefore the inequality is satisfied for all possible frequencies if we set


    T_r = 1/(2π)^2(D_xx + 2 D_xy + D_yy).

This implies that we define T_r as
 
    T_r=L_r^2/ν_opt,   where ν_opt = D_xx+2D_xy+D_yy,

where L_r = 1/ (2  π) is a relaxation length. In the hyperbolic approach proposed by Lou et al. <cit.>, based on reconstructed Discontinuous Galerkin scheme (rDG), ν_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 <cit.>. 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 1/2π 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. <cit.>. For a rectangular domain, (x, y) ∈ [0,L_x] × [0,L_y ], the following equation can be used,  
 
    L_r=2h/π/Ñ(π/Ñ+4), where Ñ=1/h√((1/L^2_x+1/L^2_y)),
 
 where h is grid spacing for a given uniform grid but has a very minor impact on L_r when h ≪1.
This approach is considered in this current paper and is found important for magnetized electrons test case.











§ REVIEW OF NUMERICAL SCHEMES AND BOUNDARY CONDITIONS

In this section the procedures of cell centered upwind finite difference schemes <cit.> and boundary conditions <cit.> are briefly explained. The domain is discretized uniformly into a Cartesian grid of size  N_x × 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=[x_j-1/2,x_j+1/2]×[y_k-1/2,y_k+1/2] for 1 ≤ j ≤ N_x, 1 ≤ k ≤ N_y, where

    x_j+1/2 =  x_a + j Δ x,   
      y_k+1/2 = y_a + k Δ y

and

    Δ x = x_b - x_a/N_x,   Δ y = y_b - y_a/N_y

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

    . ∂ Q/∂ t|_j,k + . ∂ E_x/∂ x|_j,k + . ∂ E_y/∂ y|_j,k  = . S |_j,k

where ∂ E_x/∂ x| _j,k, ∂ E_y/∂ y| _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. <ref>, respectively.









 §.§ Overview of upwind finite difference schemes


For the evaluation of the flux derivatives, Lele has developed a family of compact schemes in <cit.>. 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. <cit.> and Boersma <cit.> 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 <cit.> 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 <cit.>. These schemes capture the discontinuities without spurious oscillations and recover formal order of accuracy in smooth regions. Deng et al. <cit.> 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 <cit.>. 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 <cit.> 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


    4^th order:; |∂ E_x/∂ x|_j,k   = 1/Δ x[9/8(E_j + 1/2,k - E_j - 1/2,k) -1/24(E_j + 3/2,k - E_j - 3/2,k) ]
    6^th order:;  |∂ E_x/∂ x|_j,k   = 1/Δ x[
          75/64(E_j+1/2,k - E_j-1/2,k)
          - 25/384( E_j+3/2,k -  E_j-3/2,k) . 
       . + 3/640(E_j+5/2,k - E_j-5/2,k)
        ],

where E_j+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 <cit.> 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.<cit.> proposed a family of cell interface-and-node-to-node differencing schemes for the computation of the first derivatives which are robust than Eq.(<ref>) 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+1/2 in Eq.(<ref>) which requires knowledge of the solutions Q^R and Q^L, shown in Fig. <ref>. The solutions of Q^R and Q^L are obtained by upwind interpolation.   The Q_j+1/2^(L) is biased to left and similarly Q_j+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) <cit.> are given in Appendix A. Compact linear interpolation polynomials, denoted as U-5C, for the left and right interface are as follows,

    5^th order  compact1/2 Q^L_j-1/2+ Q^L_j+1/2 + 1/10 Q^L_j+3/2= 1/10 Q_j-1+  Q_j + 1/2 Q_j+1
    1/10 Q^R_j-1/2+ Q^R_j+1/2 + 1/2 Q^R_j+3/2= 1/2 Q_j+  Q_j+1 + 1/10 Q_j+2.

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 


    3^rd order   Q_j+1/2^(L) = 1/8(-Q_j-1 + 6Q_j + 3Q_j+1) 
    Q_j+1/2^(R) = 1/8(3Q_j + 6Q_j+1 - Q_j+2),
    
    5^th order   Q_j+1/2^(L)= 1/128(3Q_j-2 -20Q_j-1 +90Q_j+60Q_j+1-5Q_j+2) 
    Q_j+1/2^(R)= 1/128(-5Q_j-1 +60Q_j +90Q_j+1 -20Q_j+2 + 3Q_j+3).



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 <cit.> for details). We evaluate the numerical flux by the following upwind flux:

    𝐄_j+1/2,k =  1/2(𝐄^L + 𝐄^R) -
      1/2𝐏^-1| 𝐏𝐀|(𝐐^R-𝐐^L),

where |𝐀| is the flux Jacobian, given in Eq.(<ref>). In the local preconditioning method <cit.>, the numerical flux
  𝐄_j+1/2,k is constructed based on the preconditioned Jacobian. 𝐏^-1  is multiplied in the dissipation term in the above expression to cancel the effect of 𝐏 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.(<ref>) or (<ref>).

  * compute the upwind flux evaluation at the cell interface, Eq.(<ref>), and

  * evaluate the derivative by cell-interface-to-node differencing, Eq.(<ref>). 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.





 §.§ Time discretization

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

    𝐏^-1_j,k d 𝐐_j,k/ dτ = 𝐑𝐞𝐬_j,k,

where
 
    𝐑𝐞𝐬_j,k =  - . ∂ E_x/∂ x|_j,k - . ∂ E_y/∂ y|_j,k  + . S |_j,k,

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

    𝐐^(1)    =    𝐐^n + Δτ𝐑𝐞𝐬(𝐐^n),  
    𝐐^(2)    =    3/4𝐐^n + 1/4𝐐^(1) + 1/4Δτ𝐑𝐞𝐬(𝐐^(1)) ,
    𝐐^n+1    =    1/3𝐐^n + 2/3𝐐^(2) + 2/3Δτ𝐑𝐞𝐬(𝐐^(2)),

where Δτ is the pseudo time step given by the following equation


    Δτ = CFL min(x_i+1-x_i/|√(D_xx/T_r)|,y_i+1-y_i/|√(D_yy/T_r)|).


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.




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


  * Primary variable T

The physical domain is extended through ghost cells and the Dirichlet boundary condition is employed at the cell interface, say T_1/2. The 3^rd and 5^th   order accurate polynomials for the left and right boundaries are given by Eqs. (<ref>) and (<ref>).









    Left  boundary   
    T_0    = 1/3(8T_1/2-6T_1 + T_2),
    
    T_0   =128/35(T_1/2 -35/32 T_1 +35/64 T_2-7/32 T_3+5/128T_4).



    Right  boundary   
    T_N+1   = 1/3(8 T_N+1/2 + T_N-1-6T_N),
    
     T_N+1   =128/35(T_N+1/2 +5/128 T_N-3 -7/32 T_N-2-35/64 T_N-1-35/32T_N).

 

  * 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.(<ref>) and <ref>,  

    Left  boundary   
    g_0   = 3g_1 - 3g_2 +  g_3,
    
    g_0   = 5g_1-10g_2+10g_3-5g_4+g_5.



    Right  boundary   
    g_N+1    = 3g_N - 3g_N-1+g_N-2,
    
    g_N+1   = g_N-4 - 5g_N-3 + 10g_N-2 -10g_N-1+5g_N.




















§ 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_⊥ is taken as 1.



 §.§ Grid aligned test case

Example 4.1. The following anisotropic diffusion equation identified in Ref. <cit.>, with a spatial domain (x,y) ∈ [0,1] × [0,1],  has been considered. 


    ∂ T/∂ t = D_xx∂ ^2 T/∂ x^2 + D_yy∂ ^2 T/∂ y^2 + S,

where the diffusion tensor and source term are 
 
    𝐃 = [ [ D_xx    0;    0 D_yy;      ]] = [ [ D_||    0;    0  D_⊥;      ]] =  [ [ 10^γ    0;    0    1;      ]],     S = 50.5 sin(π x)sin(π y),

where γ is varied from 0 to 9 and the exact solution for this test case is given by, 
 
    T_exact(x,y) =1/2π^2sin(π x)sin(π y), x,y ∈ [0,1],

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 β = 0^∘, shown in Fig. <ref>. Grid refinement study is carried out for all the upwind schemes for two different relaxation times ν=1 and ν_opt. It can be seen from Fig. <ref> and Table <ref> 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 ν_opt is much lower than ν = 1. Fig. <ref> shows that the error remains the same for increasing value of D_||, which is varied from 10^0 to 10^9, for ν_opt whereas it increases for ν =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 ν_opt can make the scheme independent of anisotropy and also the angle of misalignment. The solution contours obtained by U-5E on a 32 × 32 grid and computed values at the geometric center for various schemes are shown in Figs. <ref> and <ref> respectively.












 §.§ Constant angle of misalignment 

In the following test cases the applicability of the approach to full diffusion tensor, with a constant angle of misalignment β shown in Eq.(<ref>), is analyzed. For all problems in this subsection, we set 
D_⊥=1 and D_|| = 10^γ, where γ varies again 0 to 9.

























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

    T_exact(x,y) = xy(sin(π x)sin(π y))^ω, x,y ∈ [0,1],

where ω = 10 and the angle of misalignment β is set to a constant value of 30^∘. The source term is calculated by plugging Eq.(<ref>) into Eq. (<ref>). The simulations are conducted with grid refinements from 8×8 to 128×128. Figure <ref> shows the results obtained with ν_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 ν_opt makes the schemes independent of the anisotropy. 








Figures <ref> and <ref> show the solution contours and order of accuracy for an anisotropy of 10^9 at a constant β of 30^∘ respectively. All the schemes can capture the temperature peak without any numerical diffusion.



It can be easily confirmed from the Figs. <ref>, <ref>, and <ref> 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, β, is varied from 0^∘ to 90^∘. It can be seen from Fig. <ref> that the results obtained by ν_opt are independent of anisotropy whereas the error increases with anisotropy for ν=1.



Example 4.2b. In this test case, the following temperature distribution from Ref. <cit.> is considered to test the effect of steep temperature gradients,

    T_exact(x, y) = 1/1 - e^- 2^1/2 ω(e^- 2^ω(x^2 + y^2)^ω - e^- 2^1/2 ω), x,y ∈ [-0.5, 0.5],

where  D_|| = 1, 10^3, 10^6, and 10^9 and ω is taken as 6, and the angle of misalignment is set to a constant value of 45^∘ for the anisotropy of 10^9. Finite-volume schemes proposed in Ref. <cit.> 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. <ref> that the hyperbolic scheme can also handle situations with severe anisotropy as well as sharp gradients without losing accuracy. Table <ref> show that the third order of accuracy is obtained for the primary variable, T, for all the schemes. 








 §.§ 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 ν_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] × [0, 1] for both the test cases. The preconditioned formulation  employed in this study 
is necessary for variable diffusion tensor problems <cit.>, where ν_opt is a function of space. Again, we set D_⊥=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((x-0.5)^2+(y-0.5)^2/0.01)  and the diffusion tensor given by Lou et al. <cit.> is modified by including extreme anisotropy through D_||,

    D = [ D_||(y^2 + (x + 1)^2)     - x y
    
     - x y    (y + 1)^2] = [ 10^γ(y^2 + (x + 1)^2)     - x y
    
     - x y    (y + 1)^2],









where γ =  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. <ref>. Once again, the compact scheme is more accurate than explicit interpolation. Solution contours computed on a grid size of 96×96 is shown in Fig. <ref>. 


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

    T(x,y) = (xy+(ax+by)(x^2+y^2)^3/2),

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

    D = [D_||cos^2 (β ) + D_⊥sin^2 (β )   1/2(D_|| - D_⊥) sin (2 β )
    
    1/2(D_|| -D_⊥) sin (2 β )    D_||sin^2 (β ) + D_⊥cos^2 (β )],

where  D_|| =10^9 and β = arctan(x+y).
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. <ref>. All the schemes are 3rd order accurate because of the boundary conditions as shown in Table <ref>. Solution contours computed on a grid size of 64×64 is shown in Fig. <ref>.








 §.§ 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 <cit.>, 
where ν_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., 

    0 = ∂_x (𝐃∂_x T)+∂_y (𝐃∂_y T) - S,

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

    𝐃 = [ [ D_xx    0;    0 D_yy;      ]] = [ [ D_||(1+ T^2)            0;            0        1+T^2;              ]] =  [ [ 10^γ(1+T^2)           0;           0       1+T^2;             ]]
       T_exact= sin(π x) sin(π y), 
       S = π^2(3 D_||sin^2 (π x  )sin^2 (π y  ) - 2 D_||sin^2 (π y  ) + D_|| + 3 sin^2 (π x  )sin^2 (π y  ) - 2 sin^2 (π x  ) + 1) 
       sin (π x  )sin (π y  ).

The problem reduces to isotropic diffusion, which is considered in Ref. <cit.>, for D_|| = 1. Table <ref> 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, 


    𝐃   = [ [ D_xx D_xy; D_yx D_yy;      ]] = Θ^-1[ [       D_||          0;          0 D_⊥(1+g+h);            ]] Θ,  where Θ =  [
       [  cosβ -sinβ;  sinβ  cosβ ]]





and g and h are temperature gradients ∂ T/∂ x and ∂ T/∂ y respectively. Such temperature dependence is proposed by Lafleur et al. <cit.> for magnetized plasmas in Hall thrusters. The exact solution for this problem is the same as that of the [sec:4.2a]Example 4.2a. Table <ref> 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 [sec:4.2a]Example 4.2a.
























 §.§  Magnetized electron test case:

Final numerical problem is the magnetized electron test case, which has anisotropic diffusion, considered in <cit.>. 
The hyperbolized mass and momentum equations, readers can refer the aforementioned references for further details, are given by
   
    n_ e/T_ e∂ϕ/∂ t-∇·(n_ eu⃗_ e)   =-n_ eν_ion,
    1/ν_col∂/∂ t(n_ eu⃗_ e)
          - n_ e[μ]∇ϕ+[μ]∇(n_ eT_ e)   =-n_ eu⃗_ e,

 where  n_ e, u⃗_ e = (u_x, u_y), ϕ, ν_ion, and T_ e, are the electron number density, electron velocity, space potential,  ionization collision frequency, and electron temperature respectively. The electron mobility tensor [μ], similar to [D] in Eq.(<ref>), can be expressed as follows, 
   
    [μ]=[
       [ μ_xx μ_xy; μ_xy μ_yy ]]   =Θ^-1[
       [ μ_||     ;       μ_⊥ ]]Θ,
       
       Θ=[
       [  cosθ -sinθ;  sinθ  cosθ ]], 
    
    μ_||   = e/m_ eν_col,   μ_⊥=μ_||/1+(μ_||B)^2.
        
For simplified analysis, the following values are assumed for all the test calculations,
   
    n_ e=1,T_ e=1,ν_ col=1,ν_ion=0.


 Then, the system (<ref>) can be written in the form (<ref>) with

    𝐐= [ [   ϕ; u_x; u_y ]],   𝐄_x = [ [ -u_x;  - ϕ;    0 ]],       𝐄_y = [ [ -u_y;    0;  - ϕ ]],   𝐒 = [ [    0; -û_x; -û_y ]],

where (û_x,û_y) =  [μ]^-1u⃗_ e, and 
   
    𝐏^-1 = Θ^-1P^-1Θ,   
          P^-1=(
             [        1        0        0;        0 T_r/μ_||        0;        0        0  T_r/μ_⊥ ]),
            Θ=(
       		[     1     0     0;     0  cosθ  sinθ;     0 -sinθ  cosθ ]).

Note that the relaxation time T_r has been introduced, which allows us to incorporate the dimensionally-consistent optimal length scale, i.e., 
Equation (<ref>) with Equation (<ref>), without changing steady solutions that we seek. In contrast, the formulation in the previous 
studies <cit.> corresponds to the following:
   
    𝐏^-1 
       = 
      [μ]^-1(
             [    1    0    0;    0 μ_xx    0;    0    0 μ_yy ])
           ,

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  <cit.>, 
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 𝐄_x and 𝐄_y, and 
therefore can be discretized straightforwardly by any method designed for conservative systems, including the method described in Section <ref>.
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 (<ref>) because the diffusion tensor
is not included in the first equation: it gives

    𝐏𝐀_n
    =𝐏∂(𝐄_x n_x+𝐄_y n_y)/∂ Q 
    = [
    0     -n_x     -n_y
    
    - μ_xx n_x /T_r    - μ_xy n_y /T_r    0
    
    - μ_xy n_x /T_r    - μ_xyy n_y /T_r    0
    ],

where 𝐏 is defined by Equation (<ref>).
Fortunately, it leads to the same set of eigenvalues (<ref>) for the matrix M in Equation (<ref>) when the Fourier mode is substituted
 (although the matrix itself is different).
Therefore, the relaxation time as derived in Equation (<ref>) is optimal for this system also.
The dissipation term is, however, slightly different from Equation (<ref>) because, although the eigenvalues are the same,

    Λ_𝐧 = λ[-1     0     0
    0     0     0
    0     0     1],   λ =  √(1/T_r(μ_xx n_x^2 + 2 μ_xy n_x n_y + μ_yy n_y^2)) ,

the right-eigenvectors are different:

    𝐑_n= [μ_xx n_x^2 + 2 μ_xy n_x n_y + μ_yy n_y^2/λ    0     -  μ_xx n_x^2 + 2 μ_xy n_x n_y + μ_yy n_y^2/λ
    μ_xx n_x  +  μ_xy   n_y      -n_y    μ_xx n_x  +  μ_xy   n_y
    μ_xy n_x + μ_yy n_ y     n_x    μ_xy n_x + μ_yy n_ y 
     ],

which results in the following dissipation matrix:

    P^-1 |𝐏𝐀_n|  
    =   P^-1𝐑_n|Λ_n| 𝐑^-1_n
    =
     λ[
     1       0       0 
      
    0     n_x^2/λ^2     n_x n_y/λ^2 
      
    0     n_x n_y /λ^2      n_y^2/λ^2  ],



 We consider the test case calculated by Kawashima et al. <cit.> with uniformly angled magnetic lines of force at  θ = 45^∘ from the vertical, shown in Fig. <ref>. Uniform magnetic confinement ratio, μ_||/μ_⊥ 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.(<ref>).


    ϕ̃= 
    
       1,     x= 0.0
    
        0,     x = 0.02, 2.0, 200.0
      ũ_y = 
    
       0,     y= 0.0
    
        0,     y = 0.01, 1.0, 100.0












 


The simulations are carried out by WCNS-Z scheme, as it is shown in Ref. <cit.> that linear schemes can cause spurious oscillations. The numerical results computed on a 96 × 96 grid for a mobility ratio of 10^3. Once again, WCNS-Z scheme along WENO extrapolation (see <cit.> for further details) are found suitable to reduce the spurious oscillations, as seen in Figs. <ref>. 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] × [0,100] is 
considered in Refs.<cit.>. Our scheme yields the same solution independent of the domain size.




























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

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

  * 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 <cit.> that WENO interpolation approach along with WENO extrapolation for boundary conditions can result in oscillation free results. Dimensionally consistent formulation given in <cit.> is also required for scale-invariant results.

The extension to curvilinear coordinates using WCNS<cit.> and field aligned unstructured meshes by discontinuous galerkin methods<cit.> 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+1/2 are only presented here. The WCNS numerical flux is constructed as follows:

    Q^L_j+1/2 = ∑_k=0^2ω_k Q_j+1/2^(k),

where

    Q_j+1/2^(0)   = 1/8(3Q_j-2 - 10Q_j-1 + 15Q_j),
    Q_j+1/2^(1)   = 1/8(-Q_j-1 + 6Q_j + 3Q_j+1), 
    Q_j+1/2^(2)   = 1/8(3Q_j + 6Q_j+1 - Q_j+2),

and ω_k are the nonlinear weights. These nonlinear weights are computed as follows:

    ω_k = α_k/∑_k=0^2α_k,  where α_k = C_k/(β_k + ϵ)^m,   k = 0, 1, 2,

and

    C_0 = 1/16,  
    C_1 = 10/16,  
    C_2 = 5/16,


    β_0    = 1/4( Q_i-2 - 4 Q_i-1 +3 Q_i)^2 + 13/12( Q_i-2 - 2 Q_i-1 + Q_i)^2
    β_1    = 1/4( Q_i-1 - Q_i+1)^2  + 13/12( Q_i-1 - 2 Q_i + Q_i+1)^2 
    β_2    = 1/4( 3Q_i - 4 Q_i+1 +Q_i+2)^2+ 13/12( Q_i - 2 Q_i+1 + Q_i+2)^2.

Here ϵ=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 <cit.>. The improved nonlinear weights are as follows:


    ω^z_k=α^z_k∑_k=0^2α^z_k,  α^z_k=γ_k(1+(τϵ+β_k)^p),

 
    τ=|β_0-β_2|,

where the smoothness indicators β_k's are the same as those given in Equations (<ref>), ϵ=10^-40, and τ 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. <cit.> 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.




§ REFERENCES

elsarticle-num

