High-order upwind and non-oscillatory approach for steady state diffusion, advection-diffusion and application to magnetized electrons

TL;DR

This paper develops high-order upwind and WENO schemes for steady state magnetized electron fluid equations, demonstrating their accuracy and ability to capture sharp gradients without oscillations, especially under strong magnetic confinement.

Contribution

The paper introduces a novel application of higher order upwind and WENO schemes to magnetized electron fluid equations and advection-diffusion problems, showing their effectiveness in capturing sharp gradients.

Findings

Linear schemes suffice for diffusion equations.

WENO schemes do not cause adverse effects on solutions.

WENO captures sharp gradients effectively in magnetized electron simulations.

Abstract

Steady state simulations} of magnetized electron fluid equations with strong anisotropic diffusion based on the first-order hyperbolic approach is carried out using cell-centered higher order upwind schemes, linear and weighted essentially non-oscillatory (WENO). Along with the magnetized electrons, the diffusion equation is also simulated to demonstrate the implementation and design order of the accuracy of the approach due to their similar upwind structure. We show the adequacy of linear upwind schemes for diffusion equation and the use of shock-capturing scheme like WENO does not have any adverse effect on the solution, unlike the total-variation diminishing (TVD) methods. We further extended the approach to advection-diffusion equation, and appropriate boundary conditions have obtained a consistent design accuracy of the third and fifth order. We implemented the WENO approach to advection-diffusion equation by using the split hyperbolic method to demonstrate the advantage of non-oscillatory schemes to capture sharp gradients in boundary layer type problems without spurious oscillations. Finally, numerical results for magnetized electrons simulations indicate that with increasing strength of magnetic confinement it is possible to capture sharp gradients without oscillations by WENO scheme.

Tables

Table 1. Table 1: L 2 subscript 𝐿 2 L_{2} errors and order of convergence of u 𝑢 u , primary variable, by 3rd, 5th order explicit and compact and WENO schemes for 1D diffusion equation.

Number	Upwind-3E		Upwind-5E		Upwind-5C		WENO-5Z
of points	error	order	error	order	error	order	error	order
24	1.54E-02		1.83E-03		1.26E-03		2.64E-03
48	2.11E-03	2.87	9.70E-05	4.24	1.10E-04	3.52	7.70E-05	5.10
96	2.63E-04	3.01	4.02E-06	4.59	4.84E-06	4.50	3.99E-06	4.27
192	3.23E-05	3.02	1.33E-07	4.91	1.61E-07	4.91	1.33E-07	4.90
384	4.00E-06	3.02	4.22E-09	4.98	5.07E-09	4.99	4.22E-09	4.98

Table 2. Table 2: L 2 subscript 𝐿 2 L_{2} errors and order of convergence of gradient variable by 3rd, 5th order explicit and compact and WENO schemes one-dimensional diffusion problem, Example 4.1.1 .

Number	Upwind-3E		Upwind-5E		Upwind-5C		WENO-5Z
of points	error	order	error	order	error	order	error	order
24	5.82E-03		4.66E-04		1.64E-04		1.44E-03
48	6.91E-04	3.07	2.79E-05	4.06	2.92E-05	2.48	2.73E-05	5.72
96	6.48E-05	3.42	8.23E-07	5.08	8.86E-07	5.04	8.23E-07	5.05
192	5.83E-06	3.47	1.94E-08	5.40	2.10E-08	5.40	1.94E-08	5.40
384	5.31E-07	3.46	4.40E-10	5.47	4.71E-10	5.48	4.40E-10	5.46

Table 3. Table 3: L 2 subscript 𝐿 2 L_{2} errors and order of convergence of primary variable, u 𝑢 u , by 3rd order explicit, 5th order explicit and compact and WENO schemes one-dimensional diffusion problem, Example 4.1.1 .

Number	Upwind-3E		Upwind-5E		Upwind-5C		WENO-5Z
of points	error	order	error	order	error	order	error	order
16	7.01E-04		5.46E-06		6.48E-06		6.59E-04
32	8.24E-05	3.09	2.63E-07	4.37	2.97E-07	4.45	6.20E-05	3.41
64	9.47E-06	3.12	1.00E-08	4.72	1.09E-08	4.77	6.60E-06	3.23
128	1.10E-06	3.10	3.44E-10	4.86	3.64E-10	4.90	7.23E-07	3.19
256	1.31E-07	3.07	1.13E-11	4.93	1.14E-11	5.00	8.19E-08	3.14

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} errors and order of convergence of potential by 3rd order explicit, 5th order explicit and compact and WENO schemes for Neumann boundary test case.

Number	Upwind-3E		Upwind-5E		Upwind-5C		WENO-5Z
of points	error	order	error	order	error	order	error	order
${16}^2$	2.01E-03		3.23E-05		4.11E-05		3.22E-05
${32}^2$	2.69E-04	2.90	1.42E-06	4.50	1.68E-06	4.61	1.42E-06	4.50
${64}^2$	3.29E-05	3.03	5.61E-08	4.67	6.12E-08	4.78	5.61E-08	4.67
${128}^2$	3.93E-06	3.06	1.99E-09	4.82	2.06E-09	4.89	1.99E-09	4.82
${256}^2$	4.73E-07	3.06	6.63E-11	4.91	6.68E-11	4.95	6.62E-11	4.91

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for two-dimensional Poisson problem, Example 4.1.4 , by 3rd order explicit, 5th order explicit and compact and WENO schemes.

Number	Upwind-3E		Upwind-5E		Upwind-5C		WENO-5Z
of points	error	order	error	order	error	order	error	order
${16}^2$	2.466E-01		3.06E-02		2.65E-02		6.15E-02
${32}^2$	3.221E-02	2.94	1.16E-03	4.72	7.32E-04	5.18	1.30E-03	5.56
${64}^2$	4.006E-03	3.01	3.42E-05	5.09	1.35E-05	5.76	3.52E-05	5.21
${128}^2$	4.985E-04	3.01	1.01E-06	5.09	2.30E-07	5.87	1.01E-06	5.12
${256}^2$	6.218E-05	3.00	3.03E-08	5.05	4.01E-09	5.84	3.01E-08	5.07

Table 6. Table 6: L 2 subscript 𝐿 2 L_{2} errors and order of convergence by 3rd order explicit, 5th order explicit and compact and WENO schemes for anisotropic diffusion problem given in Example 4.1.5 .

Number	Upwind-3E		Upwind-5E		Upwind-5C		WENO-5Z
of points	error	order	error	order	error	order	error	order
${16}^2$	4.99E-04		7.47E-06		5.68E-06		7.61E-06
${32}^2$	6.23E-05	3.00	1.74E-07	5.42	8.61E-08	6.04	1.92E-07	5.31
${64}^2$	7.83E-06	2.99	4.43E-09	5.30	1.35E-09	5.99	4.42E-09	5.44
${128}^2$	9.81E-07	3.00	1.23E-10	5.17	2.29E-11	5.88	1.23E-10	5.17
${256}^2$	1.23E-07	3.00	3.68E-12	5.06	4.22E-13	5.77	3.68E-12	5.06

Table 7. Table 7: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for Example 4.2.2 by using WENO scheme along with WENO and Lagrange extrapolation techniques for boundary conditions.

Number	Advection				Advection-Diffusion
of points	WENO-5Z-L		WENO-5Z-W		WENO-5Z-L		WENO-5Z-W
	error	order	error	order	error	order	error	order
$8^2$	8.36E-03		7.67E-03		4.55E-02		3.99E-02
${16}^2$	1.36E-03	2.62	1.32E-03	2.53	7.83E-03	2.54	1.35E-02	1.56
${32}^2$	1.41E-04	3.27	1.41E-04	3.23	1.30E-03	2.59	4.92E-03	1.45
${64}^2$	1.56E-05	3.18	1.56E-05	3.18	1.48E-04	3.14	1.75E-03	1.49
${128}^2$	1.73E-06	3.17	1.73E-06	3.17	1.50E-05	3.30	7.51E-04	1.22

Table 8. Table 8: L 2 subscript 𝐿 2 L_{2} errors and order of convergence for Example 4.2.3 by using WENO scheme along with WENO and Lagrange extrapolation techniques for boundary conditions.

Number	WENO-5Z-W		WENO-5Z-L
of points	error	order	error	order
${20}^2$	2.21E-03		1.91E-03
${40}^2$	1.07E-03	1.04	2.00E-04	3.26
${80}^2$	4.60E-04	1.22	2.34E-05	3.09
${160}^2$	1.99E-04	1.21	2.59E-06	3.18
${320}^2$	9.06E-05	1.14	2.68E-07	3.27

Table 9. Table 9: Error in transverse electron flux for μ | | / μ ⟂ {\mu}_{||}/{\mu}_{\perp} =1000 and θ 𝜃 \theta =45 ∘ .

Number	WENO-5Z-W		MFAM
of points	${\Gamma }_e$	Error $\%$	${\Gamma }_e$	Error $\%$
96$\times$48	-0.01861	1.12	-0.01878	0.21
192$\times$96	-0.01878	0.21	-0.01879	0.16
384$\times$192	-0.01881	0.05	-0.01882	0.0
768$\times$384	-0.01882	0.0	-0.01882	0.0

Table 10. Table 10: Relative l 2 subscript 𝑙 2 l_{2} error for MFAM and WENO-5Z-W for μ | | / μ ⟂ {\mu}_{||}/{\mu}_{\perp} =1000.

Number of points	$\tilde{\varphi }$	${\tilde{u}}_x$	${\tilde{u}}_y$
96$\times$48	3.28E-03	1.01E-02	1.03E-02
192$\times$96	2.62E-03	9.92E-03	1.02E-02
384$\times$192	1.73E-03	8.71E-03	8.92E-03
768$\times$384	9.82E-04	8.12E-03	8.35E-03

Equations

$$\nabla \cdot \left(n_{\mathrm{e}}{\vec{u}}_{\mathrm{e}}\right)=n_{\mathrm{e}}{\nu }_{\mathrm{ion}},$$

$$n_{\mathrm{e}}\left[\mu \right]\nabla \varphi -\left[\mu \right]\nabla \left(n_{\mathrm{e}}{T}_{\mathrm{e}}\right)=n_{\mathrm{e}}{\vec{u}}_{\mathrm{e}},$$

$$\displaystyle\left[\mu\right]=\left[\begin{array}[]{cccc}\mu_{x}&\mu_{\rm c}\\ \mu_{\rm c}&\mu_{y}\end{array}\right]$$

$${\mu }_{||}$$

$$\nabla \cdot \left(n_{\mathrm{e}}\left[\mu \right]\nabla \varphi -\left[\mu \right]\nabla \left(n_{\mathrm{e}}{T}_{\mathrm{e}}\right)\right)=n_{\mathrm{e}}{\nu }_{\mathrm{ion}}.$$

$$\frac{n_{\mathrm{e}}}{T_{\mathrm{e}}}\frac{\partial \varphi }{\partial t}-\nabla \cdot \left(n_{\mathrm{e}}{\vec{u}}_{\mathrm{e}}\right)$$

$$\frac{1}{{\nu }_{\mathrm{col}}}\frac{\partial }{\partial t}\left(n_{\mathrm{e}}{\vec{u}}_{\mathrm{e}}\right)-n_{\mathrm{e}}\left[\mu \right]\nabla \varphi +\left[\mu \right]\nabla \left(n_{\mathrm{e}}{T}_{\mathrm{e}}\right)$$

$${\tilde{n}}_{\mathrm{e}}=\frac{n_{\mathrm{e}}}{n_{\mathrm{e}}^{\ast }},{\tilde{T}}_{\mathrm{e}}=\frac{T_{\mathrm{e}}}{T_{\mathrm{e}}^{\ast }},\tilde{\varphi }=\frac{\varphi }{{\varphi }_{\mathrm{a}}},{\left(\tilde{x},\tilde{y}\right)}^T=\frac{1}{{\lambda }_{\mathrm{m}}^{\ast }}{\left(x,y\right)}^T,$$

$$\tilde{t}=\frac{1}{{\lambda }_{\mathrm{m}}^{\ast }}\sqrt{\frac{2e{T}_{\mathrm{e}}^{\ast }}{m_{\mathrm{e}}}}t,{\widetilde{\stackrel{⃗}{u}}}_{\mathrm{e}}=\frac{{\vec{u}}_{\mathrm{e}}}{c_{\mathrm{s}}^{\ast }}=\frac{{\vec{u}}_{\mathrm{e}}}{\sqrt{\frac{\gamma e{T}_{\mathrm{e}}^{\ast }}{m_{\mathrm{e}}}}},{\tilde{\nu }}_{\mathrm{col}}={\tau }_{\mathrm{m}}^{\ast }{\nu }_{\mathrm{col}},{\tilde{\nu }}_{\mathrm{ion}}={\tau }_{\mathrm{m}}^{\ast }{\nu }_{\mathrm{ion}}.$$

$${\tilde{n}}_{\mathrm{e}}=1,{\tilde{T}}_{\mathrm{e}}=1,{\tilde{\nu }}_{\mathrm{col}}=1,{\tilde{\nu }}_{\mathrm{ion}}=0.$$

$$\frac{\partial \tilde{\varphi }}{\partial \tilde{t}}-\sqrt{\frac{\gamma }{2}}\tilde{\nabla }\cdot {\overrightarrow{\stackrel{~}{u}}}_{\mathrm{e}}$$

$$\frac{\partial {\overrightarrow{\stackrel{~}{u}}}_{\mathrm{e}}}{\partial \tilde{t}}-\frac{1}{\sqrt{2\gamma }}\left[\tilde{\mu }\right]\tilde{\nabla }\tilde{\varphi }$$

$${\hat{\mathbf{P}}}^{-1}\frac{\partial Q}{\partial \tilde{t}}+\frac{\partial {\mathbf{E}}_{\mathbf{x}}}{\partial \tilde{x}}+\frac{\partial {\mathbf{E}}_{\mathbf{y}}}{\partial \tilde{y}}=S,$$

$$P=\left(\begin{array}[]{ccc}\sqrt{\frac{2}{\gamma}}&0&0\\ 0&\frac{\sqrt{2\gamma}}{\tilde{\mu}_{||}}&0\\ 0&0&\frac{\sqrt{2\gamma}}{\tilde{\mu}_{\perp}}\end{array}\right),\hskip 20.0pt\Theta=\left(\begin{array}[]{ccc}1&0&0\\ 0&\cos\theta&\sin\theta\\ 0&-\sin\theta&\cos\theta\end{array}\right).$$

$$\frac{\partial \tilde{\varphi }}{\partial \tilde{t}}-\nabla \cdot {\overrightarrow{\stackrel{~}{u}}}_{\mathrm{e}}$$

$$\frac{\partial {\overrightarrow{\stackrel{~}{u}}}_{\mathrm{e}}}{\partial \tilde{t}}-\nabla \tilde{\varphi }$$

$$\mathbf{[g]}=\frac{\sqrt{2\gamma}}{{\tilde{\mu}_{||}}}\left[\begin{array}[]{cc}\cos^{2}\theta+\frac{{\tilde{\mu}_{||}}}{{\tilde{\mu}_{\perp}}}\sin^{2}\theta&(\frac{{\tilde{\mu}_{||}}}{{\tilde{\mu}_{\perp}}}-1)\cos\theta\sin\theta\\ (\frac{{\tilde{\mu}_{||}}}{{\tilde{\mu}_{\perp}}}-1)\cos\theta\sin\theta&\sin^{2}\theta+\frac{{\tilde{\mu}_{||}}}{{\tilde{\mu}_{\perp}}}\cos^{2}\theta\end{array}\right].$$

$${\mathbf{\hat{PJ}}}_{x}=\left(\begin{array}[]{ccc}0&-1&0\\ -1&0&0\\ 0&0&0\end{array}\right),\hskip 20.0pt{\mathbf{\hat{PJ}}}_{y}=\left(\begin{array}[]{ccc}0&0&-1\\ 0&0&0\\ -1&0&0\end{array}\right),\quad\lambda=\pm 1,0.$$

$$\frac{\partial u}{\partial t}-\nu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$$

$$\displaystyle{\color[rgb]{0,0,0}\frac{\partial u}{\partial\tau}}-\nu\frac{\partial p}{\partial x}-\nu\frac{\partial q}{\partial y}$$

$$\displaystyle{\color[rgb]{0,0,0}\frac{\partial p}{\partial\tau}}-\frac{\partial u}{\partial x}$$

$$\displaystyle{\color[rgb]{0,0,0}\frac{\partial q}{\partial\tau}}-\frac{\partial u}{\partial y}$$

$$\{\begin{array}{l}\begin{array}{rl}\text{cancelto}\frac{\partial u}{\partial \tau }-\nu \frac{\partial p}{\partial x}-\nu \frac{\partial q}{\partial y}& =-S,\\ \text{cancelto}\frac{\partial p}{\partial \tau }-\frac{\partial u}{\partial x}& =-p,\\ \text{cancelto}\frac{\partial q}{\partial \tau }-\frac{\partial u}{\partial y}& =-q.\end{array}\end{array}\to$$

$${\color[rgb]{0,0,0}\frac{\mathbf{\partial Q}}{\partial{\tau}}}+\frac{\mathbf{\partial E_{x}}}{\partial{x}}+\frac{\mathbf{\partial E_{y}}}{\partial{y}}=\mathbf{S},$$

$$\mathbf{Q}=\left[\begin{array}[]{ccc}u\\ p\\ q\end{array}\right],\quad\mathbf{E_{x}}=\left[\begin{array}[]{ccc}-\nu p\\ -u\\ 0\end{array}\right],\quad\mathbf{E_{y}}=\left[\begin{array}[]{ccc}-\nu q\\ 0\\ -u\end{array}\right],\quad\mathbf{S}=\left[\begin{array}[]{ccc}-S\\ -p\\ -q\end{array}\right],$$

$${\bf J}_{x}=\nu\left(\begin{array}[]{ccc}0&-1&0\\ -1&0&0\\ 0&0&0\\ \end{array}\right),\quad{\bf J}_{y}=\nu\left(\begin{array}[]{ccc}0&0&-1\\ 0&0&0\\ -1&0&0\\ \end{array}\right),\quad\lambda=\pm\sqrt{\nu},0.$$

$$|{\bf J}_{dx}|=\bf R_{x}|\Lambda|L_{x}=\sqrt{\nu}\left(\begin{array}[]{ccc}1&0&0\\ 0&1&0\\ 0&0&0\\ \end{array}\right).$$

$${\hat{E}}_{i+\frac{1}{2},j}$$

$$\displaystyle=\frac{1}{2}\left(\begin{array}[]{ccc}-(p_{R}+p_{L})\\ -(u_{R}+u_{L})\\ 0\\ \end{array}\right)-\frac{\sqrt{\nu}}{2}\left(\begin{array}[]{ccc}1&0&0\\ 0&1&0\\ 0&0&0\\ \end{array}\right)\left(\begin{array}[]{ccc}u_{R}-u_{L}\\ p_{R}-p_{L}\\ q_{R}-q_{L}\\ \end{array}\right).$$

$${\hat{E}}_{i+\frac{1}{2},j}=\frac{1}{2}({\mathbf{E}}_{\mathbf{L}}+{\mathbf{E}}_{\mathbf{R}})-\frac{1}{2}{\hat{\mathbf{P}}}^{-1}|\hat{\mathbf{P}}{\hat{\mathbf{J}}}_{\mathbf{x}}|({\mathbf{Q}}_{\mathbf{R}}-{\mathbf{Q}}_{\mathbf{L}}).$$

$$\mathbf{E}(Q)=\frac{\partial \mathbf{E}}{\partial \mathbf{Q}}Q=J(\mathbf{Q})Q.$$

$$\mathbf{E}$$

$${\mathbf{E}}^{+}$$

$${\mathbf{E}}^{-}$$

Full text

1 Introduction

An upwind hyperbolic-equation-system approach for two-dimensional electron fluid equations has been first proposed by Kawashima et al. [1] for robust calculation on a vertical-horizontal uniform mesh. While this work is a novel attempt for such problems, there are few drawbacks. The underlying scheme for the original paper was first order upwind scheme which was highly diffusive and might require extremely fine mesh which can be computationally expensive. Also, the original work was focused on improving diagonal dominance by avoiding cross-diffusion terms and speed of the computation but did not specifically address the issue of spurious oscillations due to misalignment of the grid and magnetic field that may occur with highly anisotropic diffusion problems on non-aligned meshes. However, mesh convergence carried out by evaluating the transverse electron flux for both the hyperbolic approach and field-aligned approach indicated that the numerical diffusion, which can be inferred as spurious oscillations, is reduced on very fine meshes. It has been proven that the hyperbolic approach attains good efficiency and accuracy in the presence of strong magnetic confinement. Kawashima et al. [2] had later extended their approach to non-isothermal system by including energy equation. They considered $3^{rd}$ order TVD scheme to reduce the numerical diffusion, but it will be shown in this paper that the TVD scheme will reduce to third order linear scheme and can still have considerable numerical errors. TVD scheme may not be able to capture the sharp gradients without oscillations. The issues that are addressed in this paper are:

• •

spurious oscillations due to misalignment of the grid and magnetic field

• •

implementation of essentially non-oscillatory schemes to reduce numerical diffusion

• •

implementation of boundary conditions

The primary influence for the magnetized electron fluid simulations is the upwind formulation for diffusion equation introduced by Nishikawa [3] based on residual-distribution (RD) method. 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 the diffusion equation is computed as a hyperbolic system. This radical approach has been shown to offer several advantages over conventional methods, like accelerated convergence for steady state solution, higher order of accuracy for both primary and gradient variables. The original approach of Nishikawa has been further extended to edge-based finite volume schemes [4] and Lee et al. have introduced cell-centered finite volume approach [5]. This approach was also extended to advection-diffusion equation by Nishikawa [6, 7], to time-dependent problems by Mazaheri and Nishikawa [8] and Navier-Stokes equation [9]. In addition to the magnetized electron simulations, in this paper, we extend the upwind formulation introduced for the diffusion and advection-diffusion equations to higher order and high-resolution cell centered schemes on uniform meshes. The similar upwind structure of diffusion and magnetized electron fluid equations and the ability of WENO scheme to capture the sharp gradients in advection-diffusion without oscillations provide essential insights for the plasma simulations.

The accuracy of the finite volume method (FVM) depends on the accuracy of the numerical fluxes computed at the cell interfaces. These numerical fluxes are evaluated based on the flux Jacobian matrices and variables defined at the cell-centers. For a hyperbolic equation the upwind fluxes at the cell interfaces, say $\hat{E}_{i+\frac{1}{2}}$, are constructed either by a Godunov type approximate Riemann solver(Roe [10], HLLC [11], etc.) or Boltzmann type solver, also known as flux vector splitting method (AUSM [12], Steger-Warming [13] etc.). According to the Godunov’s theorem [14]: “For simulations of flows involving discontinuities and sharp gradients, Gibbs phenomenon or spurious oscillations appear in the solutions near the discontinuities if the computations are carried out by linear numerical schemes that are greater than first-order accurate”.

In the last three decades, several high-order accurate schemes were developed to satisfy this criterion, and among these methods, TVD and WENO are most popular. The total-variation diminishing monotone upstream scheme for conservation laws (TVD-MUSCL) [15] is one of the popular high-order space accuracy schemes to compute the fluxes because of its simplicity and robustness. A well-known drawback of the slope limiters that are associated with these methods is that they tend to ‘clip’ smooth extrema of the flow and the accuracy necessarily degenerates to first order. On the other end, WENO schemes first introduced by Jiang and Shu [16] to capture discontinuities without spurious oscillations and are also able to achieve an arbitrarily high formal order of accuracy in smooth flows. Various versions, like WENO-5M [17], WENO-3YC [18] and WENO-5Z [19], of the WENO scheme, are proposed over the years to improve the accuracy and reduce their dissipative nature. Tan and Shu [20] had also proposed high-order boundary conditions based on the WENO idea to prevent numerical oscillations contaminating the solution due to shockwaves near the boundary. Another well-known high-order method is the family of central compact schemes developed by Lele [21]. Unfortunately, the compact central schemes cannot be used in the current ‘upwind’ formulation for diffusion and also cannot capture the discontinuities without oscillations. Pirozzoli [22] had proposed a conservative compact upwind method in combination with WENO scheme in-order to capture discontinuities. Recently, Ghosh and Baeder [23] has developed a class of upwind biased compact-reconstruction finite difference WENO schemes called CRWENO. This concept was extended in [24] where a positivity-preserving fifth-order finite volume compact-WENO (FVCW) scheme was developed. On the other side, an alternative approach is developed by Deng et al. [25] known as weighted compact non-linear schemes (WCNS) with similar discontinuity capturing abilities of WENO. These schemes are more flexible that we can interpolate not only fluxes [26], but also conservative variables [27], primitive variables or variables that are projected to the characteristic fields [28]. The other advantages of WCNS are that they can be used with flux difference splitting methods like Roe and HLLC and still maintain high-order of accuracy and also have good freestream and vortex preservation capabilities on curvilinear grids [29].

In this paper, we present the adequacy of linear upwind schemes for diffusion equation, and the use of shock-capturing scheme like WENO does not have any adverse effect on the solution, unlike TVD methods. While the third order explicit and fifth order explicit, compact and WENO reconstructions are presented in this paper, the schemes can be extended to arbitrarily high-order [30, 31, 32]. We can obtain a significant increase in efficiency because the high-order method can attain better results than the low order method on a coarse mesh. We further extended the approach to advection-diffusion equation, and appropriate boundary conditions have obtained a consistent design accuracy of third and fifth order. We implemented the WENO approach to advection-diffusion equation by using the split hyperbolic method to demonstrate the advantage of non-oscillatory schemes to capture sharp gradients in boundary layer type problems without spurious oscillations. Through detailed analysis, it is shown that the WENO scheme can also provide a robust approach in capturing sharp gradients in strongly anisotropic diffusion problems on a non-magnetic field aligned mesh. The objective is to incorporate the advantages of hyperbolic approach, higher order schemes and WENO methodology to develop a robust and efficient method that can also be extended to other equations like incompressible and compressible Navier - Stokes, non-neutral plasma simulations, multi-phase flows, etc. in future.

The rest of the paper is organized as follows. Hyperbolic approach and upwind formulation of electron fluid equations, diffusion equation, and advection-diffusion equation are described in Section 2. The cell-centered explicit and compact upwind schemes along with WENO and the implementation of boundary conditions are presented in Section 3. Several test cases in one and two-dimensional problems for diffusion equation and advection-diffusion equation are presented in Section 4. These numerical experiments validate our numerical schemes and corroborate the high-order accuracy, oscillation-free performance and implementation of boundary conditions. Finally, the simulations magnetized electron fluids are discussed in Section 4 and Section 5 summarizes our conclusions and provide suggestions for future work.

2 Governing equations

2.1 Magnetized electron fluid equations

The system of equations, electron mass and momentum equations, for the magnetized electron fluid in quasi-neutral flow described by Kawashima et al. [1] are given by

[TABLE]

[TABLE]

where $n_{\rm e}$, $\vec{u}_{\rm e}$, $\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]$ can be expressed as follows,

[TABLE]

where $e$, $B$, and $\nu_{\rm{col}}$ are the elemental charge, magnetic flux density, and electron-neutral total collision frequency, respectively. Also, $\Theta$ is the rotation matrix, and $\theta$ is the angle between the magnetic lines of force and the grid, shown in Fig. 1. The electron flux in parallel ($||$) and perpendicular ($\perp$) directions of the magnetic lines of force are described by using the electron mobility $\mu_{||}$ and $\mu_{\perp}$ respectively.

For steady state condition, the following plasma diffusion equation is obtained by substituting momentum equation 2 in to continuity equation 1

[TABLE]

For a given electron number density and temperature distribution, this equation can then be solved for plasma potential, $\phi$, using a general Poisson’s equation solver [33, 34]. Equation (4) is usually solved by a central scheme [35, 36] and it is a well known fact that the central scheme can lead to numerical oscillations if sharp gradients are present in the flow. This approach also suffers from poor iterative convergence due to the large disparity between $\mu_{\rm||}$ and $\mu_{\rm\perp}$. To avoid these difficulties Kawashima et al. [1] have proposed hyperbolic approach by introducing pseudo-time terms in the continuity and momentum equations. This approach is similar to the artificial compressibility method utilized in solving incompressible Navier-Stokes [37], where a pseudo-time derivative of the pressure is added to the continuity equation. Similarly, a pseudo-time derivative of plasma potential is added in the continuity equation and the corresponding equations are as follows,

[TABLE]

For the simulations, the equations are expressed in non-dimensional form. The non-dimensionalization procedure of the equations 5 and the definitions of dimensionless quantities are same as that of Ref. [1] which are as follows:

[TABLE]

[TABLE]

where $\tau_{\rm m}$, $v_{\rm{e,th}}$, $c_{\rm s}$ and $\gamma$ are the mean free time, electron thermal velocity, electron acoustic velocity, and specific heat ratio respectively. Also, the electron number density $n_{\rm e}^{*}$, electron temperature ${T_{\rm e}^{*}}$, plasma potential ${\phi_{a}}$, and mean free path $\lambda_{\rm m}^{*}$ are the reference values. For simplified analysis the following values are assumed for all the test calculations,

[TABLE]

Finally, the simplified nondimensional system of electron fluids for the Eqs. (5) can be expressed as follows:

[TABLE]

2.1.1 Preconditioned system

Due to the significant discrepancy between parallel, $\mu_{||}=10^{3}$, and orthogonal, $\mu_{\perp}=1$, mobilities, the condition number of the system increases significantly and degrades the convergence performance. Convergence can be made independent of mobility by altering the eigenvalues of the system such that all of them are of the same order. A popular approach is to multiply the system of equations by a preconditioning matrix to normalize the eigenvalues [38]. The system can be written in the preconditioned form as,

[TABLE]

where $\mathbf{\hat{P}}$ = $\Theta^{-1}P\Theta$ is a preconditioning matrix which is derived based on the idea of electron mobility tensor rotation

[TABLE]

Finally, the preconditioned system of equations can be rewritten as follows:

[TABLE]

where, $\mathbf{[g]}$ is a tensor which is written as follows:

[TABLE]

The flux Jacobian matrices and the corresponding eigenvalues for the preconditioned system of equations can be expressed as follows,

[TABLE]

Therefore, the Jacobian matrices are not affected by the magnetic field as they are included in the source terms and also the eigenvalues are significantly simplified which improves the convergence speed. The construction of upwind fluxes will be explained in the next subsection along with the diffusion equation.

2.2 Diffusion equation and construction of hyperbolic scheme

In this subsection, the hyperbolic approach for diffusion equation, first described by Nishikawa [3], is briefly explained here. Consider the following diffusion equation in two dimensions,

[TABLE]

where $\nu$ is the diffusion coefficient and S is the source term. In the original idea of Nishikawa [3], new variables representing the gradients of the primary variable are introduced, and the original diffusion equation has been converted into a system of three coupled first-order equations. By defining new variables $p$ = $\frac{\partial u}{\partial x}$ and $q$ = $\frac{\partial u}{\partial y}$, and introducing pseudo-time terms we can obtain the following first order hyperbolic system:

[TABLE]

It is important to state the fact that the first-order system reduces to the diffusion equation at the steady state, i.e. the pseudo-time terms will be zero at steady state,

[TABLE]

The Eqs. (16) can be represented in vector form as,

[TABLE]

where the conservative variables, fluxes in x and y-direction and source terms are,

[TABLE]

respectively. The Jacobian matrices, $\mathbf{J_{x}}$ and $\mathbf{J_{y}}$, in the x- and y-directions and the corresponding eigenvalues, $\lambda$, are given by

[TABLE]

We can notice that the above Jacobian matrices are the same as that of the hyperbolic form of magnetized electron equations derived in the earlier section, Eq. (14). The upwind fluxes are the same for both set of equations, and by changing the source terms and by using appropriate boundary conditions the magnetized electron fluid equations and diffusion can be simulated by the same code. In steady state, the diffusion equation will reduce to Laplace’s or Poisson’s equation. By the replacing the primary variable $u$ with $\phi$, space potential, the hyperbolic formulation can be used to model many physics problems that are time independent. By using the right, $\mathbf{R}$, and left, $\mathbf{L}$, eigenvectors of Jacobian matrices, and the diagonal eigenvalue matrix, $\mathbf{|\Lambda|}$, the absolute Jacobian or the Roe matrix in x-direction, $\mathbf{|J_{dx}|}$, can be written as,

[TABLE]

Finally, the upwind flux in the x-direction can be expressed as Eq. (22),

[TABLE]

where the left and right fluxes, $\mathbf{E_{L}}$ and $\mathbf{E_{R}}$, are defined at the cell interfaces which are evaluated by the interpolation polynomials that are discussed later. Similarly, the upwind flux for the magnetized electrons will be constructed by following the standard procedure of local preconditioning method [39], i.e., the preconditioned Jacobian is multiplied by $\mathbf{P}^{-1}$ to cancel the effect of $\mathbf{P}$,

[TABLE]

2.2.1 Homogeneity and flux vector splitting

Interestingly, the hyperbolic form of the diffusion equation 16 satisfy the homogeneity property,

[TABLE]

The proof of this property is easy to notice. By multiplying the Jacobian matrices shown in Eq. (20) by the vector Q we can reproduce the flux vector E(Q). This remarkable property of the hyperbolic form can be used to solve the diffusion equation by employing Flux vector splitting schemes. For example, in the original work of Kawashima et al. [1] the upwind fluxes are computed by using Steger-Warming [13] flux vector splitting scheme, shown in Eq. (25)

[TABLE]

with,

[TABLE]

where ${\bf J}^{+}$ and ${\bf J}^{-}$ in x-direction are expressed as follows,

[TABLE]

and $\Lambda^{+}$ is the diagonal matrix of positive eigenvalues, and $\Lambda^{-}$ is the diagonal matrix of negative eigenvalues. The Jacobian matrices satisfy, ${\bf J}^{+}+{\bf J}^{-}={\bf J}$ and the reconstruction for $\mathbf{E^{+}}$ uses a biased stencil with one more point to the left, and that for $\mathbf{E^{-}}$ uses a biased stencil with one more point to the right, to obey correct upwinding. In the current upwind formulation for diffusion equation, there is no difference between flux vector splitting and Roe solver. For all the simulations in this paper, Roe fluxes are used unless otherwise stated.

2.3 Advection-Diffusion equation and construction of hyperbolic scheme

In this section, we consider the advection-diffusion equation,

[TABLE]

where where $u$ is the solution variable, (a, b) is a constant advection vector, and $\nu$ is a constant positive diffusion coefficient. Similar to the diffusion equation new variables $p$ = $\frac{\partial u}{\partial x}$ and $q$ = $\frac{\partial u}{\partial y}$ are defined and pseudo-time derivatives are introduced to obtain the following first order hyperbolic system:

[TABLE]

The Eqs. (29) can be represented in vector form as,

[TABLE]

where the conservative variables, fluxes in x and y-direction and source terms are,

[TABLE]

respectively. The split hyperbolic formulation introduced by Nishikawa in Ref. [7] is considered in this paper. The advection, $E^{a}_{x}$, and diffusion, $E^{d}_{x}$, fluxes in x-direction are separated as

[TABLE]

The absolute Jacobian matrix for the diffusive fluxes would be the same as in Eq. (21) and the corresponding absolute Jacobian matrix for the advection flux is straightforward

[TABLE]

and finally, the upwind flux can be constructed as,

[TABLE]

The upwind fluxes can also be constructed by the unified advection-diffusion approach proposed by Nishikawa in Ref. [6] and the numerical schemes presented here can be used in a straight forward manner and is currently beyond the scope of the current paper.

3 Numerical method

In this section, the procedures of WENO and linear upwind schemes are briefly explained. From the Eq. (34) we can see that the computation of $E_{i+\frac{1}{2}}$ requires knowledge of the values $q_{R}$ and $q_{L}$ at the cell interface. The values of $q_{R}$ and $q_{L}$ can be obtained by upwind interpolation to the same order of accuracy. The ${Q}_{i+\frac{1}{2}}^{(L)}$ should be biased to left and similarly ${Q}_{i+\frac{1}{2}}^{(R)}$ has to be biased to the right for correct upwinding, shown in Fig. 2.

3.1 Weighted Essentially Non-Oscillatory Schemes

In the WENO scheme, the fifth order upwind-biased interpolation is nonlinearly weighted from three different third order interpolations on sub-stencils, $S_{0}$, $S_{1}$ and $S_{2}$, shown in Fig. 2. For simplicity, the interpolation polynomials to the left side of the cell interface at $x_{i+\frac{1}{2}}$ are only presented here. The three third order interpolation formula of variable $Q$ is given by

[TABLE]

where ${\bar{Q}}_{i+\frac{1}{2}}^{(k)}$ are approximated values at cell interfaces from different sub-stencils and $Q_{i}$ are the values at cell centers. In WENO literature, the variable $Q$ can either be fluxes [26, 40], conservative variables [41], primitive variables or characteristic variables [42]. In this current work the conservative variables $(Q)$, which are also the fluxes $(E)$, are directly interpolated. The three third order upwind approximation polynomials in Eqn 35, are chosen dynamically by a nonlinear convex combination which adapts either to a higher order approximation in smooth regions of the solution, or to a lower-order spatial discretization that avoids interpolation across discontinuities and provides the necessary numerical dissipation for shock capturing. The fifth order WENO scheme can be expressed as

[TABLE]

where $\omega_{k}$ are the nonlinear weights which are given by,

[TABLE]

where $m$, $\gamma_{k}$ and $\beta_{k}$ are a positive integer, ideal linear weights, and smoothness indicators, respectively. $\epsilon=10^{-6}$ is a small constant to prevent division by zero. The non-linear weights of the convex combination are based on local smoothness indicators $\beta_{k}$, which measure the sum of the normalized squares of the scaled $L^{2}$ norms of all derivatives of the lower order polynomials. The basic weighting strategy is to assign small weights to those lower order polynomials whose underlying stencils contain discontinuities so that an essentially non-oscillatory solution is obtained. The traditional smoothness indicators for fifth order upwind interpolation are given by Jiang, and Shu [16] denoted as WENO-JS.

[TABLE]

where $k$ is the polynomial degree of $p_{j}(x)$. Evaluating of each $k$, we obtain the following equations

[TABLE]

Borges et al. [19] proposed a new approach, denoted as WENO-5Z in this paper, for the nonlinear weights obtained by WENO-JS as they are known to lose accuracy at critical points and are also excessively dissipative in smooth regions. The improved non-linear weights are as follows:

[TABLE]

where the smoothness indicators $\beta_{k}$’s are the same as those given in Eqs. (39), $\epsilon=10^{-40}$, and $\tau_{5}$ is the smoothness indicator of the large stencil given by,

[TABLE]

The variable $p$ is used to tune the dispersive and dissipative properties of the scheme. It is reported by Borges et al. [19] that the scheme becomes more dissipative when $p$ is increased. In this paper, $p=2$ is employed in all test problems for diffusion and advection-diffusion equations and $p=1$ for magnetized electron fluids. Through spectral and approximate dispersion relation (ADR) analysis, Jia et al. [43] found that the anti-dissipation of the WENO-5Z scheme is less than that of the WENO-JS scheme. They also demonstrated that the WENO-5Z scheme is not only less dissipative and dispersive but also relatively more accurate and safer. These properties of the scheme were found to be useful for magnetized electron fluid simulations discussed later.

3.2 Linear and compact interpolation

For a smooth function, the fifth order WENO scheme is theoretically equivalent to the optimal fifth order linear upwind scheme, denoted as U-5E in this paper, that is, the nonlinear weights, $\omega_{k}$, are equal to the ideal linear weights, $\gamma_{k}$. The fifth order interpolation formulas for left and right interfaces are given by,

[TABLE]

[TABLE]

Similarly, third order extrapolation formulas, same as that of in Eq. (35), denoted as U-3E in this paper, are given by,

[TABLE]

Finally, we have also implemented the compact upwind interpolation, denoted as U-5C in this paper, given by equations (45). Compact schemes are a family of interpolation schemes which are implicit in space and therefore requires an inversion of a tridiagonal matrix. They are characterized by high spectral resolution and have significantly lower dispersion errors compared to that of non-compact schemes. Implementation of the boundary conditions is same as that of the non-compact schemes.

[TABLE]

Finally, the derivative $(\partial E/\partial x)$ is computed as follows for all the schemes:

[TABLE]

In general, the integrals of the fluxes are discretized using a high-order Gaussian quadrature with suitable Gaussian integration points over the faces of the control volume[44, 45] to achieve higher order accuracy in, third order or more, multidimensional finite volume method. In the present hyperbolic approach, the fluxes are no longer non-linear and are same as that of the conservative variables. Due to this, finite volume method can also obtain higher order accuracy by using point values and reconstruction of the fluxes at Gaussian quadrature points is not necessary. Therefore, we cannot distinguish the difference between cell-centered finite volume and finite difference schemes for the flux computation. In fact, we also computed the fluxes by WCNS in [46] and the results are identical. In this paper, we focus only on the standard WENO schemes and not a comparison of WENO and WCNS. The implementation of interpolation polynomials beyond $5^{th}$ order are straightforward and are not presented in this paper.

3.3 Total variation diminishing schemes

Several TVD schemes are proposed in the literature for shock-capturing and preventing Gibbs oscillations. Two of the popular approaches are considered here. First one is the van Leer’s standard TVD-MUSCL [15] reconstruction for the left and right interfaces, ${\bar{Q}^{L}}_{i+\frac{1}{2}}$ and ${\bar{Q}^{R}}_{i+\frac{1}{2}}$, are given by,

[TABLE]

where $\kappa$ is a free parameter which is set to 1/3 for the third-order limiter and

[TABLE]

where $\epsilon$ is a small value $10^{-16}$. In this paper we considered minmod and van Leer limiters which are given by

[TABLE]

Second approach is that of Cockburn and Shu [47] which can be considered as Generalized MUSCL and is popular in discontinuous Galerkin methods. For $k=\frac{1}{3}$, Eq. (47) will reduce to third order linear interpolation formula, Eq. (44), if limiters are not considered. Now we define,

[TABLE]

Then, by using the definition of minmod limiter we can obtained a generalized limiter for 3 points

[TABLE]

where the minmod function defined as

[TABLE]

3.4 Implementation of boundary conditions

In this section, the implementation of boundary conditions for the all the equations is described. Two different approaches are employed for the numerical boundary conditions, weak boundary condition, and ghost-cells. 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.

3.4.1 Weak boundary formulation

The Dirichlet boundary condition $u$ = $u_{0}$ at, say, $x=0$ can be implemented through the numerical flux by weak boundary condition as described by Nishikawa and Roe [48]. For a one-dimensional problem, the upwind flux for a is given by,

[TABLE]

[TABLE]

where ($u_{R}$, $p_{R}$) are given by a higher order interpolation from the interior of the domain, and the left state ($u_{L}$, $p_{L}$) is specified by the boundary condition:

[TABLE]

Since the value of gradient variable is not known we set $p_{L}$ = $p_{R}$. This approach is consistent with the characteristic condition at $x=0$. We can specify only one condition in hyperbolic approach only one wave enters the domain as shown in Fig. 3. For a Neumann boundary condition, say $p=p_{0}$, the value of $p$ is specified instead of $u$:

[TABLE]

3.4.2 Lagrange type extrapolation for ghost cells

In this approach, an additional ”ghost cells” are introduced by extending the physical domain. Unlike the weak boundary condition, the Dirichlet boundary condition is employed at the cell interface as shown in Fig. 4. The number of ghost-cells depends on the interior interpolation scheme and cell-interface to cell-center Lagrange extrapolation formula given by Eq. (57) is used to compute the values in ghost cells:

[TABLE]

where r is the order of the Lagrange polynomial and the constants $C_{rj}$ are obtained by,

[TABLE]

For Dirichlet boundary condition, where the value of $u$ is prescribed at the cell interface, say $u_{\frac{1}{2}}$, one can extrapolate the ghost cell values by using third order accurate polynomial given by Eq. (59).

[TABLE]

For the gradient variable, $p$, the values in ghost cells are extrapolated by using the interior points by extrapolation, cell-center to cell-center [49, 50]. Such extrapolation is also consistent with the outgoing characteristics at the boundary as shown in Fig. 3. For example, At the boundary, say $x=0$, the values of ghost cell $p_{0}$ can be approximated by 3rd order extrapolation given by Eq. (60),

[TABLE]

Lagrange extrapolation can also be approximated by a (s-k)th order Taylor expansion

[TABLE]

where $u_{R}^{(k)}$ is a (s-k)th order approximation of $\frac{\partial^{k}u}{\partial x^{k}}$ at x = 1. If $u(x)$ is smooth near the boundary, $u_{R}^{(k)}$ can be obtained by

[TABLE]

where $P_{s-1}(x)$ is a Lagrange polynomial of degree at most $s-1$.

3.4.3 WENO type extrapolation for ghost points

In situations where there are sharp gradients close to the boundary, the Lagrange extrapolation may result in severe oscillations. To overcome such difficulties a more robust WENO type extrapolation is also proposed by Tan and Shu [49, 20]. Analogous to the idea of WENO, three polynomials are constructed for the extrapolation given by,

[TABLE]

Using the standard WENO procedure, the equation 62 can be rewritten as,

[TABLE]

where $\omega_{r}$ are the typical nonlinear weights given by, $\omega_{r}=\frac{\alpha_{r}}{\sum_{r=0}^{2}\alpha_{s}}$ and $\alpha_{r}=\frac{d_{r}}{(\epsilon+\beta_{r})^{3}}$, where the linear weights $d_{r}$ are chosen as

[TABLE]

and the smoothness indicators $\beta_{r}$ are obtained by,

[TABLE]

Readers can refer [20] for the explicit expressions of smoothness indicators. In this paper $3^{rd}$ order WENO extrapolation is considered for all the simulations.

3.5 Time discretization

After discretizing the spatial derivative, the set of ordinary differential equations obtained are,

[TABLE]

where the operator $Res(Q)=E^{\prime}_{x}$ and $E^{\prime}_{x}$ is approximated by linear upwind interpolations or WENO interpolation. For time integration the following third order TVD Runge-Kutta method [16] is used

[TABLE]

A fourth order non-TVD Runge-Kutta scheme can also be employed for the computations. Numerical results obtained by Eq. (3.5) are only presented in this paper.

4 Numerical Results and Discussion

In this section, various test cases in one and two dimensions for diffusion equation, advection-diffusion equation, and electron fluid equations are tested by using the upwind schemes and the boundary conditions described in earlier section. All the simulations are carried out on a uniform Cartesian mesh, and TVD-RK is used for time integration for all the problems.

4.1 Diffusion equation

For the diffusion equation, the numerical solutions are also computed with the standard $2^{nd}$ order central scheme using successive over relaxation for comparison. For the linear upwind schemes, depending on the order of the interior scheme the corresponding $r^{th}$ order extrapolation polynomial is used for the numerical boundary conditions. Only Lagrange extrapolation and weak boundary procedures are implemented for the hyperbolic form of diffusion equation as WENO extrapolation is found to be unstable.

Example 4.1.1. To investigate the implementation and accuracy of the numerical schemes the following one-dimensional diffusion equation including a source term is considered with the domain size of $x$ $\in$ [0,1].

[TABLE]

where $\nu=1$ and Dirichlet boundary conditions, $u_{1}=2,u_{N}=1$, are considered. The exact steady state solution for $u$ is given by:

[TABLE]

where the constants and $A$ and $C$ are assumed to be 10.0 and 3.0 respectively.

The equivalent first order hyperbolic equation system can be written as follows:

[TABLE]

The conservative variables are ${\bf{Q}}=[\ u,\ p]$ and the flux vector is ${\mathbf{E}}=[-p,-\ u]$. The Jacobian matrix and the corresponding eigenvalues for the above system of equations is,

[TABLE]

The right and left eigenvectors can be used to compute the absolute flux Jacobian

[TABLE]

The simulations are conducted with grid refinements from $N$ = 24 to 384 by second-order central, and all the upwind schemes including TVD-MUSCL scheme. For upwind schemes, the numerical solution is computed by an explicit time-marching until the residuals are dropped below $10^{-12}$ in $L_{1}$ norm with a constant CFL $=0.65$. Table 1 shows the $L_{2}$ error for the velocity and the order of accuracy for the U-3E, U-5E,U-5C and WENO-5Z schemes respectively. We can see that design order of accuracy is obtained for all the schemes.

Table 2 shows the $L_{2}$ error for the gradient variable, $p$, and the order of accuracy for the U-3E, U-5E,U-5C and WENO-5Z schemes respectively. We can see that design order of accuracy is obtained for the gradient variable as well for all the schemes. $L_{2}$ error convergence results are shown in Fig. 5(a) and Fig. 5(b) for the solution and the gradient variables respectively.

Fig. 6 shows the solutions contours of various schemes in comparison with the analytical solution for both primary and gradient variables. Fig. 6(c) shows the comparison of ghost cell approach and weak boundary implementation, and we can observe that ghost cell approach represents the solution more accurately on coarse meshes.

Weak boundary implementation was unstable beyond $3^{rd}$ order accuracy whereas the ghost cell approach was found to be stable until 6th order accuracy. We note that the ghost cell approach may be complicated to implement on an unstructured mesh compared to the weak formulation and also the difference between these approaches is minimal on finer meshes. Based on these observations only ghost cell approach has been employed for the all the test cases in the next subsections. WENO extrapolation is not considered for diffusion equation as it is found to be unnecessary.

In Fig. 7(a), we can observe the difference between the shock-capturing schemes, TVD-MUSCL and Generalized-MUSCL. TVD-MUSCL formulation has not contaminated the solution and it is reduced to the linear third order scheme, regardless of the limiter. On the other hand, the Generalized-MUSCL approach resulted in unnecessary oscillations and is only first order accurate, shown in Fig. 7(b), and it may not be appropriate for the diffusion equation in hyperbolic form.

Fig. 8(a) shows the solution obtained by weighted essentially non-oscillatory schemes, WENO-JS, WENO-M and WENO-Z respectively. As expected, WENO-Z and WENO-M schemes have better accuracy than WENO-JS. WENO-M and WENO-Z gave similar results but it is well known that WENO-M scheme is more computationally intensive. Based on this analysis WENO-Z is considered for all the simulations in following sections. Order of accuracy for WENO schemes is shown in Fig. 8(b).

From this analysis, shock-capturing schemes can be considered unnecessary for diffusion equation in hyperbolic approach. TVD schemes contaminated the solution depending on the type of the scheme considered whereas the WENO scheme did not show any unnecessary oscillations for all the variations. It was also mentioned by Nishikawa [3] that no discontinuity capturing mechanism is required for the upwind formulation of the diffusion equation.

Example 4.1.2. First test case is the following two-dimensional diffusion equation considered by Nishikawa [3]

[TABLE]

with a spatial domain of $[0,1]\times[0,1]$, where $\nu$=1 and the following Dirichlet boundary conditions are considered:

[TABLE]

The exact steady state solution is given by,

[TABLE]

The simulations are conducted with grid refinements from 16$\times$16 to 256 $\times$ 256 by second-order central, and all the upwind schemes. For upwind schemes, the numerical solution is computed by an explicit time-marching until the residuals are dropped below $10^{-12}$ in $L_{1}$ norm with a constant CFL $=0.5$. The exact solution and numerical solution contours computed by the U-5E scheme are shown in Fig. 9(a) and Fig. 9(b). Computed values of $u$ for various schemes along the geometric center line along the horizontal axis are shown in 10(a). $L_{2}$ error convergence results are shown in Fig. 10(b) for the primary variable and we can observe that the design order of accuracy is obtained for all the schemes.

Table 3 shows the $L_{2}$ error for the primary variable, $u$, and the order of accuracy for the U-3E, U-5E,U-5C and WENO-5Z schemes respectively. We can see that design order of accuracy is obtained for the velocity for all the linear upwind schemes. For the WENO scheme, the 5th order boundary conditions are found to be unstable, and only 3rd order boundary conditions are considered. TVD-MUSCL has once again reduced to 3rd order linear scheme and is not discussed henceforth for diffusion problems.

Example 4.1.3 In this test case, the verification of Neumann boundary condition which is also implemented as Dirichlet boundary condition in the hyperbolic formulation for diffusion equation is considered. The following two-dimensional Laplace equation with the spatial domain of [0, 1] $\times$ [0, 1], Eq. (77), has been considered.

[TABLE]

which has the following Neumann and Dirichlet boundary conditions

[TABLE]

where $L=1$ and $M=1$ and the analytical solution is given by Eq. (79)

[TABLE]

The simulations are conducted with grid refinements from 16$\times$16 to 256 $\times$ 256 for the test case. For the upwind schemes, the numerical solution is again computed by an explicit time-marching until the residuals are dropped below $10^{-12}$ in $L_{1}$ norm with a constant CFL $=0.5$. Computed values of $\phi$ for various schemes along the geometric center line along the horizontal axis are shown in Fig. 11 for a 32 $\times$ 32 grid. Design accuracy is obtained for all the schemes, and the implementation of Neumann boundary condition is successfully verified through this test case. Unlike the previous test case, WENO scheme also shows 5th order accuracy. Table 4 show the $L_{2}$ norms for all the schemes and the advantage of high-order methods over lower-order methods can be seen as they need less number of computational cells to get a solution with the same accuracy. This advantage enables high-order methods to use coarse meshes, in comparison with the lower-order methods.

Example 4.1.4. For the previous test cases and the given numerical grids, we could not distinguish the numerical results by different methods. Hence, we have only presented the numerical results of these schemes. In this example, we consider the following Poisson equation, i.e., with a source term, given by Eq. (80)

[TABLE]

where the domain is [0, 1] $\times$ [0, 1] and $\phi=0$ at all the boundaries. The exact solution for this test case is given by

[TABLE]

In this test case, the advantage of 5th order schemes can be observed. The simulations are conducted with grid refinements from 16$\times$16 to 256 $\times$ 256 for the test case. For the upwind schemes, the numerical solution is deemed to have reached a steady state when the residuals are dropped below $10^{-12}$ in $L_{1}$ norm and constant CFL $=0.5$ is used. In Fig. 12(a) the computed results at geometric center on a 32 $\times$ 32 grid is shown and we can observe that the third order explicit scheme has deviated considerably from the exact solution in comparison with the other schemes. We can see from the $L_{2}$ norms shown in the Fig. 12(b) that even though third order accuracy is obtained for the upwind scheme, U-3E, the solution is inferior to standard second order case on coarse meshes. Such results are also observed in edge based methods proposed by Nishikawa in [Ref. [4], Figs. 4 and 5]. The source term may have a significant effect on the solution accuracy for the upwind formulation. Table 5 shows the $L_{2}$ error for the primary variable, $\phi$, and the order of accuracy for the U-3E, U-5E, U-5C and WENO-5Z schemes respectively. Compact schemes show better accuracy for problems with source term with increasing grid size in comparison with the explicit scheme.

Example 4.1.5. In this test case the following anisotropic diffusion equation identified by Kuzmin et al. [51], with a spatial domain $(x,y)\in[0,1]\times[0,1]$, has been considered.

[TABLE]

where the diffusion tensor and source term are given by

[TABLE]

The exact steady state solution for this test case is given by,

[TABLE]

, which is also imposed as Dirichlet boundary conditions. Again, grid refinement study is carried out for all the upwind schemes using explicit time stepping. The solution contours obtained on a 32 $\times$ 32 grid shown in Figs. 13(a) and 13(b) are indistinguishable and also unlike the slope limiting approach of Kuzmin et al. the current scheme does not pollute the smooth extrema and also satisfies the discrete maximum principle. Table 6 shows the design order of accuracy for the primary variable, $u$, for the U-3E, U-5E,U-5C and WENO-5Z schemes respectively. Once again, compact schemes show better accuracy for problems with increasing grid size in comparison with the explicit scheme. Computed values at geometric center and $L_{2}$ convergence errors for various schemes are shown in Figs. 14(a) and 14(b) respectively. Convergence rates for anisotropic diffusion problem by using explicit time stepping are extremely slow and can be significantly improved by implicit time marching approach [52].

4.2 Advection-Diffusion equation

For the advection-diffusion equation, the numerical solutions are also computed on a uniform cartesian mesh. For the linear upwind schemes, depending on the order of the interior interpolation the corresponding $r^{th}$ order extrapolation polynomial is used for the ghost boundary conditions. For WENO scheme both Lagrange and WENO extrapolations are considered for boundary conditions.

Example 4.2.1. The following one dimensional-advection diffusion equation is considered with the spatial domain of [0,1] and $u(0)=0$ and $u(1)=1$

[TABLE]

where, $s(x)=\pi Re[a\cos(\pi x)+\pi\nu sin(\pi x)$ and the Reynolds number, $Re$ = $\frac{a}{\nu}$. The exact steady state solution is given by [6]

[TABLE]

The exact solution is a smooth sine curve for low Reynolds numbers, which are diffusion dominant, and develops a sharp gradient close to the boundary for high Reynolds numbers when advection becomes the dominant effect. The main objective of this test case to verify the design accuracy and implementation of the numerical schemes and also to show the advantage of WENO scheme in capturing the sharp gradients.

The simulations are conducted on a uniform grid with grid refinements from $N$ = 8 to 256 by all the upwind schemes. The numerical solution is computed by an explicit time-marching until the residuals are dropped below $10^{-10}$ in $L_{1}$ norm. A constant value of CFL is considered depending on the order of the scheme. Simulations are carried out for three Reynolds numbers: $Re=0.1,1,10$. We can see that the design order of accuracy is obtained for all the schemes for all the Reynolds numbers in Fig. 16. The 5th order compact scheme is one order of magnitude more than the design accuracy for $Re=0.1$. Fig. 15 shows the solutions contours of various schemes in comparison with the analytical solution for the primary variable, $u$. From Fig. 15(b), it is obvious to notice that the Generalized MUSCL scheme can pollute the smooth profile whereas TVD-MUSCL reduces to linear 3rd order scheme. For Re = 1 there were unnecessary oscillations in the profile, and for the Reynolds number of 10, even though there are no oscillations, the sharp gradient has been completely cut-off. Numerical results obtained by linear upwind schemes were, in fact, better than the results obtained by Generalized MUSCL scheme.

To obtain the solution without oscillations in the case of high-Re, extremely fine grids are required to meet the well-known requirement on the mesh Reynolds-number, i.e., $Re<2$, as explained by Nishikawa [6]. Computations are carried out by Nishikawa on substantially coarser grids by grid stretching to meet the criteria. Mazaheri and Nishikawa [8] had merely refined the grid to satisfy the mesh Reynolds-number requirement by using a highly efficient implicit solver. In this test case, we incorporated the non-oscillatory schemes, WENO, to capture the sharp gradients without oscillations on uniform meshes which might not satisfy the mesh Reynolds-number requirement. Simulations carried out for the Reynolds number of 500 on a uniform grid size of 256 has been shown in Fig. 17. Generalized MUSCL scheme did not show any oscillations but completely cut-off the sharp gradient and thereby failing to capture it. On the other hand, TVD-MUSCL has reduced oscillations compared to that of linear 3rd order scheme. Oscillations are also observed in both 5th order schemes. However, the amplitude of the oscillations is much smaller for the compact scheme than the explicit scheme as expected. Finally, sharp gradient developed at $x=1$ is easily captured by WENO scheme without spurious oscillations as shown in Fig. 17(d). Also, it is important to note that the linear schemes can also give oscillation free solutions on dense grids by satisfying the mesh-Reynolds number criteria. These important inferences were found to be useful for the simulation of magnetized electrons which can develop sharp gradients with increasing strength of magnetic confinement.

Example 4.2.2. The two-dimensional advection-diffusion equation is solved by split flux approach for two test cases.

[TABLE]

The first test case is an advection-dominated problem considered by Nishikawa and Liu [53] which has the following exact solution,

[TABLE]

where $\xi$ = ax + by, $\eta$ = bx - ay and the viscosity, $\nu=0.01$. The advection vector $(a,b)$ is given by $(1.0,0.0)$. This exact solution is given by Eq. (88) is very smooth and therefore can be used to verify the accuracy of numerical schemes. 2. 2.

For the second test case the following exact solution is considered,

[TABLE]

where the viscosity coefficient is $\nu=0.1$, and the advection vector $(a,b)$ is given by $(7.0,4.0)$ respectively. This test case was originally proposed by Nishikawa and Roe in [54], where both advective and the diffusive terms are equally important which may exist in the middle of the boundary layer.

Computations are carried out with grid refinements from 8$\times$8 to 128 $\times$ 128 for both the test cases. For all the numerical schemes, the simulations are carried out until the residuals are dropped below $10^{-14}$ in $L_{1}$ norm and CFL is dependent on the numerical scheme used, which is in between from 0.3 to 0.7. Both WENO extrapolation denoted as WENO-5Z-W, and Lagrange extrapolation, WENO-5Z-L, are implemented for the numerical boundary conditions for the WENO scheme, which are 3rd order accurate.

Figs. 18(a) and 18(b) show the velocity contours obtained by the upwind scheme, U-5E, for both the test cases on a 48 $\times$ 48 grid. Design accuracy is obtained for all the linear upwind schemes, shown in Fig. 19(a) and the implementation of boundary conditions for advection dominated problem are verified through this test case. First the first test condition, WENO scheme is third accurate for both WENO and Lagrange type extrapolations which is unlike the one-dimensional test case where the WENO scheme was fifth order accurate. The second test case was found to be more challenging than the advection-dominant flow. Design accuracy is obtained for all the linear upwind schemes, shown in Fig. 19(b). For the WENO scheme, the WENO extrapolation technique for the boundary conditions was observed to be first order accurate whereas the Lagrange extrapolation is 3rd order accurate similar to the advection-dominated problem. Table 7 show the $L_{2}$ norms for both the test cases for both WENO implementations.

Example 4.2.3. In this test case, we solve the steady state problem of the two-dimensional advection-diffusion equation, with boundary layers along x = 1 and y = 1, considered by Chou and Shu [55] which has the following exact solution given by

[TABLE]

The viscosity coefficient $\nu$ is taken as 0.05, and the exact solution is imposed on the boundaries. Linear upwind schemes are not considered for this problem as they will produce spurious oscillations. The numerical results are shown in Table 8 indicate that the third order accuracy is obtained for WENO-5Z-L whereas the accuracy of WENO-5Z-W is reduced to first order. This test case demonstrates that the WENO scheme can resolve the boundary layers by using first-order hyperbolic approach.

4.3 Magnetized electron fluid simulation

Magnetized electron fluid simulations can have all the features that are discussed in earlier sections like sharp gradients in the flow field, strongly diffusion dominated flow with anisotropic diffusion and also have the convection aspect of the flow. A key distinction between our present methodology and the earlier efforts based on first order upwind and third order accurate monotonicity preserving TVD discretization is in the implementation of high-order finite-volume WENO discretization. The gain in accuracy and reduction in numerical oscillations are significantly more pronounced for the WENO formulations when compared with TVD discretization. Before we discuss the results of the magnetized electron fluids we briefly summarized the inferences from the diffusion and advection-equations to better understand the issues faced:

• •

Both diffusion and advection-diffusion equations, linear upwind schemes are sufficient for problems with smooth solutions.

• •

WENO schemes are more appropriate for capturing the sharp gradients than TVD type scheme.

• •

Ghost cell approach is preferable over the weak boundary. WENO extrapolation can reduce to first order in some instances and yet provide robust results.

4.3.1 Steady state results and analysis

We consider the test calculation considered by Kawashima et al. [1] with uniformly angled Magnetic lines of force at $45^{\circ}$ from the vertical, shown in Fig. 20. For a Cartesian mesh, the condition of 45*∘* magnetic lines of force can give significant false diffusion [56]. False diffusion will occur due to the oblique flow direction and non-zero gradient in the direction normal to the flow. As a consequence of false diffusion, nonphysical local extrema can occur in regions in which the gradients of the solution are steep and not aligned with the orientation of grid, and the discretization method will be unable to capture them properly. Results were obtained for three different values of magnetic confinements, ${\mu}_{||}/{\mu}_{\perp}$=100, 500 and 1000, which are uniform 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. (91).

[TABLE]

Unlike the test cases of diffusion and advection-diffusion where the numerical schemes are compared with the analytical solutions, the linear and WENO upwind schemes are compared with the Magnetic Field Aligned Mesh (MFAM). MFAM is described by Mikellides et al. [57] and Kawashima et al. [1], which eliminates the false diffusion by aligning the mesh with the magnetic field. The simulations are conducted for both MFAM and the upwind schemes with grid refinements of $NX,NY$ = [96,48] to [768,384] for all the values of ${\mu}_{||}/{\mu}_{\perp}$. For the upwind methods, the numerical solution is computed by explicit time marching with a constant CFL of $0.32$. Regardless of the interior scheme, 3rd order accurate boundary conditions are considered for all the upwind schemes for stable computation. The steady-state solution is deemed to be reached if the $L_{1}$ norms are below $10^{-12}$. As far the MFAM approach, successive over-relaxation is used to numerically relax the solution until the $L_{2}$ norms are below $10^{-10}$ by second-order central discretization.

Figs. 21 and 22 show the velocity streamlines computed by using different schemes for mobility ratios ${\mu}_{||}/{\mu}_{\perp}$=10, 500 and 1000. Velocity streamlines computed by MFAM are shown in Figs. 22(a), 22(b) and 22(c) respectively and we can see that as ${\mu}_{||}/{\mu}_{\perp}$ is increasing the plasma is more confined to the center. For ${\mu}_{||}/{\mu}_{\perp}=100$ all the linear schemes and WENO has similar results compared to that of the MFAM whereas the Generalized MUSCL approach is showing unphysical streamlines even for such small anisotropies. From the Figs. 21(c) and 22(f), we can see that there are significant unphysical vortices in Upwind-3E for ${\mu}_{||}/{\mu}_{\perp}=1000$. Results obtained by TVD-MUSCL are same as that of Upwind-3E scheme, similar to the results observed for isotropic diffusion problems. However, the fifth order schemes, Upwind-5C, Upwind-5E, and WENO-5Z-L, were able to reduce the unphysical oscillations significantly and the difference is between them is very small. The reason for the unphysical oscillations can be understood from the velocity contour plots shown in Figs. 23 and 24.

We see from a series of Figs. $\ref{fig:3E-100-V}-\ref{fig:WENO-1000-V}$, that there can be serious numerical oscillations in the solutions of the flow field with increasingly strong magnetic confinement. The main reason is not only the anisotropic diffusion itself but also because of the sharp gradients that are not aligned with the grid. In Fig. 24(c), MFAM approach do not show such behavior through an alignment of field lines with the grid and thereby captures the sharp gradient without numerical oscillations. If the angle between the field lines and the mesh is aligned, even the central scheme could also be easily implemented. As the gradient becomes more and more skewed, it diverges much more from the grid lines, thereby generating larger oscillations, as seen in Fig. 23(c). WENO-5Z was able to capture the gradients without significant oscillations as the disturbances in the y-velocity are minimal when the solution has reached the steady state. Steady state solutions obtained for linear upwind schemes indicate that the disturbances behind the strong gradient have polluted the entire domain but similar to the boundary layer problem the solutions may be reduced with finer grids.

The oscillations due to sharp gradients can be drastically reduced on finer meshes by WENO-5Z. Even though MFAM approach did not show any spurious oscillations, such alignment may not be practical in many simulations with multiple sharp gradients. Furthermore, the boundary conditions and coding can be challenging to implement. Based on these simulations we can say that the non-oscillatory approach of WENO can be a reasonable and viable alternative.

Similar patterns are also observed in potential distributions shown in Fig. 25. The maximum principle is violated as a consequence of these spurious oscillations as we can see the values of space potential are more than the theoretical bounds defined by the boundary conditions. Although the values are relatively smaller, it should be noted that the values are non-dimensionalized and the actual error can be much higher. On coarse grids, the entire domain can be polluted with unphysical extrema. The grid size of 192 $\times$ 96 is chosen in Fig. 25 is to reiterate the fact that the spurious oscillations can be reduced by the linear upwind schemes whereas they remain the same for the Generalized-MUSCL. Spurious oscillations in potential are comparatively smaller in the compact scheme than that of explicit schemes.

As explained earlier, the steady state is considered to be reached when the $L_{1}$ norms attain negligibly small values in the order of 10*-12*, that are defined by Equation (92)

[TABLE]

where $Q$ denotes each of the conservative variables $(\tilde{\phi},\tilde{u}_{\rm x},\tilde{u}_{\rm y})$. The convergence history of space potential for various schemes is shown in Fig. 26 when the test problem is calculated on a 192 $\times$ 96 grid. The residuals are monotonically decreasing and reached the limit as mentioned earlier for all the linear upwind schemes and WENO scheme. For the Generalized MUSCL scheme, the residual does not converge beyond four orders of magnitude.

4.3.2 Effect of boundary conditions

Simulations carried out by WENO-5Z along with Lagrange extrapolation, WENO-5Z-L, for ghost cells are free of spurious values in space potential if ${\mu}_{||}/{\mu}_{\perp}\leq 500$. But with increasingly higher anisotropies small overshoots and undershoots, similar to those that are observed in linear upwind schemes, are also found in the cells next to the boundary. By keeping the objective of the uniform grid, WENO extrapolation method, WENO-5Z-W, is used for computing the values in ghost-cells.

The WENO-5Z-W scheme was able to adaptively reduce to lower accurate order boundaries and thereby prevented the unphysical oscillations in the flow field, as shown in Fig. 27. WENO extrapolation approach was also tested for smaller anisotropies, and the results indicated that the performances are similar to that of WENO-5Z-L and no unphysical extrema are generated. The minimum grid size required for WENO-5Z-W scheme to prevent unphysical extrema in space potential distribution is 16 $\times$ 16 for ${\mu}_{||}/{\mu}_{\perp}$ $\leq 500$ and 96 $\times$ 96 if ${\mu}_{||}/{\mu}_{\perp}$ $>500$.

Test calculations are also carried out for the cases when the magnetic field lines are aligned at an angle $60^{\circ}$ from the vertical and even for complicated shapes constructed using the magnetic stream function. These results further confirm the fact that the WENO scheme has better capabilities in capturing the sharp gradients without spurious oscillations even for complex magnetic field shapes.

4.3.3 Mesh convergence and accuracy for electron fluids

The computational accuracy of anisotropic diffusion by HES approach is compared with that of MFAM by using the transverse electron flux which is defined as the total electron flux from right to left boundary given by Eq. (93),

[TABLE]

where $\Omega_{L}$ is the left boundary of the calculation field. The grid convergence of the calculated electron transverse flux is evaluated only for WENO-5Z-W and MFAM schemes as for the other schemes the solution is polluted due to spurious oscillations and might not be appropriate. Table 9 shows the difference between computed transverse electron fluxes by both the approaches. The error reduces from 1.12$\%$ to 0$\%$ on mesh refinement.

Another criteria considered is the relative numerical error, $|l_{2}|$ error, which is defined as

[TABLE]

where Q is each of the conservative variable. The relative errors are calculated in comparison with the field aligned mesh results. The Table 10 show the relative $l_{2}$ errors computed for field aligned mesh and WENO-5Z-W and the results show that the solution seems to improve with finer grid resolution.

5 Conclusions

Higher order linear and non-linear schemes (WENO) are proposed for the simulation of magnetized electrons in quasi-neutral plasmas by the hyperbolic method. We first implemented the approach for diffusion equation, which has similar upwind structure, in hyperbolic form to verify the design accuracy and the implementation of boundary conditions. Then the schemes are implemented to simulate advection-diffusion equation to capture the sharp gradients in boundary layer type problems without spurious oscillations. Finally, the schemes are applied to anisotropic diffusion electron fluid equations to reduce the spurious oscillations due to sharp gradients significantly. The current methodology can be extended to the energy equation for the complete simulation of electron flow in magnetized plasmas. Implicit time stepping can substantially accelerate the steady-state convergence of all the problems considered here and will be discussed in future work. The critical findings of the paper are summarized as follows

High-order and high-resolution methods are implemented successfully for the diffusion equation in the hyperbolic form on uniform meshes. Design order of accuracy is obtained for all the schemes for all the test problems considered. 2. 2.

Shock-capturing schemes are found to be unnecessary for simple diffusion equation, and through various test cases, the inapplicability of certain TVD schemes is explained. Even though the WENO scheme is a shock-capturing scheme the steady-state solutions are not contaminated and are similar to the linear upwind schemes. 3. 3.

Ghost cell approach is found to be more accurate and stable than the weak boundary implementation, especially on coarse meshes. Linear upwind schemes are consistently stable with the corresponding higher order boundary conditions whereas WENO schemes were stable only with $3^{rd}$ order boundary conditions due to their inherent non-linearity. 4. 4.

Weighted essentially non-oscillatory schemes are implemented to capture the sharp gradients without spurious oscillations for boundary layer type problem for advection-diffusion equation. For smooth solutions design order of accuracy is obtained for all the higher order methods. 5. 5.

The difference between $5^{th}$ order explicit and compact schemes is indistinguishable for diffusion equation. As for the advection-diffusion equation, the compact schemes are slightly more accurate than explicit schemes. The advantage of compact schemes may be more pronounced if the current approach is extended to time-dependent problems and also to hyperbolic compressible and incompressible Navier-Stokes equations proposed by Nishikawa [58] and are currently being studied. 6. 6.

Significant improvement is achieved by utilizing higher order linear upwind methods and WENO schemes for the simulation of magnetized electron fluids as the results are much closer to the field-aligned mesh. For small anisotropies, linear upwind schemes as well as TVD-MUSCL can be appropriate. 7. 7.

Similar to the boundary layer problem WENO approach is found to be more suitable and robust to reduce the spurious oscillations associated with the sharp gradients with increasing anisotropic diffusion in electron fluid equations. Boundary conditions based on WENO extrapolation are found more appropriate to prevent unphysical extrema for anisotropies higher than 500. 8. 8.

For very high anisotropic diffusion problems, ${\mu}_{||}/{\mu}_{\perp}$ of the order $10^{4}-10^{9}$ that can be seen in practical applications like tokamaks and space propulsion devices, all the schemes would result in spurious oscillations on coarse meshes, say 96 $\times$ 96, and reducing such oscillations is the subject of our future work.