Asymptotic-Preserving scheme for a strongly anisotropic vorticity equation arising in fusion plasma modelling

TL;DR

This paper develops an efficient asymptotic-preserving numerical scheme to accurately simulate the evolution of electric potential in the highly anisotropic peripheral region of tokamak plasmas, crucial for fusion reactor performance.

Contribution

It introduces a novel asymptotic-preserving method tailored for strongly anisotropic vorticity equations in plasma physics, addressing computational challenges without high costs.

Findings

Successfully handles strong anisotropy and non-linear boundary conditions

Maintains accuracy across different asymptotic regimes

Reduces computational costs compared to traditional methods

Abstract

The electric potential is an essential quantity for the confinement process of tokamak plasmas, with important impact on the performances of fusion reactors. Understanding its evolution in the peripheral region - the part of the plasma interacting with the wall of the device - is of crucial importance, since it governs the boundary conditions for the burning core plasma. The aim of the present paper is to study numerically the evolution of the electric potential in this peripheral plasma region. In particular, we are interested in introducing an efficient Asymptotic-Preserving numerical scheme capable to cope with the strong anisotropy of the problem as well as the non-linear boundary conditions, and this with no huge computational costs. This work constitutes the numerical follow-up of the more mathematical paper by C. Negulescu, A. Nouri, Ph. Ghendrih, Y. Sarazin, "Existence and uniqueness of the electric potential profile in the edge of tokamak plasmas when constrained by the plasma-wall boundary physics".

Tables

Table 1. Table 1. Model parameters and geometrical configuration of Case (M) and Case (P) .

	model parameters			geometric configuration
case study	$\eta$	$\Gamma$	$\nu$	$a$	$b$	$L_z$	$l$	$L_r$
Case (M)	$[0,\mathrm{\hspace{0.17em}1}]$	$1$	$1$	$1$	$2$	$3$	$1$	$2$
Case (P)	$[0,\mathrm{\hspace{0.17em}1}]$	$3$	${10}^{-3}$	$1$	$19$	$20$	$1$	$2$

Table 2. Table 2. Case (M). L 2 superscript 𝐿 2 L^{2} -norm of the error of the numerical solution with respect to the exact solution, ∥ ϕ η − ϕ e ​ x ( M ) ∥ 2 subscript delimited-∥∥ superscript italic-ϕ 𝜂 subscript superscript italic-ϕ 𝑀 𝑒 𝑥 2 \left\lVert\phi^{\eta}-\phi^{(M)}_{ex}\right\rVert_{2} , and estimated order of convergence, p 𝑝 p , of the ( A ​ P ) η subscript 𝐴 𝑃 𝜂 \left(AP\right)_{\eta} and ( S ​ P ) η subscript 𝑆 𝑃 𝜂 \left(SP\right)_{\eta} schemes, for several grid sizes ( Δ ​ z = Δ ​ r = h Δ 𝑧 Δ 𝑟 ℎ \Delta z=\Delta r=h ) and η = 0 ,  10 − 14 ,  10 − 8 ,  1 𝜂 0 superscript 10 14 superscript 10 8 1 \eta=0,\,10^{-14},\,10^{-8},\,1 . The symbol ✗ denotes the lack of a numerical solution due to blow-up of the numerical scheme.

			$\eta =0$		$\eta ={10}^{-14}$
			AP scheme		SP scheme		AP scheme
$h$	$M_z\times M_r$	$M$	$L^2$-norm error	$p$	$L^2$-norm error	$p$	$L^2$-norm error	$p$
$1/4$	$12\times 9$	$80$	$8.0960\times {10}^{-1}$	$2.258$	$1.3123\times {10}^{-1}$	$2.119$	$8.0960\times {10}^{-1}$	$2.258$
$1/8$	$24\times 17$	$288$	$1.6925\times {10}^{-1}$	$2.051$	$3.0213\times {10}^{-1}$	$-$	$1.6925\times {10}^{-1}$	$2.051$
$1/16$	$48\times 33$	$\mathrm{1\hspace{0.17em}088}$	$4.0852\times {10}^{-2}$	$2.013$	✗	$-$	$4.0852\times {10}^{-2}$	$2.013$
$1/32$	$96\times 65$	$\mathrm{4\hspace{0.17em}224}$	$1.0123\times {10}^{-2}$	$2.003$	✗	$-$	$1.0123\times {10}^{-2}$	$2.003$
$1/64$	$192\times 129$	$\mathrm{16\hspace{0.17em}640}$	$2.5251\times {10}^{-3}$	$1.995$	✗	$-$	$2.5252\times {10}^{-3}$	$1.989$
$1/128$	$384\times 257$	$\mathrm{66\hspace{0.17em}048}$	$6.3342\times {10}^{-4}$	$1.999$	✗	$-$	$6.4528\times {10}^{-4}$	$1.976$
$1/256$	$768\times 513$	$\mathrm{263\hspace{0.17em}168}$	$1.5852\times {10}^{-4}$	$-$	✗	$-$	$1.6403\times {10}^{-4}$	$-$

Equations

$$-{\partial }_t{\partial }_r^2\varphi -\frac{1}{\eta }{\partial }_z^2\varphi +\nu {\partial }_r^4\varphi =\mathcal{S}-\frac{1}{\eta }{\partial }_z\mathcal{F},t\ge 0,(r,z)\in \Omega ,$$

$${\partial }_r\varphi (0,r,z)={\partial }_r{\varphi }_{in}(r,z),(r,z)\in \Omega ,$$

$${\partial }_r\varphi (t,r,z)={\partial }_r^3\varphi (t,r,z)=0,t\ge 0,(r,z)\in \Sigma ,$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{z}\phi(t,r,a)=\eta(1-e^{\Lambda-\phi(t,r,a)})+{\mathcal{F}}(t,r,a)\,,\quad t\geq 0\,,\quad(r,z)\in\Gamma_{a}\,,\\[8.53581pt] \displaystyle\partial_{z}\phi(t,r,b)=-\eta(1-e^{\Lambda-\phi(t,r,b)})+{\mathcal{F}}(t,r,b)\,,\quad t\geq 0\,,\quad(r,z)\in\Gamma_{b}\,.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}N+\nabla\cdot(N\,{\mathbf{u}_{e}})-\nabla\cdot(D_{N}\,\nabla_{\perp}N)=S_{N}\,,\\[8.53581pt] \displaystyle\partial_{t}\Gamma+\nabla\cdot(\Gamma\,{\mathbf{u}_{i}})-\nabla\cdot(D_{\Gamma}\,\nabla_{\perp}\Gamma)=-\nabla_{||}P\,,\\[8.53581pt] \displaystyle\partial_{t}W+\nabla\cdot(W\,{\mathbf{u}_{i}})-\nabla\cdot(D_{W}\,\nabla_{\perp}W)=\nabla\cdot[N({\mathbf{u}^{i}_{\nabla B}}-{\mathbf{u}^{e}_{\nabla B}})]+j_{||}\ {\mathbf{b}}\,,\end{array}\right.$$

$$W:=\nabla \cdot \left[\frac{1}{|\mathbf{B}{|}^2}({\nabla }_{\perp }\varphi +\frac{1}{N}{\nabla }_{\perp }N)\right],j_{||}:=-\frac{1}{{\eta }_{||}}{\nabla }_{||}\varphi +\frac{1}{{\eta }_{||}N}{\nabla }_{||}N.$$

$${\mathbf{u}}_{\mathbf{E}}:=\frac{\mathbf{E}\times \mathbf{B}}{|\mathbf{B}{|}^2}{\mathbf{u}}_{\nabla B}^{i,e}:=\pm \frac{2T_{i,e}}{e|\mathbf{B}|}\frac{(\mathbf{B}\times \nabla \mathbf{B})}{|\mathbf{B}{|}^2},$$

$$(R)\,\,\,\left\{\begin{array}[]{l}\displaystyle-\partial^{2}_{z}\phi=-\partial_{z}{\mathcal{F}}\,,\quad t\geq 0\,,\quad(r,z)\in\Omega\,,\\[8.53581pt] \displaystyle\partial_{z}\phi(t,r,a)={\mathcal{F}}(t,r,a)\,,\quad t\geq 0\,,\quad(r,z)\in\Gamma_{a}\,,\\[8.53581pt] \displaystyle\partial_{z}\phi(t,r,b)={\mathcal{F}}(t,r,b)\,,\quad t\geq 0\,,\quad(r,z)\in\Gamma_{b}\,,\end{array}\right.$$

$${\varphi }^{\eta }=p^{\eta }+\eta q^{\eta },\forall \eta >0,$$

$$\left\{\begin{array}[]{l}\displaystyle-\partial_{z}^{2}p^{\eta}=-\partial_{z}{\mathcal{F}}\,,\quad t\geq 0\,,\quad(r,z)\in\Omega\,,\\[8.53581pt] \displaystyle\partial_{z}p^{\eta}={\mathcal{F}}\quad\textrm{on}\quad\Gamma_{a}\cup\Gamma_{b}\,,\end{array}\right.$$

$$q_{|{\Gamma }_q}^{\eta }\equiv q_{\star }\text{on}{\Gamma }_q:=\{(r,z)\in \Omega /z=\frac{a+b}{2},r\in [0,L]\}\Rightarrow p_{|{\Gamma }_q}^{\eta }={\varphi }_{|{\Gamma }_q}^{\eta }-\eta q_{\star }.$$

$$(AP)_{\eta}\,\,\,\left\{\begin{array}[]{l}\displaystyle-\partial_{t}\partial^{2}_{r}\phi^{\eta}-\partial^{2}_{z}q^{\eta}+\nu\partial^{4}_{r}\phi^{\eta}={\mathcal{S}}\,,\quad t\geq 0\,,\quad(r,z)\in\Omega\,,\\[8.53581pt] \displaystyle-\partial^{2}_{z}\phi^{\eta}=-\eta\partial^{2}_{z}q^{\eta}-\partial_{z}{\mathcal{F}}\,,\quad t\geq 0\,,\quad(r,z)\in\Omega\,,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle q^{\eta}_{|\Gamma_{q}}\equiv q_{\star}\\[8.53581pt] \displaystyle\partial_{z}q^{\eta}_{|\Gamma_{a}}=(1-e^{\Lambda-\phi^{\eta}(t,r,a)})\\[8.53581pt] \displaystyle\partial_{z}q^{\eta}_{|\Gamma_{b}}=-(1-e^{\Lambda-\phi^{\eta}(t,r,b)})\\[8.53581pt] \displaystyle\partial_{z}\phi^{\eta}_{|\Gamma_{a}}=\eta(1-e^{\Lambda-\phi^{\eta}(t,r,a)})+{\mathcal{F}}(t,r,a)\\[8.53581pt] \displaystyle\partial_{z}\phi^{\eta}_{|\Gamma_{b}}=-\eta(1-e^{\Lambda-\phi^{\eta}(t,r,b)})+{\mathcal{F}}(t,r,b)\,.\end{array}\right.$$

$$(AP)_{0}\,\,\,\left\{\begin{array}[]{l}\displaystyle-\partial_{t}\partial^{2}_{r}\phi^{0}-\partial^{2}_{z}q^{0}+\nu\partial^{4}_{r}\phi^{0}={\mathcal{S}}\,,\quad t\geq 0\,,\quad(r,z)\in\Omega\,,\\[8.53581pt] \displaystyle-\partial^{2}_{z}\phi^{0}=-\partial_{z}{\mathcal{F}}\,,\quad t\geq 0\,,\quad(r,z)\in\Omega\,,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle q^{0}_{|\Gamma_{q}}\equiv q_{\star}\\[8.53581pt] \displaystyle\partial_{z}q^{0}_{|\Gamma_{a}}=(1-e^{\Lambda-\phi^{0}(t,r,a)})\\[8.53581pt] \displaystyle\partial_{z}q^{0}_{|\Gamma_{b}}=-(1-e^{\Lambda-\phi^{0}(t,r,b)})\\[8.53581pt] \displaystyle\partial_{z}\phi^{0}_{|\Gamma_{a}}={\mathcal{F}}(t,r,a)\\[8.53581pt] \displaystyle\partial_{z}\phi^{0}_{|\Gamma_{b}}={\mathcal{F}}(t,r,b)\,,\end{array}\right.$$

$$t_k:=k\Delta t,k=0,\dots ,N_t,\Delta t:=T/N_t.$$

$$-{\partial }_r^2\left(\frac{{\varphi }^{n+1}-{\varphi }^n}{\Delta t}\right)-\frac{1}{\eta }{\partial }_z^2{\varphi }^{n+1}+\nu {\partial }_r^4{\varphi }^{n+1}={\mathcal{S}}^{n+1}-\frac{1}{\eta }{\partial }_z{\mathcal{F}}^{n+1},$$

$$-{\partial }_r^2{\varphi }^{n+1}-\frac{\Delta t}{\eta }{\partial }_z^2{\varphi }^{n+1}+\nu \Delta t{\partial }_r^4{\varphi }^{n+1}=\Delta t{\mathcal{S}}^{n+1}-\frac{\Delta t}{\eta }{\partial }_z{\mathcal{F}}^{n+1}-{\partial }_r^2{\varphi }^n,$$

$$\{\begin{array}{rl}& {\partial }_r{\varphi }^{n+1}={\partial }_r^3{\varphi }^{n+1}=0\text{on}{\Sigma }_0\cup {\Sigma }_l\cup {\Sigma }_{L_r};\\ & \text{periodic boundary conditions for ϕn+1 on Γ0∪ΓLz};\\ & {\partial }_z{\varphi }_a^{n+1}=\eta \left(1-e^{\Lambda -{\varphi }_a^{n+1}}\right)+{\mathcal{F}}_a^{n+1}\text{on}{\Gamma }_a;\\ & {\partial }_z{\varphi }_b^{n+1}=-\eta \left(1-e^{\Lambda -{\varphi }_b^{n+1}}\right)+{\mathcal{F}}_b^{n+1}\text{on}{\Gamma }_b;\end{array}$$

$$\{\begin{array}{rl}& -{\partial }_r^2{\varphi }^{n+1}-\Delta t{\partial }_z^2{q}^{n+1}+\nu \Delta t{\partial }_r^4{\varphi }^{n+1}=\Delta t{\mathcal{S}}^{n+1}-{\partial }_r^2{\varphi }^n,\\ & -{\partial }_z^2{\varphi }^{n+1}+\eta {\partial }_z^2{q}^{n+1}=-{\partial }_z{\mathcal{F}}^{n+1},\end{array}$$

$$\{\begin{array}{rl}& {\partial }_r{\varphi }^{n+1}={\partial }_r^3{\varphi }^{n+1}=0\text{on}{\Sigma }_0\cup {\Sigma }_l\cup {\Sigma }_{L_r};\\ & \text{periodic boundary conditions for both ϕn+1 and qn+1 on Γ0∪ΓLz};\\ & {\partial }_z{\varphi }_a^{n+1}=\eta \left(1-e^{\Lambda -{\varphi }_a^{n+1}}\right)+{\mathcal{F}}_a^{n+1}\text{on}{\Gamma }_a;\\ & {\partial }_z{\varphi }_b^{n+1}=-\eta \left(1-e^{\Lambda -{\varphi }_b^{n+1}}\right)+{\mathcal{F}}_b^{n+1}\text{on}{\Gamma }_b;\\ & {\partial }_z{q}_a^{n+1}=\left(1-e^{\Lambda -{\varphi }_a^{n+1}}\right)\text{on}{\Gamma }_a;\\ & {\partial }_z{q}_b^{n+1}=-\left(1-e^{\Lambda -{\varphi }_b^{n+1}}\right)\text{on}{\Gamma }_b;\\ & q^{n+1}\equiv q_{\star }\text{on}{\Gamma }_q,\end{array}$$

$$e^{\Lambda -{\varphi }_{\alpha }^{n+1}}=e^{\Lambda -{\varphi }_{\alpha }^n}{e}^{{\varphi }_{\alpha }^n-{\varphi }_{\alpha }^{n+1}}\simeq e^{\Lambda -{\varphi }_{\alpha }^n}\left(1+{\varphi }_{\alpha }^n-{\varphi }_{\alpha }^{n+1}\right)\left(\alpha =a,b\right),$$

$$\{\begin{array}{rl}& {\partial }_z{\varphi }_a^{n+1}-\eta e^{\Lambda -{\varphi }_a^n}{\varphi }_a^{n+1}=\eta \left[1-e^{\Lambda -{\varphi }_a^n}\left(1+{\varphi }_a^n\right)\right]+{\mathcal{F}}_a^{n+1}\text{on}{\Gamma }_a,\\ & {\partial }_z{\varphi }_b^{n+1}+\eta e^{\Lambda -{\varphi }_b^n}{\varphi }_b^{n+1}=-\eta \left[1-e^{\Lambda -{\varphi }_b^n}\left(1+{\varphi }_b^n\right)\right]+{\mathcal{F}}_b^{n+1}\text{on}{\Gamma }_b,\end{array}$$

$$\{\begin{array}{rl}& {\partial }_z{\varphi }_a^{n+1}-\eta e^{\Lambda -{\varphi }_a^n}{\varphi }_a^{n+1}=\eta \left[1-e^{\Lambda -{\varphi }_a^n}\left(1+{\varphi }_a^n\right)\right]+{\mathcal{F}}_a^{n+1}\text{on}{\Gamma }_a,\\ & {\partial }_z{\varphi }_b^{n+1}+\eta e^{\Lambda -{\varphi }_b^n}{\varphi }_b^{n+1}=-\eta \left[1-e^{\Lambda -{\varphi }_b^n}\left(1+{\varphi }_b^n\right)\right]+{\mathcal{F}}_b^{n+1}\text{on}{\Gamma }_b,\\ & {\partial }_z{q}_a^{n+1}-e^{\Lambda -{\varphi }_a^n}{\varphi }_a^{n+1}=\left[1-e^{\Lambda -{\varphi }_a^n}\left(1+{\varphi }_a^n\right)\right]\text{on}{\Gamma }_a,\\ & {\partial }_z{q}_b^{n+1}+e^{\Lambda -{\varphi }_b^n}{\varphi }_b^{n+1}=-\left[1-e^{\Lambda -{\varphi }_b^n}\left(1+{\varphi }_b^n\right)\right]\text{on}{\Gamma }_b.\end{array}$$

$$\begin{array}{rl}{\left({\partial }_z^{2}u\right)}_{i,j}=& a_i^{-}{u}_{i-1,j}-\left(a_i^{-}+a_i^{+}\right)u_{i,j}+a_i^{+}{u}_{i+1,j},\\ {\left({\partial }_r^{2}u\right)}_{i,j}=& b_j^{-}{u}_{i,j-1}-\left(b_j^{-}+b_j^{+}\right)u_{i,j}+b_j^{+}{u}_{i,j+1},\\ {\left({\partial }_r^{4}u\right)}_{i,j}=& \left(b_{j-1}^{-}{b}_j^{-}\right)u_{i,j-2}-b_j^{-}\left(b_{j-1}^{-}+b_{j-1}^{+}+b_j^{-}+b_j^{+}\right)u_{i,j-1}+\\ & +\left(b_{j-1}^{+}{b}_j^{-}+{\left(b_j^{-}+b_j^{+}\right)}^2+b_j^{+}{b}_{j+1}^{-}\right)u_{i,j}+\\ & -b_j^{+}\left(b_j^{-}+b_j^{+}+b_{j+1}^{-}+b_{j+1}^{+}\right)u_{i,j+1}+\left(b_j^{+}{b}_{j+1}^{+}\right)u_{i,j+2},\end{array}$$

$$a_i^{-}=\frac{2}{\Delta z_i^{-}\left(\Delta z_i^{-}+\Delta z_i^{+}\right)},a_i^{+}=\frac{2}{\Delta z_i^{+}\left(\Delta z_i^{-}+\Delta z_i^{+}\right)}i=1,\dots M_z,$$

$$b_j^{-}=\frac{2}{\Delta r_j^{-}\left(\Delta r_j^{-}+\Delta r_j^{+}\right)},b_j^{+}=\frac{2}{\Delta r_j^{+}\left(\Delta r_j^{-}+\Delta r_j^{+}\right)}j=1,\dots M_r,$$

$$\Delta z_i^{-}=z_i-z_{i-1},\Delta z_i^{+}=z_{i+1}-z_i,\Delta r_j^{-}=r_j-r_{j-1},\Delta r_j^{+}=r_{j+1}-r_j.$$

$${\varphi }_{ex}^{(M)}(r,z):=\sin \left(2\pi z\right)\cos \left(\pi r\right)+\Lambda +\eta \sin \left(2\pi z\right),$$

$${\mathcal{S}}^{(M)}(r,z):=4{\pi }^2\sin \left(2\pi z\right)+\nu {\pi }^4\sin \left(2\pi z\right)\cos \left(\pi r\right),$$

$${\mathcal{F}}^{(M)}(r,z):=2\pi \cos \left(2\pi z\right)\cos \left(\pi r\right)+2\pi \eta .$$

$${\mathcal{S}}^{(P)}(r,z):=2\cdot 10^{-3}\exp \left(-20L_r{\left(r-\frac{3}{4}{L}_r\right)}^2\right),$$

Full text

1. Introduction

The subject matter of the present paper is related to the magnetically confined fusion plasmas with the objective to contribute by some means to the improvement of the numerical schemes used for the simulation of the plasma evolution in a tokamak.

Succeeding to produce energy via thermonuclear fusion processes in a tokamak, is strongly dependent on the aptitude to confine the plasma in the core of the tokamak and at the same time on the ability to control the plasma heat-flux on the wall of the device. These two requirements (sine qua non) are incommoded by the turbulent plasma transport occurring in such high temperature environments. So, one of the most important research fields at the moment is the comprehension, the prediction as well as the control of the turbulent plasma flow in a tokamak, in particular in the edge region of such a device. What we mean with “edge” or “peripheral” region of the tokamak is the region constituted firstly of the open magnetic field lines of the SOL (Scrape-off Layer), intercepting the wall (at the limiter), and secondly of the closed field line region, nearby the separatrix.

The study of peripheral tokamak plasmas is crucial for several reasons. First of all, this region imposes boundary conditions for the core plasma, and has thus to be treated with care. In particular, it is very important for the confinement properties of the reactor, to understand the turbulences occurring there. Secondly, it is the SOL-plasma which enters into contact with the wall of the reactor, such that controlling this peripheral plasma can be important for the protection of the reactor.

In this peripheral region, fluid models are often employed to describe the plasma evolution. This is due to the befalling smaller temperatures, which induce smaller particle mean-field paths such that closure relations, such as the Braginskii closures [2], are valid and permit to obtain macroscopic fluid models from more accurate kinetic ones. The TOKAM3X code [20, 21], which is at the basis of the present work, is founded on such a fluid approach. The strong magnetic configuration has as an effect that the underlying conservation laws are strongly anisotropic. The plasma dynamics parallel to the magnetic field lines is very fast, whereas the dynamics in the perpendicular direction (with respect to the magnetic field) is more constraint. Multiple scales are thus present in the governing conservation laws, which render the mathematical as well as numerical study very challenging. Specific numerical algorithms have hence to be designed in order to capture accurately the physics of interest, without unnecessary small time and space steps. This is of primary importance, in order to reduce computational burden in such huge simulations.

The aim of the present paper is hence to contribute to the advance of the TOKAM3X code by proposing an efficient numerical scheme for a singularly-perturbed equation, which seems to be for the moment one of the troublesome points of this code. Indeed, the vorticity equation, which permits to compute the electric potential $\phi$, is strongly anisotropic, leading numerically to an ill-conditioned problem. Standard numerical techniques require refined meshes in order to account for the strong anisotropy and to get accurate results. This leads necessarily to excessive computational costs and multi-scale numerical schemes could be a welcome alternative. The Asymptotic-Preserving scheme we propose in this paper is based on a microscopic-macroscopic decomposition of the unknown function (electric potential $\phi$) separating in this manner the fast and rapid scales in the problem and permitting to transform the original singularly-perturbed problem into a regularly-perturbed one. This last one is then simply discretized by a standard finite-difference method. This procedure allows the choice of meshes independent on the perturbation parameter $\eta$, which describes the anisotropy, the numerical results remaining uniformly accurate for all ranges of $\eta\in[0,1]$.

Several types of AP-schemes have been introduced in literature, for various types of singularly-perturbed problems. To mention only some examples, we refer the interested reader to [3, 4, 8, 10, 12, 11, 17] and references therein. Briefly, Asymptotic-Preserving schemes are efficient procedures for solving singularly-perturbed problems. They consist in trying to mimic, at the discrete level, the asymptotic behaviour of the problem-solution as the perturbation parameter $\eta$ tends towards zero. This property renders the AP-scheme uniformly accurate, with respect to $\eta$. The AP-scheme proposed in the present paper is based on the experience of the authors acquired with the study of highly anisotropic elliptic [5, 6, 7] or parabolic equations [14, 15].

The structure of this paper is the following: In Section 2 we shall present in detail the mathematical problem and the numerical difficulties one can encounter for its resolution. Section 3 presents and investigates mathematically an Asymptotic-Preserving reformulation for the underlying singularly-perturbed problem, reformulation which suites better for the limit of vanishing resistivity parameter, i.e. $\eta\rightarrow 0$. Section 4 deals then with the discretization of the AP-method we introduced above. And finally, in Section 5 we shall first validate the AP-scheme, and then apply it in the case of a real thermonuclear plasma experiment.

2. The mathematical model

Let us now present the singularly-perturbed, non-linear problem, describing the evolution of the electric potential $\phi$ in the peripheral tokamak region, represented by the domain $\Omega\subset\mathbb{R}^{2}$, which is sketched in Fig. 1. This equation reads

[TABLE]

and is completed with the initial condition

[TABLE]

for some given function $\phi_{in}$. The imposed boundary conditions are the no-slip boundary conditions on $\Sigma\coloneqq\Sigma_{0}\cup\Sigma_{l}\cup\Sigma_{L_{r}}=\{(r,z)\in\partial\Omega\,\,/\,\,r=0,\,r=l\,\,\textrm{or}\,\,r=L_{r}\}$

[TABLE]

periodic boundary conditions on $\Gamma_{0}\cup\Gamma_{L_{z}}=\{(r,z)\in\partial\Omega\,\,/\,\,z=0\,\,\textrm{or}\,\,z=L_{z}\}$, that means in the core of the plasma, and the nonlinear sheath boundary conditions on the limiters $\Gamma_{a}\cup\Gamma_{b}=\{(r,z)\in\partial\Omega\,\,/\,\,z=a\,\,\textrm{or}\,\,z=b\}$

[TABLE]

The source term is composed of a stiff and a non-stiff part, denoted respectively by $\partial_{z}{\mathcal{F}}$ and ${\mathcal{S}}$, whereas $\eta>0$, $\nu>0$ and $\Lambda>0$ are some given constants. The parameter $\Lambda$ is the sheath floating potential and $\eta$ represents the parallel resistivity of the plasma and is supposed to be very small $0<\eta\ll 1$, resulting in a singularly-perturbed problem, which is very challenging to solve numerically. The domain $\Omega$ covers both, closed flux surfaces (periodic region) and open flux surfaces, touching the limiter (nonlinear boundary conditions), describing thus the so-called SOL region (Scrape-off Layer) of a tokamak plasma.

Our model problem (2.1)–(2.4) is extracted from the TOKAM3X model [20, 21] which describes the dynamics of a magnetically confined tokamak edge-plasma via a fluid approach. At the basis of the TOKAM3X code are the balance equations for the electrons and ions, coupled to the Poisson equation for the electrostatic potential. Under some suitable assumptions, such as for example the quasi-neutrality condition, the low mass ratio $m_{e}/m_{i}$ assumption, the drift-approximation, and so on, one obtains a fluid model, based on the particle balance equation, the parallel momentum equation, the charge balance ($\nabla\cdot j=0$) and the parallel Ohm’s law, describing the evolution of respectively the electron density $N$, the parallel ion momentum $\Gamma$, the plasma potential $\phi$ and the parallel current $j_{||}$, parallel with respect to the imposed strong magnetic field. Introducing the vorticity quantity $W$, the set of equations at the foundation of the TOKAM3X code, read

[TABLE]

with the vorticity $W$ and the parallel current $j_{||}$ defined as

[TABLE]

In this system ${\mathbf{B}}$ is the magnetic field, with direction ${\mathbf{b}}\coloneqq{{\mathbf{B}}\over|{\mathbf{B}}|}$, $S_{N}$ is a source term, $P$ is the static pressure defined as $P\coloneqq N\,(T_{i}+T_{e})$ and $\eta_{||}$ is the normalized parallel collisional resistivity of the plasma. Furthermore ${\mathbf{u}_{i,e}}={u^{i,e}}_{||}\,{\mathbf{b}}+{\mathbf{u}^{i,e}}_{\perp}$ denote the particle macroscopic velocities, with the perpendicular parts ${\mathbf{u}^{i,e}}_{\perp}={\mathbf{u}}_{\mathbf{E}}+{\mathbf{u}^{i,e}_{\nabla B}}$ characterized in terms of the electric and curvature drift velocities, given by

[TABLE]

and where $T_{i,e}$ are the respective particle temperatures, considered in the present case as constant. The diffusion coefficients $D_{N,\Gamma,W}$ account for the collisional transport and/or diffusive transport (neoclassical and anomalous) and model the turbulences at small scales. System (2.5) allows to study the isothermal, electrostatic, 3D turbulences arising in the SOL. For its detailed derivation from the underlying conservation laws, we refer the interested reader to the references [20, 21]. The numerical simulations as well as a detailed analysis of the obtained results are also presented there. The authors remark that there are still some numerical difficulties to be further inquired in the vorticity equation, and this due to the strong anisotropy of the problem. In order to address these numerical difficulties, we decided to extract in this paper the vorticity equation from the whole system (2.5) and simplify it, keeping only the terms which cause complications. A simplified version of this potential equation is given by our problem (2.1)–(2.4) (see previous work [18] for its derivation from (2.5)).

The rigorous mathematical study, meaning the study of the well-posedness of this problem $\forall\eta>0$ (existence, uniqueness and stability of a weak solution), has been considered in [18]. We observed in that work that the primary mathematical difficulties (for fixed $\eta>0$) arise on one hand from the degeneracy in time of the problem (regularization) and on the other hand from the non-linear boundary conditions. These difficulties will naturally also occur in the numerical treatment of this problem.

The main goal of the present paper is to propose a numerical follow-up of the latter more mathematical paper, in the aim to design an efficient numerical scheme for the resolution of problem (2.1)–(2.4). Typical computational challenges within this task come from:

• •

the high anisotropy of the problem, described by the small resistivity parameter $\eta\ll 1$, and introduced into the problem by the strong magnetic field, which confines the plasma;

• •

the non-linearity of the boundary conditions, describing the plasma-wall interactions (Bohm-conditions).

Concerning the first point, the anisotropy, we have to deal with a typical singularly-perturbed problem when $\eta\rightarrow 0$. Such kind of problems are particularly difficult to treat numerically, as the equations change type as the perturbation parameter $\eta$ tends towards zero, leading (usually) to ill-posed problems in the limit. In the present case, the so-called “reduced-model” has the form

[TABLE]

associated with the other boundary conditions on $\Sigma$ and $\Gamma_{0}\cup\Gamma_{L_{z}}$ as well as the initial conditions. One can now remark that this reduced model is ill-posed, as it admits either no solution (if the initial condition is not well-prepared) or an infinite amount of solutions, as one can add to any solution another $r$-dependent function, satisfying the boundary conditions on $\Sigma$ and the initial condition. On the discrete level, the $(R)$-problem will always have infinitely many solutions, as one steps over the initial condition. This distinction between well-prepared and not well-prepared initial condition is related to the creation of a boundary layer near $t=0$.

At the discrete level, all these complications will be translated into the fact that the linear system to be solved will become singular, as $\eta\rightarrow 0$, or equally, becomes ill-conditioned as $\eta\ll 1$, leading hence to erroneous results. For not too small, fixed $\eta$-values, a preconditioner could help, however in a general case, where the perturbation parameter $\eta$ varies within the simulation domain $\Omega$ and takes various orders of magnitude, the preconditioner will no longer follow and new techniques have to be employed to rescue the user.

The occurrence of all these difficulties (theoretical as well as numerical) is strongly related to the multi-scale character of our problem. Indeed, for small $\eta\ll 1$ the problem is evolving more rapidly in the $z$-direction, which represents the direction of the strong magnetic field, than in the perpendicular $r$-direction. Different space- and time-scales are hence introduced in the problem by the perturbation parameter $\eta$, and in order to be accurate, a standard numerical scheme has to take into account for all these small scales, by imposing very restrictive grid-conditions as for example meshes of order $\eta$. This can become rapidly too costly.

The aim of the next sections will be hence to introduce a multi-scale numerical scheme, based on the study of the asymptotic behaviour of the solution $\phi^{\eta}$ as $\eta$ goes to zero. A decomposition of $\phi^{\eta}$ into a macroscopic part and a microscopic one separates somehow the different dynamics in the problem. This procedure permits then the use of judicious mesh-sizes and time-steps, adapted to the physical phenomenon one wants to study and not to the perturbation parameter and permits to treat with no huge computational costs even the limiting case $\eta\equiv 0$.

3. An asymptotic-preserving reformulation for the evolution problem

This section is devoted to a reformulation of the original singularly-perturbed electric potential problem (2.1)–(2.4), denoted in the sequel by $(SP)_{\eta}$, into a problem which behaves better in the limit $\eta\rightarrow 0$.

The essence of our numerical method for the resolution of (2.1)–(2.4) is based on the following micro-macro decomposition of the unknown $\phi^{\eta}$

[TABLE]

where the macroscopic function $p^{\eta}$ is chosen to be solution of the “dominant” problem

[TABLE]

associated with the other homogeneous or periodic conditions on the remaining boundaries. For the uniqueness of the decomposition, we shall further fix $q^{\eta}$ (or equivalently $p^{\eta}$) on an interface, denoted here $\Gamma_{q}$, as follows

[TABLE]

Indeed, the decomposition (3.7) is now unique for given $\phi^{\eta}$ and fixed $\eta>0$, precisely due to the fact that we impose $q^{\eta}$ on $\Gamma_{q}$, fixing thus also $p^{\eta}$ on this interface. Problem (3.8) becomes thus a well-posed elliptic problem with unique solution $p^{\eta}$ for all $\eta>0$.

With this decomposition, the problem $(SP)_{\eta}$ transforms now into the completely equivalent system

[TABLE]

associated with the usual initial condition for $\phi^{\eta}$, the usual homogeneous resp. periodic boundary conditions on $\Sigma$ resp. $\Gamma_{0}\cup\Gamma_{L_{z}}$ and the following boundary conditions on $\Gamma_{a}\cup\Gamma_{b}$ for the unknowns $(\phi^{\eta},q^{\eta})$

[TABLE]

Note that no boundary condition for $q^{\eta}$ is needed on $\Sigma$, as no $r$-derivatives of $q^{\eta}$ occur in the system. We shall call in the following this problem the Asymptotic-Preserving reformulation of our Singularly-Perturbed problem (2.1)–(2.4), denoted simply by $(AP)_{\eta}$.

The equivalence of both problems, $(SP)_{\eta}$ and $(AP)_{\eta}$, for any $\eta>0$, is due to the uniqueness of the solution of (3.8) when imposing $p^{\eta}_{|\Gamma_{q}}=\phi^{\eta}_{|\Gamma_{q}}-\eta\,q_{\star}$. This equivalence together with the mathematical existence and uniqueness studies of the problem $(SP)_{\eta}$, considered in [18], permits to show that $(AP)_{\eta}$ is well-posed for each $\eta>0$. The essential difference between these two reformulations is perceived only in the limit $\eta\rightarrow 0$. Indeed, $(SP)_{\eta}$ becomes singular, as explained in Section 2. Let us yet formally investigate what happens with the micro-macro reformulation $(AP)_{\eta}$, when $\eta$ goes to zero. Setting formally $\eta\equiv 0$ one gets the system

[TABLE]

associated with the following boundary conditions for the unknowns $(\phi^{0},q^{0})$

[TABLE]

which is a type of saddle-point problem. The unknown $q^{0}$ can be seen here as a Lagrangian multiplier, associated to the constraint $-\partial^{2}_{z}\phi^{0}=-\partial_{z}{\mathcal{F}}$. The rigorous well-posedness of this limit problem is not the aim of the present paper, which is much more numerical, can however be a nice extension of this work, involving saddle-point theory.

Hence, the main benefit of the AP-reformulation $(AP)_{\eta}$ is that in the limit $\eta\rightarrow 0$ one gets a well-posed problem, such that we have no more to face numerical singularities, when solving the reformulation for small $\eta$-values. This big advantage shall be extensively put into light with the simulations performed in Section 5.

Before proceeding to the numerical treatment, let us mention here some words about the essence of the micro-macro decomposition (3.7), which is at the basis of our AP-reformulation. The idea behind was to eliminate the dominant, stiff operator, and this has been done by introducing a sort of separation of scales. The function $p^{\eta}$ solves the dominant operator, whereas $q^{\eta}$ incorporates the microscopic information. In the limit $\eta\rightarrow 0$ the microscopic part $q^{0}$ is still present in the homogenized limit model $(AP)_{0}$, and it is this part which renders the problem well-posed. It permits to recover the microscopic information, which was lost in the reduced model (2.6).

4. The numerical discretization via finite differences

In contrast to previous papers on highly anisotropic elliptic or parabolic problems [5, 6, 7, 14, 15], we made here the choice to solve the model equations outlined in Eqs. (2.1)–(2.4) and (3.10)–(3.11) by means of a numerical approach based on finite difference approximations, instead of relying on the finite element method. The reason for this choice is the fact that we would like to provide the team developing the TOKAM3X code, based on a discretization of the balance equations via the finite volume method, with a technique directly applicable in that context in a way as straightforward as possible.

For the investigation of the consistency and accuracy of the numerical scheme proposed here, a Cartesian uniform mesh with constant mesh size along both longitudinal and radial directions ($\Delta r\equiv\Delta z=const$) has been adopted (see Section 5.2). This constant uniform mesh may however be abandoned to introduce mesh refinement near the borders of the domain, or where steep gradients of the solution are to be expected. Such a non-uniform Cartesian mesh has been adopted in the second part of our investigation, when the analysis of a problem setting inspired by a real physical plasma application is proposed. This will be presented in Section 5.3.

The semi-discretization in space is not the most critical part in the construction of an AP scheme, though. In fact, special care must be paid to the time-discretization, in particular when decisions have to be taken concerning which terms to take implicitly and which ones explicitly. In order not to destroy the desirable properties of our AP formulation, the time discretization is here based on the implicit Euler scheme. The interval $[0,\,T]$ is discretized in uniform steps $\Delta t$, the solution being evaluated at the time instants $t_{0},t_{1},\ldots t_{N_{t}}$ defined as

[TABLE]

Let us present now the discretizations of both formulations, the $\left(AP\right)_{\eta}$-scheme formulated in Eqs. (3.10)–(3.11) respectively the $\left(SP\right)_{\eta}$-scheme formulated in Eqs. (2.1)–(2.4). In the following Section 5, we shall then compare the performances of these two formulations.

4.1. Semi-discretization in time of $\left(SP\right)_{\eta}$.

Discretizing Eq. (2.1) in time by means of the implicit Euler scheme, yields

[TABLE]

and hence for $n=0,\ldots,N_{t}-1$,

[TABLE]

where $\phi^{k}$, $\mathcal{S}^{k}$ and $\mathcal{F}^{k}$ denote the solution and the source terms $\mathcal{S}$ resp. $\mathcal{F}$, evaluated at the $k$-th time step, and $\phi^{0}$ is given by the initial condition $\phi_{in}$ (2.2).

This equation is associated with the following nonlinear boundary conditions:

[TABLE]

where we used the notation $\phi^{n+1}_{a}\equiv\phi^{n+1}\left(t,r,z=a\right)$, $\phi^{n+1}_{b}\equiv\phi^{n+1}\left(t,r,z=b\right)$ and similarly for ${\mathcal{F}}^{n+1}$.

4.2. Semi-discretization in time of $\left(AP\right)_{\eta}$.

Discretizing Eq. (3.10) in time by means of the implicit Euler scheme yields for $n=0,1,\ldots,N_{t}-1$

[TABLE]

where $q^{n+1}$ denotes the microscopic function $q$ at the $(n+1)$-th time step. Eq. (4.16) is associated with the following nonlinear boundary conditions:

[TABLE]

where $q_{\star}$ is an arbitrary constant.

4.3. Linearization of the boundary conditions.

The treatment of the nonlinear boundary conditions is the same for both formulations $\left(SP\right)_{\eta}$ or $\left(AP\right)_{\eta}$, and is based on a linearization obtained via a Taylor expansion. The following approximation has been adopted:

[TABLE]

leading to the following Robin boundary condition, substituting the boundary conditions on $\Gamma_{a}$ and $\Gamma_{b}$ for the $\left(SP\right)_{\eta}$ scheme (see Section 4.1),

[TABLE]

and the following conditions substituting the boundary conditions on $\Gamma_{a}$ and $\Gamma_{b}$ for the $\left(AP\right)_{\eta}$ scheme (see Section 4.2):

[TABLE]

4.4. Semi-discretization in space.

The computational domain $\Omega$, sketched in Fig. 1, has been discretized by means of a structured Cartesian grid. We shall denote by $M_{z}$ the maximum number of grid nodes in the $z$ (longitudinal) direction (the node located on $z=L_{z}$ is not included in the computational domain – and hence also in $M_{z}$ – because of the periodic boundary conditions on $\Gamma_{0}$ and $\Gamma_{L_{z}}$), and by $M_{r}$ the maximum number of grid nodes in the $r$ (radial) direction. Since the domain is not a square, the total number of grid points, denoted by $M$, is in general $M<M_{z}M_{r}$.

In the case of the $\left(SP\right)_{\eta}$ scheme, Eq. (4.14) is discretized in the $M$ grid points of the computational domain, leading to a system of algebraic equations $\mathbf{A}\mathbf{x}=\mathbf{b}$ with a number $M_{U}$ of unknowns matching the number of grid nodes, i.e. $M_{U}=M$, where the components of the vector $\mathbf{x}$ represent the grid values of the function $\phi$.

In the case of the $\left(AP\right)_{\eta}$ scheme, Eq. (4.16)1 is discretized on the $M$ grid points of the computational domain, and Eq. (4.16)2 is discretized on the $M-M_{r}$ grid nodes of the computational domain except the interface $\Gamma_{q}$ (due to the fact that the value of the function $q$ is prescribed on $\Gamma_{q}$, as explained in Section 3), leading to $M_{U}=2M-M_{r}$ unknowns. The components of the vector $\mathbf{x}$ represent in this case the $M$ grid values of the function $\phi$ and the $M-M_{r}$ grid values of the microscopic function $q$.

Standard finite difference approximations have been used to discretize the second and fourth order derivatives appearing in Eq. (4.14) and Eq. (4.16). Denoting with $u_{i,j}$ the value of the unknown variable $u$ ($u=\phi,q$) in the node located in the $i$-th column and $j$-th row of the Cartesian grid, whose $z$ and $r$ coordinates are, respectively, $z_{i}$ and $r_{j}$, we have

[TABLE]

with

[TABLE]

and

[TABLE]

The boundary conditions given in Eq. (4.19) (for the $\left(SP\right)_{\eta}$ scheme) or Eq. (4.20) (for the $\left(AP\right)_{\eta}$ scheme) are then properly used to eliminate the ghost-nodes.

This discretization of the differential equations, together with the shape of the computational domain, result in a sparse matrix $\mathbf{A}$ of the emerging linear system having a peculiar pattern. This pattern is represented in Fig. 2 for both $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ formulations of the problem, for a particularly poor discretization of the computational domain for the sake of the visualization ($M_{z}=12$, $M_{r}=8$, $M=80$; $M_{U}=80$ for the $\left(SP\right)_{\eta}$, $M_{U}=151$ for the $\left(AP\right)_{\eta}$ problem).

In the case of the $\left(SP\right)_{\eta}$ formulation, it is clearly visible that the band of the sparse matrix changes in correspondence to the change of the longitudinal size of the domain $\Omega$. This pattern is reproduced in the upper left block of the matrix pertaining to the $\left(AP\right)_{\eta}$ formulation (in this case, the unknowns numbered from $1$ to $M$ represent the grid values of the variable $\phi$, and the unknowns ranging from $M+1$ to $M_{U}$ represent the grid values of the microscopic function $q$).

The linear system associated to this sparse matrix is solved by means of the MUMPS library [16]. The error analysis, provided in Section 5, is carried out by means of an estimate of the upper bound of the error affecting the solution [1]– a metric which turns out to be particularly meaningful for sparse linear systems arising from a physical application – as well as by means of the traditional condition number [19]. The first metric is directly provided by the MUMPS library; the second is calculated with the aid of the linalg module of the SciPy software package [13].

5. Simulation results

In this section, selected numerical solutions are presented for two particular case studies. For this purpose, we have designed and implemented a code that allows to solve the model equations Eqs. (2.1)–(2.4) via both the $\left(SP\right)_{\eta}$ scheme – leading to Eqs. (4.14)–(4.15) – and the $\left(AP\right)_{\eta}$ scheme, leading to Eqs. (4.16)–(4.17).

The first case, a mathematical example denoted as Case (M), is a set-up for which a time-independent analytical solution has been constructed. The aim of this investigation is twofold. The primary aim is to validate the developed codes by comparing the numerical solution to the exact analytical solution. The second main purpose is to compare the numerical solution obtained by means of the $\left(SP\right)_{\eta}$ formulation to the one obtained with the $\left(AP\right)_{\eta}$ formulation, in order to highlight how the $\left(AP\right)_{\eta}$ approach allows to overcome the difficulties encountered when using the $\left(SP\right)_{\eta}$ approach, as the perturbation parameter $\eta\to 0$.

The second case, a physical example denoted as Case (P), is a set-up closer to a real physical scenario inspired by the Tokamak configuration of interest for the research group developing the TOKAM3X code. The main purpose here is to investigate if the $\left(AP\right)_{\eta}$-based formulation represents a viable approach to the numerical study of the problem in a practical situation.

The details about the domain size, the model parameters and other settings of the numerical algorithms for the two above-mentioned cases are detailed in Section 5.1. In Section 5.2 and in Section 5.3, selected numerical solutions for Case (M) and Case (P), respectively, are presented.

5.1. Presentation of the test cases

The mathematical example, Case (M), regards a setting of the problem for which a time-independent exact solution was explicitly constructed. The set of model parameters ($\nu$, $\Lambda$ and $\eta$), as well as the geometric configuration (dimensions $a$, $b$, $L_{z}$, $l$, $L_{r}$; see Fig. 1), are listed in Tab. 1.

The exact analytical solution $\phi$, denoted as $\phi^{(M)}_{ex}$, along with the source terms $\mathcal{S}$ and $\mathcal{F}$, denoted respectively as ${\mathcal{S}}^{(M)}$ and $\mathcal{F}^{(M)}$, are defined as follows:

[TABLE]

It is worth noticing that the parameters listed in Tab. 1 as well as the source terms $\mathcal{S}^{(M)}$, $\mathcal{F}^{(M)}$ and the solution $\phi^{(M)}_{ex}$, do not have any relevant physical meaning. The only purpose of their choice is to have at our disposal a relatively simple analytical solution to Eqs. (2.1)–(2.4) that allows us to validate the numerical code and to investigate the good properties of the asymptotic-preserving formulation of the problem.

The physical example, Case (P), regards a more physically-oriented set-up of the problem. The model parameters and the geometric configuration – also listed in Tab. 1 – completely characterize the physical systems together with the source terms ${\mathcal{S}}$ and $\mathcal{F}$, denoted in this case respectively as $\mathcal{S}^{(P)}$ and $\mathcal{F}^{(P)}$, and defined as

[TABLE]

The analysis of the configuration Case (P) is motivated by the purpose of applying the proposed numerical scheme to a real physical application as the study of a Tokamak plasma: in this respect, Case (P) resembles (despite significant simplifications and idealizations) to the actual configuration and working scenario of a Tokamak plasma of interest for the studies carried out by the team developing the numerical code TOKAM3X [9].

In both settings – Case (M) and Case (P) – the initial data $\phi_{in}\left(r,z\right)\coloneqq\phi\left(t=0,r,z\right)$ is defined as follows:

[TABLE]

5.2. Numerical investigations of Case (M)

Let us test now the performances of both, $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ formulations, in Case (M). An exact steady-state solution to Eqs. (2.1)–(2.4) is given in (5.24). In this case, the adopted grid is uniform, with constant and equal step size in each direction, $\Delta z=\Delta r=h$; a sketch of the domain with an example of its discretization is provided in Fig. 3.

5.2.1. Validation of the numerical code and order of convergence

The $L^{2}$-norm of the error of the numerical solution (with respect to the exact solution), as well as the order of convergence of the $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ numerical schemes, are shown in Tab. 2.

It is evident that the $\left(AP\right)_{\eta}$ numerical scheme guarantees good performances for any value of the parameter $\eta$ in the interval $[0,\,1]$; even values of $\eta$ as small as $10^{-14}$ or the limit case $\eta\equiv 0$ are not critical. Contrary to this, the $\left(SP\right)_{\eta}$ scheme can be used only for larger values of $\eta$. As Tab. 2 shows, the case $\eta=10^{-14}$ is only treatable with the $\left(AP\right)_{\eta}$ scheme, and even in the case $\eta=10^{-8}$, the $\left(SP\right)_{\eta}$ scheme does not allow to obtain reasonable numerical results (or any result at all) when the grid is fine.

The behavior of the $H^{1}$-norm of the error of the numerical solution (with respect to the exact solution) is shown in Fig. 4. Comparing these data to those shown in Tab. 2, it is easily seen that when considering the $L^{2}$-norm, the convergence is of second order in space, while it is of first order in the $H^{1}$-norm, which is completely standard.

5.2.2. Comparison of the $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ formulations

In order to compare the behavior of the $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ schemes for small values of $\eta$, it is useful to take a closer look at the linear system of algebraic equations emerging from the discretization of the model equations.

A classical metric which provides insight in the reliability and accuracy of the solution of a linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$, is the condition number of the matrix $\mathbf{A}$, defined as $\kappa\left(\mathbf{A}\right)\coloneqq\left\lVert\mathbf{A}^{-1}\right\rVert\left\lVert\mathbf{A}\right\rVert$. The larger the condition number, the closer to singular is the matrix. In this latter case, obtaining the solution to the linear system could not be feasible or, when feasible, the solution is typically affected by a large error and therefore should not be trusted as a reliable solution.

In more details, the relative effect of round-off errors on the solution of a linear system is bounded by the condition number of the matrix $\kappa_{2}\left(\mathbf{A}\right)=\left\lVert\mathbf{A}^{-1}\right\rVert_{2}\left\lVert\mathbf{A}\right\rVert_{2}$ in the following manner

[TABLE]

where $\mathbf{\tilde{x}}=\mathbf{x}+\delta\mathbf{x}$ is the computed numerical solution of $\mathbf{A}\mathbf{x}=\mathbf{b}$ and can be seen as the exact solution of the slightly perturbed linear system $(\mathbf{A}+\delta\mathbf{A})\,\mathbf{\tilde{x}}=\mathbf{b}+\delta\mathbf{b}$, perturbations being due to round-off errors.

This traditional relative error-estimate is very pessimistic, for example it does not take into account for the particular form of the matrix $\mathbf{A}$, or for the special form of the right-hand side $\mathbf{b}$. A more satisfying approach has been proposed by Arioli et al. in [1], where the sparsity of the matrix is kept in mind. The relative error of the solution to the linear system is estimated in that work as follows

[TABLE]

where the quantities $\omega_{1}$ and $\omega_{2}$ (backward errors) are defined as follows

[TABLE]

with $\mathbf{A}_{i}$ the $i$-th row of $\mathbf{A}$, $\left|\mathbf{A}\right|$ being the matrix with elements $\left|A_{ij}\right|$, $\left|\mathbf{b}\right|$ the vector with elements $\left|b_{i}\right|$, and $I_{2}$ representing the set of indices of the equations of the linear system such that $\left|\mathbf{A}\mathbf{\tilde{x}}-\mathbf{b}\right|_{i}$ is nonzero and $\left(\left|\mathbf{A}\right|\left|\mathbf{\tilde{x}}\right|+\left|\mathbf{b}\right|\right)_{i}$ is small, whereas $I_{1}$ representing the set of indices of the remaining equations (if $I_{2}$ is empty, then $\omega_{2}=0$).

Moreover, two condition numbers $\kappa_{\omega_{1}}$ resp. $\kappa_{\omega_{2}}$ of the system (not just of the matrix) are defined as follows:

[TABLE]

where $\mathbf{1}$ is the column vector with all unitary elements. In the evaluation of $\kappa_{\omega_{\iota}}$ only those equations with index $i\in I_{\iota}$ are considered ($\iota=1,2$).

In order to evaluate numerically the accuracy of the resolution of the linear system corresponding to our two problems $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$, in particular to examine the influence of the parameter $\eta$ on the obtained results, we plotted in Fig. 5 two graphs. Firstly, the standard condition number $\kappa_{2}\left(\mathbf{A}\right)=\left\lVert\mathbf{A}^{-1}\right\rVert_{2}\left\lVert\mathbf{A}\right\rVert_{2}$ of the matrix $\mathbf{A}$ associated to the $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ schemes, is plotted as a function of the perturbation parameter $\eta$. As expected, the results show that as the perturbation parameter $\eta\to 0$, the condition number $\kappa_{2}$ associated to the $\left(SP\right)_{\eta}$ increases as $\eta^{-1}$, signifying that the matrix $\mathbf{A}$ gets closer and closer to a singular-matrix. In contrast, the condition number $\kappa_{2}$ associated to the $\left(AP\right)_{\eta}$ scheme never reaches critical values and remains independent of the parameter $\eta$, clearly showing that the singularity of $\eta\to 0$ of the singularly perturbed problem (2.1) has been removed in the $\left(AP\right)_{\eta}$ formulation we propose.

Secondly, to be sure that the traditional condition number is not too pessimistic in the here treated case, we decided to consider also the new metric proposed in [1], allowing to compute a good estimate of the upper bound of the relative error of the computed solution $\mathbf{\tilde{x}}$ with respect to the exact solution $\mathbf{x}$ of $\mathbf{A}\mathbf{x}=\mathbf{b}$.

The analysis of the upper bound $\bar{\varepsilon}$ of the relative error of the computed solution, shown in Fig. 5, once more clearly points out how the solution corresponding to the $\left(SP\right)_{\eta}$-scheme rapidly becomes unreliable as $\eta\to 0$. On the other hand, the relative error affecting the solution corresponding to the $\left(AP\right)_{\eta}$-scheme is roughly constant and of the order $10^{-10}$–$10^{-8}$ also for $\eta\ll 10^{-6}$.

5.3. Numerical investigations of Case (P)

In this section, we discuss the numerical results obtained for a case inspired by a real physical application, namely the evolution of the electric potential in the peripheral plasma region of a tokamak reactor. The list of model parameters and the geometric configuration of this case are given in Tab. 1, while a sketch of the discretized computational domain is provided in Fig. 6. In contrast to the mathematical case previously discussed in Section 5.2, the computational grid is not uniform in this case, being more refined (i.e. with smaller step sizes $\Delta z$ and $\Delta r$) in the proximity of the borders $\Gamma_{l}$, $\Gamma_{a}$, $\Gamma_{b}$. A graphical representation of the source terms $\mathcal{S}^{(P)}$ and $\mathcal{F}^{(P)}$, given in Eqs. (5.27)–(5.28), is provided in Fig. 7.

5.3.1. Error estimates of the numerical solution

To start, we computed also in this physical case the condition number $\kappa_{2}$ and the upper bound of the relative error $\bar{\varepsilon}$, introduced in Section 5.2.2, for the $\left(AP\right)_{\eta}$ scheme and compared it to those obtained for the $\left(SP\right)_{\eta}$ scheme.

These metrics, plotted in Fig. 8, clearly show that – as expected – also in this physical case the $\left(SP\right)_{\eta}$ scheme does not provide reliable numerical results for values of the parameter $\eta$ below $10^{-8}$, while the $\left(AP\right)_{\eta}$-scheme always gives accurate numerical results independently on the parameter $\eta$. The results are naturally obtained with a fixed grid and $\eta$ varying in $[0,1]$, clearly underlying the Asymptotic-Preserving property of the numerical scheme we introduced in this work.

5.3.2. Selection of numerical results

Let us now study in more details the numerical results obtained for the physical Case (P).

In Fig. 9 and Fig. 10, we plotted the profiles of the solution $\phi^{\eta}$ as a function of $z$ at three different values of the radial coordinate, resp. as a function of $r$ at three different values of the longitudinal coordinate. The profiles obtained by means of the $\left(SP\right)_{\eta}$ scheme are compared to those obtained with the $\left(AP\right)_{\eta}$ scheme, for several values of the parameter $\eta$. It is easily seen that as $\eta$ decreases, the profiles obtained with the $\left(AP\right)_{\eta}$ scheme converge to the limit profile obtained with $\eta=0$. This correct, expected behavior is not observed in the profiles calculated by means of the $\left(SP\right)_{\eta}$ scheme. In this latter case, in fact, when $\eta$ is small enough ($\eta=10^{-9}$, in the case under investigation), the results are clearly not showing the expected trend. This behavior is naturally also pointed out by the $\eta$-evolution of the upper bound estimate of the relative error affecting the computed numerical solution of the discretized algebraic system associated to $\left(SP\right)_{\eta}$ (see Fig. 8). The results obtained by means of the $\left(SP\right)_{\eta}$ scheme are therefore not reliable as the parameter $\eta$ approaches zero.

Looking in more details at Fig. 9 and 10 and reminding the reduced problem (2.6), one can do a very nice observation. As stated in Section 2, the solution to (2.6) is not unique, as one can add an arbitrary $r$-dependent function $\psi(r)$ to one solution $\phi(r,z)$, in order to get another solution. The $z$-evolution of all these solutions is however well-defined. This specific feature can be now observed in Fig. 9 and 10. Indeed, for small $\eta$-values, instead of solving $\left(SP\right)_{\eta}$ the computer will solve the reduced problem (2.6), and will hence arbitrarily fix a function $\psi(r)$. Fig. 9 shows precisely this behaviour, as the $z$-evolution is identical for both, $\left(SP\right)_{\eta}$ and $\left(AP\right)_{\eta}$ schemes, however $\left(SP\right)_{\eta}$ has some difficulties below $\eta=10^{-8}$ to find the right constant, for fixed $r$. In Fig.  10, one notices immediately the wrong $r$-evolution for the $\left(SP\right)_{\eta}$-scheme as $\eta$ gets smaller and smaller.

The computations shown in Fig. 9 and Fig. 10 have been carried out with $M_{z}=1\,280$, $M_{r}=129$, $M=156\,992$ ($M_{U}=156\,992$ for the $\left(SP\right)_{\eta}$ scheme; $M_{U}=313\,855$ for the $\left(AP\right)_{\eta}$ scheme).

In Fig. 11, the evolution of the solution $\phi^{\eta}$ at four different time instants ($t=10$, $100$, $200$ and $t=\infty$) are represented. The initial data, as previously mentioned, is $\phi_{in}(z,r)=\Lambda$. The field $\phi^{\eta}$, obtained for $\eta=10^{-12}$ using the $\left(AP\right)_{\eta}$ scheme, evolves towards the asymptotic solution represented in Fig. 11(d) for large times. The last of the shown fields, corresponding to $t=\infty$, has been obtained with a time-independent version of the numerical code, and as such it represents the steady state solution.

In Fig. 12 a selection of the longitudinal and radial profiles of the solution $\phi^{\eta}$ are plotted for several time instants, including those presented in Fig. 11. From these profiles it is easy to observe that the solution converges to the steady-state solution labeled with $t=\infty$. Briefly, the AP-scheme seems to recover also very well the $t\rightarrow\infty$ asymptotics, and this for all values of $\eta\in[0,1]$. The numerical solutions graphically shown in Fig. 11 and Fig. 12 have been obtained with the $\left(AP\right)_{\eta}$ scheme, with $M_{z}=192$, $M_{r}=49$, $M=8\,280$ and $\Delta t=1$. Grid refinement did not lead to appreciable changes in the results.

All the numerical results presented in this work have been obtained working in double precision floating-point arithmetics on a single node workstation, based on an Intel(R) Core(TM) i7-4770 CPU at 3.40GHz with hyper-threading enabled. A OpenMP-based parallelization has been adopted (when convenient) for the calculation of the linear system coefficients, and the MUMPS library has been compiled as to make use of OpenMP directives. In all cases, the computational time (“wall time”) of each single computation ranges from few seconds to few minutes, such that no further form of software or hardware acceleration was needed for the purpose of the present study.

6. Conclusion

We introduced in this work an Asymptotic-Preserving numerical scheme for the resolution of the highly anisotropic vorticity equation arising in plasma modelling. Numerical simulations permitted to underline the advantages of this new scheme as compared to standard discretizations, used in present simulation codes. In particular our AP-scheme permits to choose the grid independent on the perturbation parameter $\eta$. This property can lead to an essential gain in memory and computational time, without loss of accuracy, for $\eta$-values below a certain threshold value. This threshold value can be situated in our test cases at approximately $\eta_{*}\in[10^{-8},10^{-6}]$. If the physical parameter $\eta$ is smaller than $\eta_{*}$, then our AP-methodology could be of important benefit for tokamak simulations. For larger $\eta$-values however, standard discretizations are to be preferred, as only one equation (for the unknown $\phi$) has to be solved, whereas our AP-scheme is constituted of two equations for $(\phi,q)$. Keep nevertheless in mind that the $\eta$-value can be varying in the simulation domain, taking different orders of magnitude in various parts of the domain. In such a case, the AP-formulation takes the advantage. Anyhow, our here presented AP-formulation has still to be extended to this situation of variable $\eta(x)$, as for the moment we treated only constant $\eta$-values. Another future step would be also to implement this new AP-method in the TOKAM3X code and validate its practical use. All this will be the aim of a future work.

Acknowledgments. The authors would like first thank Patrick Tamain, Guido Ciraolo and Davide Galassi for fruitful discussions on this paper. Furthermore, we would like to acknowledge support from the Italian National Group for Mathematical Physics (GNFM/INdAM). This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.