Recent Upgrades in a 2D Turbulent Transport Solver Based on a Hybrid Discontinuous Galerkin Method for the Simulation of Fusion Plasma in Tokamak

Abstract

The simulation of fusion plasmas in realistic magnetic configurations and tokamak geometries still requires the development of advanced numerical algorithms owing to the complexity of the problem. In this context, we propose a Hybrid Discontinuous Galerkin (HDG) method to solve 2D transport fluid equations in realistic magnetic and tokamak wall geometries. This high-order solver can handle magnetic equilibrium free structured and unstructured meshes allowing a much more accurate discretization of the plasma facing components than current solvers based on magnetic field aligned methods associated with finite-differences (volumes) discretization. In addition, the method allows for handling realistic magnetic equilibrium, eventually non steady, a critical point in the modeling of full discharges including ramp up and ramp down phases. In this paper, we introduce the HDG algorithm with a special focus on recent developments related to the treatment of the cross-field diffusive terms, and to an adaptive mesh refinement technique improving the numerical efficiency and robustness of the scheme. The updated solver is verified with a manufactured solution method, and numerical tests are provided to illustrate the new capabilities of the code.

1. Introduction

Research in magnetic confinement fusion plasmas explores the possibility of producing carbon-free electric power by using fusion in deuterium–tritium plasmas heated to temperatures up to ${10}^7$–${10}^8$ K, and confined by a magnetic field in machines of toroidal shape known as tokamaks. With ITER and the promise of burning plasmas, the control of heat exhaust in high energy confinement configurations has become a topic of critical importance for the operation [1]. The difficulty to get global experimental measurements in tokamak makes complementary numerical simulations in realistic tokamak conditions a valuable asset to design optimized plasma scenarios, allowing for controlling the heat outfluxes and to prevent material damages. However, such numerical simulation remains a very challenging issue. This problem is multi-physics and multi-scales due to plasma wall interactions and turbulence. The geometry also adds a complexity in realistic configurations due to the shape of the tokamak wall and of the magnetic equilibrium. In addition, the strong anisotropy of the magnetic field components leads to a preferred orientation denoted as the parallel direction, with reference to the direction along the magnetic field lines. This leads to specific numerical issues as ill-conditioned algebraic operators to invert, and significant spurious numerical diffusion in the direction orthogonal to the anisotropy direction. Routine simulations able to provide information in acceptable timings in a tokamak of the size of ITER are still today restricted to 2D models based on averaged axisymmetric fluid–drift Braginskii equations [2,3,4]. The current code of the fusion community is generally based on first and second order finite differences or finite volumes, and so on structured meshes. Their discretization is aligned along the magnetic field lines to take advantage of the transport features and limit spurious numerical diffusion [5]. Thus, the accurate discretization of realistic tokamak wall geometries as well as plasma regions around singularities such the X-point or the tokamak center remains challenging with these codes. In addition, the simulation of transient phases of the plasma discharge when the plasma equilibrium moves is not affordable without a very expensive on the fly re-meshing of the computational domain. To overcome these limitations, we have recently considered a Hybrid Discontinuous Galerkin (HDG) method. Such discretization based on structured or unstructured meshes is magnetic equilibrium free that allows for accurate simulations of the whole vacuum chamber whatever the geometrical complexity of the tokamak wall or the magnetic equilibrium shape. It also allows for handling a non-steady magnetic equilibrium [6]—a critical point to model a full discharge including start-up and shut-down phases [7]. The ramp up and shut-down phases last about 30% of the full discharge time (a few seconds). These times are, however, long compared to characteristic turbulence times which are of the order of few micro-seconds. Indeed, the high-order accuracy of the spatial discretization allows for controlling the spurious numerical diffusion despite the strong anisotropy, as recently shown in [8] when increasing the order of interpolation p for a fixed spatial resolution. In this paper, we present an updated algorithm for solving 2D fluid–drift Braginskii equations in realistic tokamak geometries. The algorithm has been modified to handle nonlinear perpendicular diffusion terms with independent coefficients for each flow variable that allows a much more accurate description of the perpendicular transport related to turbulence. In addition, an h-refinement technique has been implemented to improve the numerical performance both in terms of memory and CPU time. This adaptive mesh refinement method can dynamically re-adjust the mesh locally according to error estimators based on the output data. The first results show it improves the global accuracy of the solution without using a global refinement of the mesh in the whole computational domain. The paper begins by introducing the physical model (Section 2) and the general features of the numerical algorithm (Section 3). The original developments are presented in Section 4 and Section 5 for the diffusive cross-field terms and the h-refinement technique, respectively. Concluding remarks and perspectives are summarized in Section 6.

2. Physical Model

The 2D computational domain mimics actual tokamaks with limiter or X-point and corresponds to the entire volume of plasma going from the core up to the wall as shown in Figure 1.

The magnetic field $\mathbf{B}$ is assigned including both closed flux surfaces in the center and open flux surfaces with field lines impacting the wall at the edge. These flux surfaces are separated by a magnetic field line in the poloidal called the separatrix in the poloidal cross-section. The strong difference of intensity between the toroidal and poloidal components $\Vert B_p\Vert <<\Vert B_t\Vert$ defines a privileged direction denoted as the parallel direction, with reference to the direction along the magnetic field lines. To take advantage of this flow anisotropy, the equations are projected along the magnetic field lines using the differential operator ${\mathbf{\nabla }}_{\Vert }=\mathit{b}·\mathbf{\nabla }$ and ${\mathbf{\nabla }}_{\perp }=\mathbf{\nabla }-\mathit{b}·\mathbf{\nabla }$, where $b=\frac{\mathit{B}}{\Vert \mathit{B}\Vert }$ is the unitary vector in the parallel direction.

2.1. Equations of the Model

The mathematical model relies on 2D fluid conservation equations based on Braginskii simplified closures [9]. Under some hypothesis and ordering detailed in Ref. [10], it corresponds to a standard model in the fusion community of advection diffusion equation that governs the transport of the mean plasma quantities as the density n is the parallel momentum $nu$, and the ion and electron total energy $E_i=\frac{3}{2}{k}_b{T}_i+\frac{1}{2}{m}_i{u}^2$ and $E_e=\frac{3}{2}{k}_b{T}_e$, respectively, with $m_i$ being the mass of the ion and $T_i$ and $T_e$ are the ion and electron temperatures, respectively. The conservation equations below correspond to a compressible adiabatic gas in the parallel direction and to an incompressible fluid in the perpendicular direction where turbulence process dominates. The system writes:

$${\partial }_{t}n+\nabla ·\left(nu\mathbf{b}\right)-\nabla ·\left(D{\nabla }_{\perp }n\right)=S_n$$

(1)

$${\partial }_t\left(m_{i}nu\right)+\nabla ·\left(m_{i}n{u}^2\mathbf{b}\right)+{\nabla }_{\Vert }\left(k_{b}n(T_e+T_i)\right)-\nabla ·\left(\mu {\nabla }_{\perp }\left(m_{i}nu\right)\right)=S_{\Gamma }$$

(2)

$$\begin{array}{cc}\hfill & {\partial }_t\left(\frac{3}{2}{k}_{b}n{T}_i+\frac{1}{2}{m}_{i}n{u}^2\right)+\nabla ·\left(\left(\frac{5}{2}{k}_{b}n{T}_i+\frac{1}{2}{m}_{i}n{u}^2\right)u\mathbf{b}\right)-nue{E}_{\Vert }\hfill \\ \hfill & -\nabla ·\left(\frac{3}{2}{k}_b(T_{i}D{\nabla }_{\perp }n+n{\chi }_i{\nabla }_{\perp }{T}_i)\right)-\nabla ·\left(-\frac{1}{2}{m}_i{u}^{2}D{\nabla }_{\perp }n+\frac{1}{2}{m}_i\mu n{\nabla }_{\perp }{u}^2\right)\hfill \\ \hfill & -\nabla ·\left(k_{\Vert i}{T}_i^{\frac{5}{2}}{\nabla }_{\Vert }{T}_i\mathbf{b}\right)+\frac{3}{2}\frac{k_{b}n}{{\hat{\tau }}_{ie}}(T_e-T_i)=S_{E_i}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & {\partial }_t\left(\frac{3}{2}{k}_{b}n{T}_e\right)+\nabla ·\left(\frac{5}{2}{k}_{b}n{T}_{e}u\mathbf{b}\right)+nue{E}_{\Vert }-\nabla ·\left(\frac{3}{2}{k}_b(T_{e}D{\nabla }_{\perp }n+n{\chi }_e{\nabla }_{\perp }{T}_e)\right)\hfill \\ \hfill & -\nabla ·\left(k_{\Vert e}{T}_e^{\frac{5}{2}}{\nabla }_{\Vert }{T}_e\mathbf{b}\right)-\frac{3}{2}\frac{k_{b}n}{{\hat{\tau }}_{ie}}(T_e-T_i)=S_{E_e}\hfill \end{array}$$

(4)

where $p_i$ and $p_e$ are the diagonal part of the ion and electron pressure stress tensor, and they are equal to $p_i=n{k}_b{T}_i$ and $p_e=n{k}_b{T}_e$$\left[{\mathrm{m}}^{-1}{\mathrm{s}}^{-2}\right]$, respectively. The constant diffusion coefficients that take into account the collisions transport and turbulent effects in the cross field direction are denoted $D,\mu ,{\chi }_i$ and ${\chi }_e$ for n, $nu$, $E_i$ and $E_e$, respectively. Their values are chosen as a compromise between estimations provided by theory or experimental measurements and numerical stability constraints. They are usually less or equal to $1{\mathrm{m}}^2{\mathrm{s}}^{-1}$. The terms $\left(k_{\Vert ,i}{T}_i^{5/2}\right)$ and $\left(k_{\Vert ,e}{T}_e^{5/2}\right)$ correspond to nonlinear parallel diffusions for ion and electron, respectively. The parallel diffusion coefficients depend on the mass of the species and are equal for the deuterium to $k_{\Vert ,i}=60\left[{\mathrm{Wm}}^{-1}{\mathrm{eV}}^{-7/2}\right]$ and $k_{\Vert ,e}=2000\left[{\mathrm{Wm}}^{-1}{\mathrm{eV}}^{-7/2}\right]$. The parameter ${\hat{\tau }}_{ie}$ is the relaxation time for the collisions coupling term between electrons and ions, $W=\frac{3}{2}\frac{k_{b}n}{{\hat{\tau }}_{ie}}(T_e-T_i)$. It is defined as:

$${\hat{\tau }}_{ie}=\frac{3\sqrt{2}}{e^4}\frac{{\epsilon }_0^2}{\Lambda }{\pi }^{\frac{3}{2}}\frac{m_i}{m_e}\sqrt{m_e}{e}^{\frac{3}{2}}\frac{T_e^{\frac{3}{2}}}{n}$$

where the Coulomb logarithm $\mathsf{\Lambda }=12$, the ionic mass $m_i=3.35·{10}^{-27}\left[\mathrm{kg}\right]$, the electronic mass $m_e=9.11\times {10}^{-31}\left[\mathrm{kg}\right]$, the vacuum permeability ${\epsilon }_0=8.85\times {10}^{-12}\left[\mathrm{C}{\mathrm{N}}^{-1}{\mathrm{m}}^{-1}\right]$ and the electron charge $e=1.60\times {10}^{-19}\left[\mathrm{C}\right]$. Finally, $S_n$,$S_{\Gamma }$, $S_{E_i},S_{E_e}$ correspond to sources’ terms.

2.2. Boundary Conditions

In the direction parallel to the magnetic field lines, the boundary conditions for the plasma are specific and correspond to the Bohm boundary conditions modeling plasma wall interactions [10]. They assume a parallel velocity of the plasma equal to or larger than the sound speed $c_s=\sqrt{\frac{k_b(T_e+T_i)}{m_i}}$ and leave free the density value at the wall that corresponds to ([6]):

$$\begin{array}{cc}\hfill & u\ge c_{s}if\mathbf{b}·\mathbf{n}>0\hfill \\ \hfill & u\le -c_{s}if\mathbf{b}·\mathbf{n}<0\hfill \end{array}$$

(5)

where n is the outer normal of the surface. For the electrons and ion energy equations, the Bohm conditions impose the parallel fluxes on the sheath transmission values, leading to:

$$\begin{array}{cc}\hfill & (n{E}_i+p_i)u-\frac{k_{\Vert i}}{m_i}{T}_i^{5/2}{\nabla }_{\Vert }{T}_i={\gamma }_{i}u{p}_i+\frac{1}{2}n{u}^3\hfill \\ \hfill & (n{E}_e+p_e)u-\frac{k_{\Vert e}}{m_i}{T}_e^{5/2}{\nabla }_{\Vert }{T}_e={\gamma }_{e}u{p}_e\hfill \end{array}$$

(6)

where ${\gamma }_i=2.5$ and ${\gamma }_e=4.5$. In the perpendicular direction to the magnetic field lines, homogeneous Neumann conditions are considered for all variables.

3. The Hybrid Discontinuous Galerkin Method

A specific Hybrid Discontinuous Galerkin (HDG) algorithm has been developed for many years [6,8,11,12], and implemented in the family of codes SOLEDGE3X [4], well-known in the international fusion community to efficiently address turbulent transport in different machines all around Europe. A complete description of the method is provided in Appendix A as well as in former papers [6,8,11,12]. In HDG, the system of Equations (1)–(4) is written in terms of conservative variables considering the vector:

$$\mathit{U}={\{U_1,U_2,U_3,U_4\}}^T={\{n,nu,n{E}_i,n{E}_e\}}^T$$

where the superscript ${\odot }^T$ stands for transpose. The discontinuous partition induces a two-step problem. In a first step, the set of conservative equations written in a weak formulation is solved element by element to express the discrete unknowns $\mathbf{U}(\mathbf{x},\mathbf{t})$ at the element nodes in terms of another approximation of the solution, called the trace solution $\hat{\mathbf{U}}$, which is defined on the borders of the element. In a second step, a global equation is set by imposing in a weak form the continuity of the fluxes across the borders of the elements to obtain $\hat{\mathbf{U}}$ in the whole mesh skeleton. Once $\hat{\mathbf{U}}$ is obtained, it is possible to recover the elementary solution $\mathbf{U}$ on each element using a local post processing. The introduction of this trace solution restricted to the skeleton of the mesh leads to a linear system of smaller size than in a classical discontinuous Galerkin method. The time discretization is fully implicit, and the nonlinear terms are linearized using a classic iterative Newton–Raphson method.

4. Implementation of Independent Nonlinear Diffusive Cross-Field Terms

In the model introduced in Section 2, the cross-field transport coefficients for $n,u,T_i,T_e$ play a fundamental role in the reliability of the solutions by modeling the perpendicular anomalous transport of particles and energy. Thus, their values directly impact the balance between the parallel and perpendicular transport which governs the plasma flow in the tokamak. With the implicit time integration scheme, the implementation of diffusion coefficients non equal for each flow variable is not straightforward. In this case, indeed, the expression of the coefficients as a function of conservative variables introduced additional nonlinear coupling between the equations as described thereafter. When assuming $D=\mu ={\chi }_i={\chi }_e$, the terms of the perpendicular dynamics in Equation (A2) depend only linearly on the unknown $\mathbf{Q}$ as follows:

$$-\nabla ·\left\{\begin{array}{c}D{\nabla }_{\perp }\left(n\right)\\ \mu {\nabla }_{\perp }\left(nu\right)\\ {\chi }_i{\nabla }_{\perp }\left(n{E}_i\right)\\ {\chi }_e{\nabla }_{\perp }\left(n{E}_e\right)\end{array}\right\}=-D_f\mathit{Q}+D_f\mathit{Q}\mathit{b}\otimes \mathit{b}$$

(7)

which represents the gradient of the conservative variable $\mathbf{U}$. When these coefficients are chosen to be non-equal, a nonlinear dependency occurs in the equations system written in conservative variables as:

$$\begin{array}{cc}\hfill & h_{\mathsf{\Gamma }}={\mathbf{Q}}_{\mathbf{t},\perp }·{\mathbf{W}}_{\mathsf{\Gamma }}={\mathbf{Q}}_{\mathbf{t},\perp }·\left[\begin{array}{c}(D-\mu )\frac{U_2}{U_1}\\ 0\\ 0\\ 0\end{array}\right]+\mu {\mathbf{Q}}_{\mathbf{2},\perp }\hfill \\ \hfill & h_{E_i}={\chi }_i{\nabla }_{\perp }\left(n{E}_i\right)={\mathbf{Q}}_{\mathbf{t},\perp }·{\mathbf{W}}_{E_i}={\mathbf{Q}}_{\mathbf{t},\perp }·\left[\begin{array}{c}(D-{\chi }_i)\frac{U_3}{U_1}(D-\mu )\frac{U_2^2}{U_1^2}\\ -(D-\mu )\frac{U_2}{U_1}\\ 0\\ 0\end{array}\right]+{\chi }_i{\mathbf{Q}}_{\mathbf{3},\perp }\hfill \\ \hfill & h_{E_e}={\mathbf{Q}}_{\mathbf{t},\perp }·{\mathbf{W}}_{E_e}={\mathbf{Q}}_{\mathbf{t},\perp }·\left[\begin{array}{c}(D-{\chi }_e)\frac{U_4}{U_1}\\ 0\\ 0\\ 0\end{array}\right]+{\chi }_e{\mathbf{Q}}_{\mathbf{4},\perp }\hfill \end{array}$$

(8)

Notice here that, for $D=\mu ={\chi }_i={\chi }_e$, Equation (8) is written as Equation (7). The linearization and integration of these nonlinear additional terms are detailed thereafter.

4.1. Linearization

The nonlinear terms of Equation (8) are written as:

$$\mathbf{h}(\mathbf{U},\mathbf{Q})={\mathbf{Q}}_{\mathbf{t}}\mathbf{W}\left(\mathbf{U}\right)$$

(9)

that linearize according to the Formula (A15) as:

$$\begin{array}{cc}\hfill \mathbf{h}({\mathbf{U}}^k,{\mathbf{Q}}^k)& =\mathbf{h}({\mathbf{U}}^{k-1},{\mathbf{Q}}^{k-1})+\frac{d}{dϵ}\mathbf{h}({\mathbf{U}}^{k-1}+ϵ\mathbf{dU},{\mathbf{Q}}^{k-1}+ϵ\mathbf{dQ}){|}_{ϵ=0}+\hfill \\ \hfill & +O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)=\hfill \\ \hfill & ={\mathbf{Q}}^{k-1}\mathbf{W}\left({\mathbf{U}}^{k-1}\right)+\frac{d}{dϵ}\left(({\mathbf{Q}}^{k-1}+ϵ\mathbf{dQ})\mathbf{W}({\mathbf{U}}^{k-1}+ϵ\mathbf{dU})\right){|}_{ϵ=0}+\hfill \\ \hfill & +O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)=\hfill \\ \hfill & ={\mathbf{Q}}^{k-1}\mathbf{W}\left({\mathbf{U}}^{k-1}\right)+\mathbf{dQ}\mathbf{W}\left({\mathbf{U}}^{k-1}\right)+{\mathbf{Q}}^{k-1}\frac{d\mathbf{W}}{d\mathbf{U}}{|}_{k-1}\mathbf{dU}+\hfill \\ \hfill & +O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)=\hfill \\ \hfill & ={\mathbf{Q}}^{k-1}\mathbf{W}\left({\mathbf{U}}^{k-1}\right)+{\mathbf{Q}}^k\mathbf{W}\left({\mathbf{U}}^{k-1}\right)-{\mathbf{Q}}^{k-1}\mathbf{W}\left({\mathbf{U}}^{k-1}\right)+\hfill \\ \hfill & +{\mathbf{Q}}^{k-1}\frac{d\mathbf{W}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k-{\mathbf{Q}}^{k-1}\frac{d\mathbf{W}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^{k-1}+O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)=\hfill \\ \hfill & ={\mathbf{Q}}^k\mathbf{W}\left({\mathbf{U}}^{k-1}\right)+{\mathbf{Q}}^{k-1}\frac{d\mathbf{W}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k+O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)\hfill \end{array}$$

(10)

where $\mathbf{dU}$ and $\mathbf{dQ}$ have been replaced by $\mathbf{dU}={\mathbf{U}}^k-{\mathbf{U}}^{k-1}$ and $\mathbf{dQ}={\mathbf{Q}}^k-{\mathbf{Q}}^{k-1}$. Then, ${\mathbf{h}}_{\mathsf{\Gamma }}({\mathbf{U}}^k,{\mathbf{Q}}^k)$, ${\mathbf{h}}_{E_i}({\mathbf{U}}^k,{\mathbf{Q}}^k)$ and ${\mathbf{h}}_{E_e}({\mathbf{U}}^k,{\mathbf{Q}}^k)$ linearize as:

$$\begin{array}{cc}\hfill & {\mathbf{h}}_{\mathsf{\Gamma }}({\mathbf{U}}^k,{\mathbf{Q}}^k)={\mathbf{Q}}^k{\mathbf{W}}_{\mathsf{\Gamma }}\left({\mathbf{U}}^{k-1}\right)+{\mathbf{Q}}^{k-1}\frac{d{\mathbf{W}}_{\mathsf{\Gamma }}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k+O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)\hfill \\ \hfill & {\mathbf{h}}_{E_i}({\mathbf{U}}^k,{\mathbf{Q}}^k)={\mathbf{Q}}^k{\mathbf{W}}_{E_i}\left({\mathbf{U}}^{k-1}\right)+{\mathbf{Q}}^{k-1}\frac{d{\mathbf{W}}_{E_i}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k+O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)\hfill \\ \hfill & {\mathbf{h}}_{E_e}({\mathbf{U}}^k,{\mathbf{Q}}^k)={\mathbf{Q}}^k{\mathbf{W}}_{E_e}\left({\mathbf{U}}^{k-1}\right)+{\mathbf{Q}}^{k-1}\frac{d{\mathbf{W}}_{E_e}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k+O({\mathbf{dU}}^2,{\mathbf{dQ}}^2)\hfill \end{array}$$

(11)

where:

$$\begin{array}{cc}\hfill \frac{d{\mathbf{W}}_{\mathsf{\Gamma }}}{d\mathbf{U}}& =\left[\begin{array}{cccc}-(D-\mu )\frac{U_2}{U_1^2}& (D-\mu )\frac{1}{U_1}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\hfill \\ \hfill \frac{d{\mathbf{W}}_{E_i}}{d\mathbf{U}}& =\left[\begin{array}{cccc}-\left((D-{\chi }_i)\frac{U_3}{U_1^2}+2(D-\mu )\frac{U_2^2}{U_1^3}\right)& 2(D-\mu )\frac{U_2}{U_1^2}& (D-{\chi }_i)\frac{1}{U_1}& 0\\ (D-\mu )\frac{U_2}{U_1^2}& -(D-\mu )\frac{1}{U_1}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\hfill \\ \hfill \frac{d{\mathbf{W}}_{E_e}}{d\mathbf{U}}& =\left[\begin{array}{cccc}-(D-{\chi }_e)\frac{U_4}{U_1^2}& 0& 0& (D-{\chi }_e)\frac{1}{U_1}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\hfill \end{array}$$

(12)

Defining now:

$$\begin{array}{cc}\hfill & {\mathbf{h}}_{\mathsf{\Gamma }}^{\mathbf{U}}={\mathbf{Q}}^{k-1}\frac{d{\mathbf{W}}_{\mathsf{\Gamma }}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k;{\mathbf{h}}_{\mathsf{\Gamma }}^{\mathbf{Q}}={\mathbf{Q}}^k{\mathbf{W}}_{\mathsf{\Gamma }}\left({\mathbf{U}}^{k-1}\right)\hfill \\ \hfill & {\mathbf{h}}_{E_i}^{\mathbf{U}}={\mathbf{Q}}^{k-1}\frac{d{\mathbf{W}}_{E_i}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k;{\mathbf{h}}_{E_i}^{\mathbf{Q}}={\mathbf{Q}}^k{\mathbf{W}}_{E_i}\left({\mathbf{U}}^{k-1}\right)\hfill \\ \hfill & {\mathbf{h}}_{E_e}^{\mathbf{U}}={\mathbf{Q}}^{k-1}\frac{d{\mathbf{W}}_{E_e}}{d\mathbf{U}}{|}_{k-1}{\mathbf{U}}^k;{\mathbf{h}}_{E_e}^{\mathbf{Q}}={\mathbf{Q}}^k{\mathbf{W}}_{E_e}\left({\mathbf{U}}^{k-1}\right)\hfill \end{array}$$

(13)

the split momentum diffusion terms are as follows:

$$\begin{array}{cc}\hfill & {\mathbf{h}}_{\mathsf{\Gamma }}={\mathbf{h}}_{\mathsf{\Gamma }}^{\mathbf{U}}+{\mathbf{h}}_{\mathsf{\Gamma }}^{\mathbf{Q}}\hfill \\ \hfill & {\mathbf{h}}_{E_i}={\mathbf{h}}_{E_i}^{\mathbf{U}}+{\mathbf{h}}_{E_i}^{\mathbf{Q}}\hfill \\ \hfill & {\mathbf{h}}_{E_e}={\mathbf{h}}_{E_e}^{\mathbf{U}}+{\mathbf{h}}_{E_e}^{\mathbf{Q}}\hfill \end{array}$$

(14)

The superscripts $\mathbf{U}$ and $\mathbf{Q}$ stand for the terms whose unknowns are ${\mathbf{U}}^k$ and ${\mathbf{Q}}^k$, respectively. These terms must be now incorporated into the matrices of the discrete linear system.

4.2. The New Discrete Linear System

From Equation (8), the diffusion terms can be actually written as the sum of two terms as ${\mathbf{S}}_d+D_f\mathbf{Q}$ with ${\mathbf{S}}_d=0$ for $D=\mu ={\chi }_i={\chi }_e$. To incorporate the new term ${\mathbf{S}}_d$ into the linear system, ${\mathbf{S}}_d$ is first written in the matrix form as:

$${\mathbf{S}}_d={\mathbf{S}}^{\mathbf{U}}+{\mathbf{S}}^{\mathbf{Q}}=\left[\begin{array}{c}0\\ {\mathbf{h}}_{\mathsf{\Gamma }}^{\mathbf{U}}\\ {\mathbf{h}}_{E_i}^{\mathbf{U}}\\ {\mathbf{h}}_{E_e}^{\mathbf{U}}\end{array}\right]+\left[\begin{array}{c}0\\ {\mathbf{h}}_{\mathsf{\Gamma }}^{\mathbf{Q}}\\ {\mathbf{h}}_{E_i}^{\mathbf{Q}}\\ {\mathbf{h}}_{E_e}^{\mathbf{Q}}\end{array}\right]$$

(15)

Focusing on the second equation of the system (A12), the local problem is as follows:

$$\begin{array}{cc}\hfill & {\left(\mathit{v},{\partial }_t\mathit{U}\right)}_{{\Omega }_i}-{\left(\mathbf{\nabla }\mathit{v},\mathit{F}-D_f\mathit{Q}+D_f\mathit{Qb}\otimes \mathit{b}-{\mathit{F}}_t\right)}_{{\Omega }_i}+\hfill \\ \hfill & +{〈\mathit{v},\left(\hat{\mathit{F}}-D_f\hat{\mathit{Q}}+D_f\hat{\mathit{Q}}\mathit{b}\otimes \mathit{b}-{\hat{\mathit{F}}}_t\right)\mathit{n}〉}_{\partial {\Omega }_i}+{\left(\mathit{v},{\mathit{f}}_{E_{\Vert }}\right)}_{{\Omega }_i}+{\left(\mathit{v},{\mathit{f}}_{E_{EX}}\right)}_{{\Omega }_i}-\hfill \\ \hfill & -{\left(\mathit{v},\mathit{g}\right)}_{{\Omega }_i}-{(\mathbf{\nabla }\mathit{v},-{\mathit{S}}^{\mathit{U}}+{\mathit{S}}^{\mathit{U}}\mathit{b}\otimes \mathit{b})}_{{\Omega }_i}-{(\mathbf{\nabla }\mathit{v},-{\mathit{S}}^{\mathit{Q}}+{\mathit{S}}^{\mathit{Q}}\mathit{b}\otimes \mathit{b})}_{{\Omega }_i}\hfill \\ \hfill & +{〈\mathit{v},\left(-{\mathit{S}}^{\hat{\mathit{U}}}+{\mathit{S}}^{\hat{\mathit{U}}}\mathit{b}\otimes \mathit{b}\right)·\mathit{n}〉}_{\partial {\Omega }_i}+{〈\mathit{v},\left(-{\mathit{S}}^{\hat{\mathit{Q}}}+{\mathit{S}}^{\hat{\mathit{Q}}}\mathit{b}\otimes \mathit{b}\right)·\mathit{n}〉}_{\partial {\Omega }_i}={\left(\mathit{v},\mathit{s}\right)}_{{\Omega }_i}\hfill \end{array}$$

(16)

Using the convention introduced in Appendix B.3, the terms with the unknown $\mathbf{U}$ are inserted in the matrix of the local problem $A_{uu}$ while the terms with the unknown $\mathbf{Q}$ are inserted in $A_{uq}$. Then, in the discrete local problem, the new matrices are as follows:

$$\begin{array}{cc}\hfill & A_{uu}\Rightarrow A_{uu}-{(\mathbf{\nabla }\mathit{v},-{\mathit{S}}^{\mathit{U}}+{\mathit{S}}^{\mathit{U}}\mathit{b}\otimes \mathit{b})}_{{\Omega }_i}+{〈\mathit{v},\left(-{\mathit{S}}^{\hat{\mathbf{U}}}+{\mathit{S}}^{\hat{\mathit{U}}}\mathit{b}\otimes \mathit{b}\right)·\mathit{n}〉}_{\partial {\Omega }_i}\hfill \\ \hfill & A_{uq}\Rightarrow A_{uq}-{(\mathbf{\nabla }\mathit{v},-{\mathit{S}}^{\mathit{Q}}+{\mathit{S}}^{\mathit{Q}}\mathit{b}\otimes \mathit{b})}_{{\Omega }_i}+{〈\mathit{v},\left(-{\mathit{S}}^{\hat{\mathit{Q}}}+{\mathit{S}}^{\hat{\mathit{Q}}}\mathit{b}\otimes \mathit{b}\right)·\mathit{n}〉}_{\partial {\Omega }_i}\hfill \end{array}$$

(17)

Let us notice that the new matrices are just related to the second equation of the system (A12), and the changes are limited to the local element by element problem. It is worth observing that, in the formulation of the global problem, the perpendicular gradient term is included in the imposition of the normal fluxes at the element boundary in Equation (A14). Moreover, it is also present in the flux vector that defines the Bohm boundary condition where the normal gradient is imposed as equal to 0. In addition, in this case, the contribution of the split diffusion term has to be considered, and the matrices for the assembling of the global problem $A_{ll},A_{lq}$ are modified in the same way by the additional terms of Equation (17).

4.3. Code Verification

The Method of the Manufactured Solution (MMS) [6] is used to verify the code with the new formulation (Equation (16)). The transport coefficients are specially set all differently from each other: $D=0.1,\mu =0.2,{\chi }_i=0.3,{\chi }_e=0.4{\mathrm{m}}^2/\mathrm{s}$. The following analytical solution is used with ${\omega }_x={\omega }_y=1$:

$$\begin{array}{cc}\hfill & n=2+sin\left(2\pi {\omega }_{x}x\right)sin\left(2\pi {\omega }_{y}y\right);E_i=20+cos\left(2\pi {\omega }_{x}x\right)sin\left(2\pi {\omega }_{y}y\right)\hfill \\ \hfill & u=cos\left(2\pi {\omega }_{x}x\right)cos\left(2\pi {\omega }_{y}y\right);E_e=10-sin\left(2\pi {\omega }_{x}x\right)cos\left(2\pi {\omega }_{y}y\right)\hfill \end{array}$$

(18)

Results of convergence plotted in Figure 2 show the expected theoretical rate of convergence in $p+1$, and thus the correct implementation of the non-equal diffusion coefficients in the solver.

4.4. Example of Simulation in the WEST Tokamak

In order to show the new capability of the code to run with different cross-field coefficients, Equations (1)–(4) are resolved in the WEST geometry (Figure 1). We assume ${\chi }_i={\chi }_e$ as it is usual in current computations of the literature, in agreement with experimental measurements carried out at the tokamak cross-section midplane [13]. For simplicity here, we choose $\nu ={\chi }_i={\chi }_e=1$, and only D, the particles diffusion, is varied in a short range between 1 and 0.6 to avoid the use of too fine meshes. A mesh of 15,591 elements with $p=6$-elements is used. These steady state simulations require a run-time of about 40 min each on 32 cpu.

In Figure 3, the 2D contours for all flow variables are compared to ones obtained with a former version of the algorithm where all cross-field coefficients had to be equal to 1. The two solutions globally agree showing that the new version of the solver is able to provide 2D plasma equilibrium in realistic geometry. As expected, the solution at $D=0.6$, however, shows some differences. The contours are sharper contours in particular at the X-point, and the density is higher in the core, of about a factor 1.6, since less matter is allowed to diffuse from it. On the parallel Mach number, the tongue of positive velocity extends towards the top to the same extent while slightly decreasing its width, as is to be expected for a lower density diffusion. Moreover, the parallel Mach number is higher at the X-points. Regarding the ion and electron temperatures now, they are globally lower for $D=0.6$, meaning that, for this value of density diffusion, the plasma in the core has a higher density but lower temperatures.

5. Spatial Adaptivity

Plasma solutions of interest for tokamak operation may exhibit large gradients both in the radial and parallel flow directions when targeting realistic conditions for the simulations, corresponding generally to small values of the cross-field diffusion coefficients [4]. This routinely leads to demanding requirements on the local spatial resolution of the mesh. In practice, failure to design a mesh that accommodates these resolution requirements results in aliasing errors in some elements of the mesh that may lead to divergence of the Newton–Raphson iterations during the convergence toward the steady state solution. With the objective of enabling a robust numerical modeling of plasma transport in the edge, an adaptive h-refinement has been implemented. The h-refinement method is based here on an oscillation indicator to target flow regions with steep gradients or discontinuities inside the domain of computation. The element size is then optimized by imposing iterative, local mesh refinements in these flow regions while keeping a coarse mesh elsewhere [14].

5.1. Refinement Process Strategy

Experience in the computation of steady-state solutions of plasma transport in the edge has led to the emergence of a strategy combining Newton–Raphson iterations, with progressive lowering of cross-field diffusion coefficients in Equations (1)–(4) in order to reach the desired value imposed by the simulation of tokamak operation (around $1$m${}^2·$s${}^{-1}$ or lower). The Newton–Raphson iterative process is led to convergence for each value of the diffusion, and the obtained solution is used as an initial condition of the Newton–Raphson iterations for the next smaller value of diffusion. The h-refinement is adopted for optimizing the mesh design, refining each element on which oscillations are detected. The procedure is stopped when the iterations reach the desired level of accuracy for the targeted diffusion coefficients’ values.

The whole process can be thus summarized as:

• Initialize the calculation with a rather coarse mesh and large values of cross-field diffusion coefficients;
• Convergence to the steady solution using Newton–Raphson iterations;

  -

  if convergence, computations are going on, lowering diffusion;

  -

  if non convergence, the refinement procedure is started;

  *

  Interpolation of the solution on the new mesh locally refined;

  *

  Convergence to the steady solution using Newton–Raphson iterations;

• Stop when diffusion coefficients reach the target values.

The mesh refinement is performed using the open-source software Mmg [15,16]. It uses a map of elemental size in which the desired element size on each vertex must be precise. This current elemental size can be defined on each vertex of an existing mesh [17] using the elemental areas $\left\{\right|{\Omega }_k\left|\right\}$ as follows:

$$h_j=\frac{{\sum }_{i\in S_j}\left|{\Omega }_i\right|{\tilde{h}}_i}{{\sum }_{i\in S_j}\left|{\Omega }_i\right|}$$

(19)

in which $S_j$ denotes the set of element indices having node j as a vertex. At the iteration n of the refinement process, a basic and straightforward formula provides a guess of the desired mesh size at the next iteration, on the element j where oscillations are detected, using the expression:

$$h_{target,j}^{(n+1)}=\frac{\left(h_j^{\left(n\right)}\right)}{\alpha }$$

(20)

where $\alpha$ ($\alpha >1$) is a control parameter to tune in order to perform the refinement. After several tests, the optimal value $\alpha =2$ has been found. In this process, it is obvious that the mesh size is decreased locally and in an isotropic way. The possibility to coarsen the mesh has not been taken into account here because in the present configurations the initial meshes are already very coarse.

The efficiency of such a refinement strategy is mainly based on the choice of a suitable mesh refinement estimator. This estimator must be well-calibrated to avoid unnecessary costly over-refinements or, on the contrary, to keep spurious undetected oscillations in the solution. Here, the estimator can be more considered as an indicator uniquely able to identify spurious oscillations in the solution, related to unresolved steep gradients or discontinuities.

5.2. Oscillation-Based Error Indicator

Adaptive mesh refinement is usually considered to converge to a numerical solution with a desired accuracy whilst using a minimal number of degrees of freedom. Adaptive mesh refinement is especially appealing in DG and HDG discretizations using $hp$-refinement as it warrants exponential convergence with the number of degrees of freedom [18]. The present refinement strategy is not driven by an accuracy criterion, but by a stability criterion to ensure the convergence of Newton–Raphson iterations towards the steady solution of Equations (1)–(4). This strategy is based upon the observation that lack of convergence mostly stems from locally insufficient spatial resolution leading to aliasing errors. These errors deteriorate the convergence of the implicit solver and the global accuracy of the solution, and even more may lead to the divergence of the computations. This problem can be overcome by increasing the resolution locally to enhance the precision of the interpolation and to damp spurious oscillations. Usually, the estimators are based on the output data of the simulation [19] to detect oscillations. The technique is inspired from shock-capturing techniques [6], although here the quantity evaluated is an oscillation rather than a discontinuity in the solution. We use a simple sensor $S_k$, defined on each element with index k defined as a function of the parallel velocity u. For a computation with a polynomial approximation of order p, this sensor consists of the norm of the local contribution of order p, normalized by the norm of the full solution on the element. It is thus defined as

$$S_k=\frac{{(u-\hat{u},u-\hat{u})}_{{\Omega }_k}}{{(u,u)}_{{\Omega }_k}}$$

(21)

where u is the solution of order p, and $\hat{u}$ is the projection of the modal expansion on the space of polynomials of order $p-1$.

5.3. Results

For simplicity, a reduced 2D fluid isothermal model is derived from Equations (1) and (2) to solve the density n and the parallel momentum $nu$ in a realistic WEST geometry (see in Ref. [6]). As in the complete model, Bohm boundary conditions are prescribed in the parallel direction to the magnetic field lines. Although simpler, this reduced model allows for evaluating most of the numerical issues. It takes into account the anisotropy in the flow dynamics between the parallel and perpendicular directions, and the balance between the transport in the two flow directions is simply modulated by varying the diffusion D ($D=\mu$). Lowering D makes the parallel transport dominant that can be very demanding for the solver, particularly in the present configuration where the mesh is not aligned along the magnetic field lines [8]. Thus, the mesh has to be successively refined when decreasing D to converge toward a plasma equilibrium as already shown in [6] for uniform meshes. In addition, there is also a geometrical complexity with a magnetic equilibrium with two X-points as well as a tokamak wall with sharp edges and corners as well as small cavities around Figure 1. This is thus an attractive configuration to test the local h-refinement technique proposed in this work.

Calculations are performed here using different meshes, automatically designed by the adaptive procedure described above. Only $p=4$-polynomials are considered.

Typical contour plots of the density and the parallel Mach number are shown in Figure 4 in the WEST poloidal cross-section for $D=0.83$m${}^2·$s${}^{-1}$. The large scale flow prediction shows very similar trends with respect to the literature [6]. The density is maximum at the core boundary, where the Dirichlet condition $n=1$ is applied, and rapidly decreased to low values in the boundary layer, called the scrape-off layer (SOL), beyond the separatrix. The parallel Mach number $u/c_s$ shows positive and negative Mach number regions, and a flow reversal around the midplane. As expected from theoretical analysis [20] and from numerical investigations [21], the solution exhibits transitions to supersonic flows in the vicinity of both divertor legs.

To show the adaptive h-refinement process, Figure 5 shows the grid refinement at three successive steps for a diffusion $D=2.63$m${}^2·$s${}^{-1}$. For each mesh, the oscillations of the solution detected by the estimator are emphasized. Starting with a relatively coarse mesh, the results show that the refinement process reduces the elements size only in the flow regions where oscillations are detected. Accordingly, the number of elements increases progressively in the poloidal cross-section with $N_e=1192$, 2221 and 2554 but much less than if a uniform refinement had been considered. At the final step, oscillations are totally damped by the increased resolution around, and the solver converges. With this procedure, the mesh is automatically designed with a number of degree of freedom which is close to being optimal.

As mentioned above, lowering the diffusion coefficient toward realistic values challenges the numerical solver by making the parallel transport dominant with oscillations if the resolution is not fine enough. The mesh must be then automatically adapted for each value of the cross-field diffusion to ensure the convergence of the algorithm. Once the solution converged, the diffusion is lowered again, and a new mesh is generated with an optimal design. This is shown in Figure 6 where the diffusion coefficient is progressively lowered by a factor of 100, and the mesh is automatically refined accordingly.

As soon as the mesh is fine enough, we can clearly expect to save on the time needed to converge. However, it is not straightforward to quantify precisely this saving. We have first compared the simulation times to convergence when lowering the diffusion coefficients (Figure 6) between simulations using the automative adaptive refinement procedure and simulations performed with a unique mesh for each value of the diffusion, corresponding to the most refined mesh designed during the automative procedure. Results are reported in

Table 1. CPU times in seconds to convergence depending on the diffusion coefficients and the corresponding meshes for simulations with and without the h–refinement technique. $N_e$ is the number of elements, nDOF is the number of degrees of freedom for $p=4$–polynomials. Without h–refinement, a unique mesh is used for each value of the diffusion, corresponding to the most refined mesh designed during the automative procedure.

D $({\mathbf{m}}^2·{\mathbf{s}}^{-1})$	${\mathit{N}}_{\mathit{e}}$	nDOF	h-Refinement (Time (s))	No h-Refinement (Time (s))	Time Saving $(\%)$
26.31	388	5820	13.96	14.20	+2%
8.32	1192	17,880	30.84	35.57	+13%
2.63	3219	49,590	126.51	131.67	+4%
0.83	6066	114,120	203.75	280.78	+28%
0.26	10,032	150,480	285.61	366.41	+23%

below, and show a savings of time up to $28\%$ as D is strictly smaller than $D=0.83{\mathrm{m}}^2·{\mathrm{s}}^{-1}$. Let us remind that target values for tokamak operation simulations are smaller than $D=1{\mathrm{m}}^2·{\mathrm{s}}^{-1}$. As expected, when the number of elements in the mesh is not high enough, there is no saving, and there is even an additional cost due to the time needed by the algorithm to design the mesh which is naturally not taken into account in the second set of simulations.

As additional information, we have compared times to converge at $D=2.63$m${}^2·$s${}^{-1}$ using the adaptive procedure described above (Figure 5) and a uniform mesh with element size equal to the size of the smallest element provided by the adaptive procedure. Doing that, the respective meshes are composed by 3219 and 98,372 elements, respectively. The corresponding times to converge are respectively equal to 126.51 s and 5488 s, which corresponds to an increase of a factor 43 when using a uniform mesh. Naturally, this is only informative since uniform meshes are rarely used, but the time to accurately design a mesh for each value of the diffusion coefficient when lowering it can be long and impossible to estimate depending on the user’s skills. The automatic design of the mesh, which does not require any adjustment by hand during the iterative process is clearly a great advantage of this procedure.

6. Conclusions

This paper presents a high-order solver based on the Hybrid Discontinuous Galerkin method to perform plasma simulations in tokamak. It solves a 2D fluid transport model for the density, parallel momentum, and the total energy for a deuterium plasma. This model is relevant with those currently implemented in fluid codes used in the fusion community. The main features of this solver are the use of unstructured meshes together with a high-order spatial approximation which allows for misaligning the discretization from the magnetic field, unlike what is required in lower-order numerical schemes in order to control the spurious numerical diffusion due to the strong anisotropy of the flow. Thus, realistic tokamak wall geometries as well as magnetic equilibriums of complex shape and eventually unsteady can be accurately treated.

The code development is still in progress. In this paper, we have generalized the treatment of the cross-field diffusion terms. The possibility to handle diffusion coefficients chosen independently for each variable is a real improvement in the modeling of the cross-field turbulent transport. To progress toward better numerical performance, the first steps of an h-refinement technique have been introduced to optimize the mesh design and save on CPU time and memory. Involving an error indicator based on spurious oscillations related to aliasing error, the mesh is refined locally and automatically around steep gradients of the solution that allows for efficiently damping the oscillations. This technique allows for saving CPU time, and clearly improves the stability and the robustness of the algorithm.

This work is thus a step forward in the development of a very efficient and accurate numerical solver able to solve a 2D transport fluid model in realistic tokamak configurations relevant for the operation.