Advances in the flux-coordinate independent approach

doi:10.1016/j.cpc.2016.12.014

Computer Physics Communications

Volume 213, April 2017, Pages 111-121

https://doi.org/10.1016/j.cpc.2016.12.014 Get rights and content

Abstract

The flux-coordinate independent approach (FCI) offers a promising solution to deal with a separatrix and X-point(s) in diverted tokamaks. Whereas the discretisation of perpendicular operators (with respect to magnetic field) is straight forward, the major complexity lies in the discretisation of parallel operators, for which field line tracing and interpolation is employed. A discrete version for the parallel diffusion operator was proposed in Stegmeir et al. (2016), which maintains the self-adjointness property on the discrete level and exhibits only very low numerical perpendicular diffusion/pollution. However, in situations where the field line map is strongly distorted this scheme revealed its limits. Moreover, the appearance of small scale corrugations deteriorated the convergence order with respect to spatial resolution (Held et al., 2016). In this paper we present an extension to the scheme where the parallel gradient is reformulated via a combination of integration and interpolation. It is shown that the resultant scheme finally combines many good numerical properties, i.e. it is self-adjoint on the discrete level, it has very low numerical perpendicular diffusion, it can cope with strongly distorted maps and exhibits optimal convergence. Another subtle issue in the FCI approach is the treatment of boundary conditions, especially where magnetic field lines intersect with material plates. We present a solution based on ghost points, whose value can be set in a flexible way according to Taylor expansion around the boundary.

Introduction

The geometry of diverted magnetic fusion devices poses a challenge to the numerical treatment of the plasma edge and scrape-off layer (SOL). Dynamics in tokamaks is usually strongly anisotropic leading to structures which are strongly elongated along magnetic field lines $(k_{∥} ≪ k_{⊥})$ . This flute mode character is usually exploited computationally in numerical codes via employing field-aligned coordinate systems [1] and sparsifying the computational grid along the resulting parallel coordinate. However, field/flux-aligned coordinates become singular on the separatrix/X-point(s) and codes based on these coordinates have to provide some workaround or special treatment for the separatrix/X-point(s). Investigation of phenomena, where the X-point plays a crucial role, might therefore become questionable, and in the worst case even numerical artefacts might arise. The flux-coordinate independent approach (FCI) [2], [3], [4], [5] offers a solution to this dilemma: The simulation domain is spanned with a cylindrical grid, which is well defined everywhere in the region of interest, and the discretisation of perpendicular operators turns out to be straight forward. For the discretisation of parallel operators a field line map is used: Parallel operators are discretised via a finite difference along magnetic field lines, for which field line tracing towards neighbouring poloidal planes is performed and required values on the field line are obtained by interpolation. Finally, the flute mode character can be exploited computationally by sparsifying the grid along the toroidal direction. The FCI approach allows a high flexibility in geometry and is used in few codes, like FENICIA [3], GRILLIX [5], FELTOR [6], BOUT++ [7]. It has been successfully applied to hyperbolic and parabolic problems.

The main complexity of the approach lies in the discretisation of parallel operators, and a major critical issue is numerical perpendicular diffusion/pollution caused by the interpolation which couples distinct magnetic field lines. The highly anisotropic dynamics, e.g. the ratio of parallel to perpendicular heat conductivity may reach levels of $χ_{∥} / χ_{⊥} \sim 10^{10}$ , implies that even a small directional error of discrete parallel operators may overwhelm the real slow perpendicular dynamics. Following the method of support operators [8], [9] numerical schemes for the parallel diffusion operator were derived in [5], [10], which conserve the self-adjointness property on the discrete level and exhibit a highly reduced level of numerical diffusion as compared to a naive discretisation. Two types of schemes, one based on interpolation and another one based on integration were derived. However, it turned out that the schemes which are based on interpolation exhibit erroneous corrugations, especially in situations where the field line map is strongly distorted, i.e. at low toroidal resolutions and in presence of strong magnetic shear. Moreover, a deteriorated convergence behaviour was recently found from numerical tests in complicated geometries with $\nabla \cdot b \neq 0$ [6], with $b$ the magnetic field unit vector. In this paper we present an extension to these schemes, where the underlying parallel gradient is reformulated via a combination of interpolation and integration. We show that the resultant scheme combines many good numerical properties: It is self adjoint on the discrete level, it has very low numerical perpendicular diffusion/pollution, it can cope even with strongly distorted field line maps and it exhibits the expected second order convergence with respect to toroidal resolution. Moreover, it is also very practical as it does not require substantial additional effort in implementation and it does not increase the computational cost. This new numerical scheme constitutes a generalisation of and significant improvement to [5], [10], since many crucial problems are resolved without any drawback in its practical application.

Furthermore, we discuss the treatment of boundaries within the FCI approach. Especially at the intersection of magnetic field lines with material surfaces (limiter, divertor) subtle numerical problems may arise, as e.g. also information from outside the boundary is required for interpolation and the distances from grid points to the boundary may vary strongly. We present in this paper a solution which is based on ghost points, whose values are set according to a Taylor expansion along magnetic field lines around the boundary. This offers a simple and flexible way to deal with different kinds of boundary conditions.

All developed methods were implemented in the code GRILLIX, with which the presented numerical tests were carried out.

Section snippets

Parallel operators within FCI approach

A general introduction into the FCI approach, which is reviewed here only shortly, can be found e.g. in [2], [3], [4], [5] and we focus on the discretisation of parallel operators. We consider in the following an axisymmetric tokamak configuration which is spanned by a cylindrical grid $(R_{i}, Z_{j}, φ_{k})$ , which is Cartesian with grid spacing $h$ within poloidal planes $k$ . Based on the assumption of a strong toroidal field $(B_{t o r} ≫ B_{p o l})$ perpendicular operators can be approximated with a stencil which remains

Numerical tests

For the following tests and examples a 3rd order bipolynomial interpolation was used to obtain values at points which do not coincide with grid points (e.g. map points). Therefore, 4×4 grid points centred around the considered point are used to obtain the interpolating polynomial. The following numerical schemes were investigated:

•
D-3: Naive scheme for the parallel diffusion operator according to Eq. (12). The discrete parallel gradient is thereby computed by pure interpolation ( $X = 0$ or Eq. (1)

Boundaries within FCI

A subtle issue within the FCI approach is the treatment of boundaries, especially where magnetic field lines intersect with material plates, since these may have an immediate effect on the whole simulation domain. In general the intersection of field lines with boundaries can be identified during the field line tracing procedure.

An obvious problem arises if a map point is close to the boundary such that the standard interpolation stamp, (e.g. 4×4centred around the map point for third order

Conclusions

The flux-coordinate independent approach offers a viable solution to deal with complex geometries of magnetic fusion devices, especially a separatrix and X-point(s). The main challenge thereby is the discretisation of parallel operators as the numerical grid is not aligned with magnetic field lines.

We presented a new scheme for the parallel gradient which is based on a combination of interpolation and integration, where the degree of integration is controlled by a single parameter $X$ . With

Acknowledgements

The authors would like to thank L. Krönert for contributing to this work with useful discussions and fruitful comments. A part of this work was carried out using the HELIOS supercomputer system at Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC), Aomori, Japan, under the Broader Approach collaboration between Euratom and Japan, implemented by Fusion for Energy and JAEA. This work has been carried out within the framework of the EUROfusion Consortium and

References (14)

M. Ottaviani
Phys. Lett. A
(2011)
F. Hariri et al.
Comput. Phys. Comm.
(2013)
A. Stegmeir et al.
Comput. Phys. Commun.
(2016)
M. Held et al.
Comput. Phys. Commun.
(2016)
M. Shashkov et al.
J. Comput. Phys.
(1995)
W.D. D’haeseleer et al.
F. Hariri et al.
Phys. Plasmas
(2014)

There are more references available in the full text version of this article.

Cited by (14)

Numerical evaluation of line, surface and toroidal integrals on level sets of toroidally symmetric functions
2023, Journal of Computational Physics
We investigate strategies to numerically integrate closed lines and surfaces that are implicitly defined by level sets (iso-contours) of continuously differentiable toroidally symmetric functions. The “grid-transform” approach transforms quantities given in non-surface-aligned coordinates onto a numerically constructed surface-aligned grid. Here, line and surface integrals, as well as so-called flux-surface averages, can be easily evaluated using high order integration formulas. We compare this method to ones that base on numerical representations of the delta-function.
For the grid-transform method we observe high order convergence of line, surface and volume integration. Quantitatively, the errors for line and area integration are several orders of magnitude smaller than previously reported errors for delta-function methods. Furthermore, a delta-function method based on a Gaussian representation shows qualitatively wrong results of surface integrals near O- and X-points. Contrarily, the grid transform method suffers no deterioration near O-points. However, close to X-points we observe reduced first order convergence in volume integral and derivative tests due to the diverging volume element.
Finally, we derive a toroidal integration based on toroidal summation and a smoothing kernel that assumes field-alignment of structures between the toroidal planes. The smoothing kernel can be interpreted as a partial flux-surface average. The resulting smoothed toroidal average eliminates unphysical poloidal oscillations that are otherwise present in the simple toroidal average.
Our methods can be applied to toroidal and flux-surface averages in simulations of three-dimensional plasma dynamics on non-aligned grids. Further applications include closed line and surface integrals in level set methods. Efficient implementations can be found in the freely available Feltor library.
A finite volume flux coordinate independent approach
2023, Computer Physics Communications
We present a novel family of schemes as the merging between a one-dimensional advection scheme with the flux coordinate independent approach. The scheme can be used to discretize the field-aligned Navier-Stokes equations in three dimensions. Our approach consists of three major steps: (i) the formulation of the one-dimensional scheme in a locally field-aligned coordinate system, (ii) a numerical evaluation of the surface integrals over the field-aligned finite volumes and (iii) the introduction of smoothing into the numerical transformation operators to ensure stability of the resulting scheme.
We study this approach at the example of a staggered finite volume scheme with a field-aligned cylinder as initial condition. We show superior stability and conservative properties over previous direct discretizations. In particular, the relative mass conservation is improved by several orders of magnitude. Without smoothing in the transformation operators the scheme is prone to oscillations in both parallel and perpendicular directions. In the presence of strong perpendicular gradients, additional parallel diffusion is needed to control spurious oscillations in the perpendicular planes.
We provide parallel implementations for various platforms including GPUs freely in the C++ library Feltor
Analysis of locally-aligned and non-aligned discretisation schemes for reactor-scale tokamak edge turbulence simulations
2023, Computer Physics Communications
Edge turbulence codes will play a key role in the interpretation of ITER data and for reliable predictions of EU-DEMO. At present, such codes are not yet capable of routine simulations at reactor scale, and instead focus on smaller experiments like TCV or ASDEX Upgrade. Numerical methods have to be identified that scale to reactor size, and are able to cope with the highly anisotropic turbulent structures and simultaneously with the complex magnetic geometry in the edge with X-point(s), the vessel wall and divertor targets present. Particularly for such conditions two main approaches, non-aligned discretisation schemes and a locally-aligned discretisation scheme (commonly referred to as Flux-Coordinate Independent approach (FCI)), have emerged. We analyse both schemes concerning their applicability and scalability to next generation fusion reactors. We find that the ratio of the computational cost of a non-aligned scheme compared to aligned scheme scales as $\propto {(\frac{R_{0}}{q ρ_{s}})}^{2}$ for shear Alfvén dynamics, and as ${(\frac{R_{0}}{q ρ_{s}})}^{3}$ for electron heat conduction, where $R_{0}$ is the major radius of a tokamak, q an estimate for the safety factor at the edge and $ρ_{s}$ is the drift scale representing the typical size of turbulent structures to be resolved. Locally-aligned schemes can therefore be considered strongly favourable for reactor scale tokamaks concerning computational performance. On the other hand, locally-aligned schemes suffer from a more complex treatment of boundary conditions, for which an immersed boundary approach (IBA) was recently proposed. We demonstrate numerically the validity of this method in combination with the FCI approach. Finally, we present a first attempt of an ITER edge turbulence simulation with the FCI code GRILLIX at realistic parameters and in realistic geometry.
A partially mesh-free scheme for representing anisotropic spatial variations along field lines: Conservation, quadrature, and the delta property
2023, Computer Physics Communications
Citation Excerpt :
Since there is no defined connection between nodes in separate meshes, a mapping and interpolation procedure is performed in order to compute derivatives in the direction parallel to the magnetic field by mapping nodal points to their associated parallel ‘neighbour’ points on the adjacent meshes. Refinements to the original FCI method have been made both by Hariri and Ottaviani et al. [3,4] and Stegmeir et al. [5,6], and implementations of the method have been realized in several codes, including BOUT++, a plasma fluid simulation code [7,8]. This BOUT++ implementation is a finite difference scheme, and while it has shown some success at modelling plasma filaments [9] and arbitrary wall geometry [10], an even more general framework is desirable.
A method for representing highly anisotropic fields is presented, based on a partially meshfree Galerkin formulation. A mapping function is used to provide information about the local direction of the anisotropy, with one of the global coordinates chosen to parameterize the ‘parallel’ position along the mapping in a one-to-one manner. Standard unstructured finite element meshes are used on planes of constant parallel coordinate to represent the necessary small-scale variations perpendicular to the mapping direction, with large spacings then possible between these planes because of the large-scale variations along the mapping. This greatly reduces the number of degrees of freedom required to represent fields in this space and the associated computational cost of simulations involving such fields. No mesh connectivity is defined between planes, and field aligned basis functions are constructed using the mapping function to extend the standard finite element basis into the full domain. Integration of the basis has been addressed with reference to methods developed for fully meshfree methods, and the scheme (as well as other similar element-free Galerkin schemes) is shown to be locally conservative under certain conditions. Robust convergence of several test problems is demonstrated.
A new conservative finite-difference scheme for anisotropic elliptic problems in bounded domain
2020, Journal of Computational Physics
Highly anisotropic elliptic problems occur in many physical models that need to be solved numerically. A direction of dominant diffusion is thus introduced (called here parallel direction) along which the diffusion coefficient is several orders larger of magnitude than in the perpendicular one. In this case, finite-difference methods based on misaligned stencils are generally not designed to provide an optimal discretization, and may lead the perpendicular diffusion to be polluted by the numerical error in approximating the parallel diffusion.
This paper proposes an original scheme using non-aligned Cartesian grids and interpolations aligned along a parallel diffusion direction. Here, this direction is assumed to be supported by a divergence-free vector field which never vanishes and it is supposed to be stationary in time. Based on the Support Operator Method (SOM), the self-adjointness property of the parallel diffusion operator is maintained on the discrete level. Compared with existing methods, the present formulation further guarantees the conservativity of the fluxes in both parallel and perpendicular directions. In addition, when the flow intercepts a boundary in the parallel direction, an accurate discretization of the boundary condition is presented that avoids the uncertainties of extrapolated far ghost points classically used and ensures a better accuracy of the solution. Numerical tests based on manufactured solutions show the method is able to provide accurate and stable numerical approximations in both periodic and bounded domains with a drastically reduced number of degrees of freedom with respect to non-aligned approaches.
Tokamak H-mode edge-SOL global turbulence simulations with an electromagnetic, transcollisional drift-fluid model
2024, arXiv

View all citing articles on Scopus

View full text

Advances in the flux-coordinate independent approach

Abstract

Introduction

Section snippets

Parallel operators within FCI approach

Numerical tests

Boundaries within FCI

Conclusions

Acknowledgements

Phys. Lett. A

Comput. Phys. Comm.

Comput. Phys. Commun.

Comput. Phys. Commun.

J. Comput. Phys.

Phys. Plasmas