A high-order non field-aligned approach for the discretization of strongly anisotropic diffusion operators in magnetic fusion

Highlights

• A non-aligned approach is presented for solving highly anisotropic problems.

• The HDG stabilization tensor is designed to provide convergence at each anisotropy level.

• High-order computations are used to reduce the numerical diffusion.

• An explicit measure of the numerical diffusion is proposed.

• A comparison with two finite-difference schemes based on non-aligned discretizations is presented.

Abstract

In this work we present a hybrid discontinuous Galerkin scheme for the solution of extremely anisotropic diffusion problems arising in magnetized plasmas for fusion applications. Unstructured meshes, non-aligned with respect to the dominant diffusion direction, allow an unequalled flexibility in discretizing geometries of any shape, but may lead to spurious numerical diffusion. Curved triangles or quadrangles are used to discretize the poloidal plane of the machine, while a structured discretization is used in the toroidal direction. The proper design of the numerical fluxes guarantees the correct convergence order at any anisotropy level. Computations performed on well-designed 2D and 3D numerical tests show that non-aligned discretizations are able to provide spurious diffusion free solutions as long as high-order interpolations are used. Introducing an explicit measure of the numerical diffusion, a careful investigation is carried out showing an exponential increase of this latest with respect to the non-alignment of the mesh with the diffusion direction, as well as an exponential decrease with the polynomial degree of interpolation. A brief assessment of the method with respect to two finite-difference schemes using non-aligned discretization, but classically used in fusion modeling, is also presented.

Program summary

Program Title: Laplace-HDG (Laplace Hybrid Discontinuous Galerkin)

CPC Library link to program files: http://dx.doi.org/10.17632/c3dhycyvj8.1

Licensing provisions: GPLv3

Programming language: Fortran 95

Nature of problem: Anisotropic Laplace problem in 2D with Dirichlet boundary conditions

Solution method: Hybrid discontinuous Galerkin scheme

Introduction

Anisotropic diffusion problems arise in a variety of applications in physics and engineering, such as, for example, hydrodynamic transport [1], image processing [2], [3], [4], diffusion in porous media [5], electromagnetic propagation [6], biomedical applications [7], [8], etc.

The typical mathematical model for the convection–diffusion of a scalar quantity $u$ in ${\mathbb{R}}^3$ writes as:

$$\frac{\partial u}{\partial t}+\mathit{a}\cdot \mathit{\nabla }u-\mathit{\nabla }\cdot \left(\mathcal{K}\mathit{\nabla }u\right)=f$$

where $\mathcal{K}$ is an anisotropic diffusion tensor, $\mathit{a}$ is the convection velocity and $f$ is a source term. In this work we focus on magnetized plasmas for fusion applications, where the scale separation induced by the magnetic field produces conduction coefficients that are several order of magnitudes larger in the parallel direction than in the perpendicular one. In this context, the usual strategy to solve problem (1) is based on finite-difference schemes with structured meshes aligned with the magnetic lines. The separation of the dynamics in the two directions is therefore guaranteed by the space discretization, and no spurious cross-field diffusion occurs, see [9]. Flux aligned coordinates are still used in the fusion community, see for example [10]. However, they require simplifications on the magnetic field topology: in fact, often the magnetic lines wind around a toroidal surface never closing on themselves (the magnetic surface are, in general, “irrational”). Therefore, the possibility of using completely aligned meshes in general tokamak simulations is excluded.

Hence, it is more and more common in the fusion community the use of toroidal coordinates for the spatial discretization of (1) in axisymmetric geometries mimicking the topology of fusion devices such as tokamaks (neglecting the use of non-symmetric components such as toroidal limiters or antennas). This entails the need of two kinds of discretizations: a poloidal one and a toroidal one, see Fig. 1. While the poloidal discretization carries the burden of describing the geometry of the reactor chamber and the plasma facing components, the toroidal one can be simply thought of a periodic extrusion in the toroidal direction of the poloidal discretization. Within this approach, field aligned coordinates only in the poloidal plane are still widely used for both 2D and 3D computations. However, a precise geometry description is cumbersome with this technique: in [11], for example, the tokamak chamber is described only until the last flux surface before the wall. In [12], a penalization technique is used to model the real chamber geometry. Unfortunately, penalization techniques are usually only first order accurate with respect to the geometry definition and also the imposition of the boundary conditions. Moreover, even if the use of field aligned coordinates in the poloidal plane allows to threat complicated magnetic topologies and X-point configurations, using domain decomposition techniques to avoid the problem of the singularity at the X-point, it still introduces low resolution areas due to the flux expansion at the singularity. Finally, it also prevents the simulation of evolving magnetic equilibria, which would require expensive on the fly re-meshing at each time step.

In the objective of relaxing the need of aligning the grid, several schemes based on non-aligned meshes have been proposed in the literature, mostly in the context of finite-differences (see details in [13]). Schemes based on non-aligned discretizations [14], [15] improve the computation of the fluxes by a proper choice of the stencil with respect to a naive discretization, or rely on high-order finite-difference schemes to reduce the numerical diffusion, such as in [16]. On the other hand, schemes based on aligned discretizations [17], [18], [19], [20] employ interpolations on the magnetic field lines to perform finite-difference discretizations along the anisotropy directions.

In the perpendicular plane, such aligned discretizations seem so far to be limited to structured meshes for simplicity and efficiency of implementation. This tends thus to complicate the description of realistic geometries and the imposition of the boundary condition. Let us notice however that such aligned discretizations in the parallel direction can be associated to unstructured mesh in the perpendicular plane, see for example in Ref. [21].

Switching to general non-structured non-aligned meshes (of triangles or quads) allows to exploit a huge flexibility in terms of geometry description and local refinement. An intermediate approach is proposed in [22], where a technique aimed at aligning a generic triangular mesh with respect to a given flow is developed in the framework of finite-volumes. The technique is shown to reduce the numerical diffusion errors for 2D computations. Triangular grid alignment in the framework of magnetic fusion is an active field, see for example [23], [24], and it can be seen as another opportunity to reduce the numerical errors and the computational cost, for simulations where the magnetic equilibrium is fixed and therefore the mesh may be designed to mimick the elongated structures appearing in the field direction. For this reason, a discussion on the definition of alignment for triangular meshes is proposed in Section 5. However, in this work the focus is put in reducing the diffusion for non-aligned grids by means of increasing the polynomial order of interpolation. Some effort in this sense has already been done. In [25], a spectral element method is used to solve a heat diffusion problem in plasmas, showing that high-order elements yield a given accuracy with less total degrees of freedom than lower-order elements. In [26], a second and fourth order finite volume schemes are tested for solving a temporal evolution diffusion problem, revealing that the fourth order scheme is more efficient than the second order one for reducing the numerical perpendicular diffusion. High-order finite-volumes schemes however suffer from large stencils, preventing the use of unstructured meshes. In fact, all the tests in [26] are performed using structured meshes. A well-established solution to that is to localize the high-order interpolation inside each element: the discontinuous Galerkin (DG) scheme provides this approach. In [27], three DG schemes are used to discretize the parallel diffusion operator, showing exponential decrease of the numerical diffusion as the interpolation order is increased.

Introduced by Cockburn in [28], the hybrid discontinuous Galerkin scheme (HDG) retains the interesting characteristic of classic DG schemes, such as stability and local conservation, while reducing the size of the assembled matrix produced by the spatial discretization thanks to the hybridization of the primal unknown, which allows to express the nodal values of the solution in terms of the nodal values of the trace unknown, defined on the element borders. In [29], a HDG scheme is proposed for the simulation of anisotropic and non-homogeneous media arising in the study of flow through porous media. The results are encouraging for mild anisotropy values considered, and are extended here to extreme anisotropies arising in fusion applications. In [30], a first application of the HDG scheme to an isothermal plasma-transport model in tokamaks is presented. The scheme is extended in [31] to a non-isothermal model containing parallel diffusion operators. Encouraging results in [31] stimulate the deeper investigation proposed in this manuscript on the discretization of parallel diffusion operators with non-aligned high-order schemes.

The following Section 2 introduces the anisotropic convection–diffusion problem to be solved. Following our recent work [30], a HDG discretization suited to handle highly anisotropic diffusion operators is proposed in Section 3. A discussion on the design of the stabilization parameter is presented in Section 4. In Section 5, the definition of alignment for triangular meshes is considered, with particular emphasis on the importance of high-order computations on the control of the numerical diffusion in case of non-aligned meshes. In the following three sections is investigated the impact of the mesh alignment and the element type (quadrangles or triangles) on the solution accuracy, using 2D and 3D well-designed numerical tests. Section 6 concerns the discretization of an annulus domain mimicking the poloidal section of a tokamak. In Section 7, the impact of the disalignment in the toroidal discretization is investigated and a brief assessment of the method is proposed with respect to the classic and the Günter’s finite-difference schemes, both being based on non-aligned discretization but have been used in magnetic fusion modeling. Finally, a 3D test in a realistic tokamak geometry is studied in Section 8. Concluding remarks and perspectives are given in Section 9.