Discontinuous Galerkin methods for plasma physics in the scrape-off layer of tokamaks
Introduction
In magnetic fusion devices, a large fraction of the thermal power flows to an actively cooled strike-plate through a narrow layer called the scrape-off layer (SOL) [44]. This layer lies at the boundary between closed field lines and field lines that connect to the wall, and its width is determined by the competition between the turbulent transport across the magnetic field and the very rapid transport along the field. This width determines the heat flux density on the wall and impacts wall erosion, recycling, and density control. Transport in the SOL is known to be highly intermittent and dominated by the ejection of coherent plasma filaments known as “blobs” or “streamers” [22], [35], [64]. The modeling of this turbulent transport is important in order to predict and understand the dependence of the SOL width on the plasma parameters and machine size.
The equations that govern SOL blob dynamics have been formulated and analyzed [5], [26], [35], [62], but many questions remain open. From a mathematical point of view the dynamics of existing SOL plasma models can be characterized as systems of equations that are driven by operators that are nonlinear and contain multiple signatures (e.g. elliptic, parabolic, dispersive, hyperbolic operators, etc.) displaying weak regularity features in the coupling. As the consistency, stability, and accuracy of numeric methods are strongly constrained by the behavior of the underlying mathematical models, the difficulties posed by these complicated SOL plasma dynamics raise interesting and challenging technical questions.
As a potential solution to the difficulties that can arise, we explore the use of a discontinuous Galerkin (DG) finite element method for modeling blob dynamics in the SOL. The discontinuous Galerkin method is a high order numerical method, that has been found to provide formidable accuracy, stability, and robustness in many areas of nonlinear dynamics [15], [16], [25], [37]. Moreover, DG methods are unified in the sense that they are well-suited for rigorous analysis of both the physics as well as the numerics and mathematics of the system. They demonstrate high order convergence rates [2], [4], are well-established as candidates for computationally optimal adaptive technologies (e.g. hp-adaptivity and r-adaptivity) [21], [52], are extremely scalable (especially on state-of-the-art architectures with high thread parallel arithmetic intensity, such as GPUs and coprocessors) [1], [34], [73], [75], and are often noted as being remarkably flexible for accurately modeling large categories of coupled systems of initial-boundary value problems with strongly nonlinear forcings. This aspect of DG methods makes them particularly appealing for studying scrape-off layer dynamics, as the system of PDEs in question in the scrape-off layer is often highly variable, requires great flexibility in representation, possess weak regularity features, and may involve complicated geometries (e.g. in the presence of magnetic chaos) with nonlinear boundary forcings that are befitting to finite element methods in general.
The application of DG methods to problems arising in plasma physics and nuclear engineering have established a high benchmark. Dawson, Wheeler and Proft studied neutron transport in [20]. In [31], [32] Hesthaven and Jacobs addressed DG plasma building block models, along with providing high resolution insight into the basic aspects of applying DG methods to plasma problems. Warburton and Karniadakis developed a turbulent DNS solver for unsteady viscous MHD in [72], and Cheng, Li, Qiu, and Xu [14] have recently studied specialized positivity preserving DG methods in the context of ideal MHD. In plasma kinetic theory, Heath, Gamba, Morrison, Cheng, and Michler [13], [28] have developed some impressive high dimensional Vlasov–Poisson schemes, where Rossmanith and Seal have similarly implemented DG schemes for Vlasov–Poisson, though in the Semi-Lagrangian setting [61]. Edge turbulence models have recently worked their way into some studies of Peterson and Hammet [59], while some of the most well-established DG results in computational plasma physics are those from the multifluid models of Shumlak, Loverich, Hakim, Srinivasan, and Meier et al. [45], [46], [48], [67], [68], which have not only been thoroughly validated, but thoroughly benchmarked against established codes in the community.
The present paper represents an interdisciplinary effort to apply general discontinuous Galerkin methods to the modeling of plasma dynamics in tokamak scrape-off layers. This work draws heavily from a number of disciplines, in the context of trying to serve as a powerful staging ground for further and deeper analysis of tokamak reactors. In this direction, the code architecture we have developed for this study has been given the moniker ArcOn. The code itself is explicit in time and consists of three coupled equations that require two local discontinuous solves, and a global linear solve in the elliptic subsystem (to be described below). It has been written (primarily) in C++, and utilizes a number of external libraries, including: deal.II [7], p4est [12], and PETSc [6]. The development of ArcOn has required the modification of deal.II, which required a fundamental update/extension to core functionality in order to provide inter-software support for periodic boundary data over massively parallel distributed meshes. We mention this process briefly here, as it was quite time-consuming and introduced a nontrivial technical challenge in the computer science aspect of the code development. These modifications were performed in consultation with project leaders, and have subsequently been openly distributed to their respective communities for broad use.
All computations below were run in parallel on the Texas Advanced Computing Center's (TACC's) 10 PetaFLOPS Dell Linux Cluster based Stampede system, comprised of 6400+ Dell PowerEdge server nodes, each outfitted with 2 Intel Xeon E5 (Sandy Bridge) processors and an Intel Xeon Phi Coprocessor (MIC Architecture). For this study, relatively small meshes were used at low polynomial order, that ranged from tens of thousands of degrees of freedom, to tens of millions. All runs were performed in parallel on between 16 to 256 cores.
An outline of the paper is as follows. In Section 2 we present the physics-based SOL model for filamentary blob transport. Here we provide motivation for the system formulation, and heuristically discuss its features, scope, and limitations. In Section 3 we review the formulation of the discontinuous Galerkin method applied to the mathematical model presented in Section 2. We address some salient features of ArcOn, including the ability to use both multiple solvers (e.g. SIPG, NIPG, etc.) and various basis functions, including both nodal and modal finite elements. We also describe the upwinding flux schemes used in the approximate Riemann-type solvers, and the strong stability preserving temporal discretization scheme employed. In Section 4 we discuss the parameters of a classical verification study by way of a manufactured solution, and perform numerical experiments that demonstrate the optimal (to super-optimal in places) convergence order in each operator subsystem of ArcOn. In Section 5 we discuss some of the basic stability features of the DG method, and how regularization procedures can be used to ensure numerical stability. We present here a basic blob example adjoined to a new artificial diffusion scheme, and run using a CFL-based variable timestepping algorithm. Finally in Section 6 we show a modon solution to the physics of the system, briefly analyze and validate some of the cogent numerical properties of blobs, and demonstrate the effectiveness of ArcOn in solving turbulent plasma dynamics at saturation.
Section snippets
The governing equations
We use a two-dimensional model for the SOL turbulence that describes the evolution of the density n, the vorticity U, and the electrostatic potential ϕ. The latter plays the role of a stream-function for the ion velocity, which is given by the perpendicular drift velocity (here and is the electric field). The direction is chosen to lie along the magnetic field, . The model is based on the assumption that due to rapid transport along the magnetic field, the
The discretization procedure
Let us discretize our spatial domain Ω. Consider the open set with physical boundary ∂Ω, given such that . Let denote the partition of the closure of the polygonal mesh of Ω, which we denote , into a finite number of polygonal elements denoted , such that , for the number of elements in . Here and below the mesh diameter h is chosen to satisfy for the distance function and elementwise face vertices when the
Numerical tests
The numerical residual for any solution component (e.g. we use ϕ here) is given by , where ϕ is the solution to (2.3) and its discrete approximate counterpart. This measure will be enough to determine the convergence order and many aspects of the numerical behavior of the solution relative to its analytic counterpart. Notice that since ArcOn utilizes both modal and nodal basis functions, there is a potential ambiguity that arises in the error measures. Namely, the initial state
Regularization processes
As the system of equations (2.1), (2.2), (2.3) is convection-dominated, dispersive, gradient driven and nonlocal, it presents its own set of numerical complications and nuances. For example, the cross product terms induce convective nonlinearities that are driven by dispersive modes. As a result, in a strict mathematical sense the natural function space that support solutions to such a system is highly irregular, and in all likelihood local in its signature behavior (e.g. ). While this
Modon solutions
Physical modons (or dipole drift vortices) emerge as important solutions in a number of fields, such as geophysical fluid dynamics [42], coastal dynamics [65] and atmospherics [27], [57]. Intuitively these objects are often discussed in tandem with their counterpart, monopolar vortices, which in the parlance of atmospheric science correspond to phenomena such as cyclones. In plasma physics applications, monopolar and dipolar drift vortices (modons) [49] are often described as the coherent
Conclusion
We have presented a new architecture called ArcOn for studying the dynamics of filamentary blobs transported through the scrape-off layer of tokamaks. In this regard we have implemented a fully discontinuous Galerkin method for solving (2.1), (2.2). In contrast to mixed form finite element methods, for example, which use continuous Galerkin projections [53] to recover solutions to Poisson, our present approach preserves the discontinuous function spaces throughout the computation, thus
Acknowledgements
The first author would like to thank Phillip Schmitz, John Evans, Jed Brown, Timo Heister, Ammar Hakim, Paul Bauman, Georg Stadler, Logan Moon, Kyle Mandli, Andy Terrel, Clint Dawson, and Irene Gamba for fun, insightful, and incredibly helpful discussions. This work was supported by the U.S. Department of Energy under Contract No. DE-FG02-04ER-54742.
References (76)
- et al.
Discontinuous Galerkin for high performance computational fluid dynamics (hpcdg)
- et al.
Discontinuous Galerkin methods for elliptic problems
An interior penalty finite element method with discontinuous elements
SIAM J. Numer. Anal.
(1982)- et al.
Unified analysis of discontinuous Galerkin methods for elliptic problems
SIAM J. Numer. Anal.
(2001/02) Convective transport in the scrape-off layer of tokamaks
Phys. Plasmas
(Jun. 2005)- et al.
PETSc users manual
(2013) - et al.
deal.II differential equations analysis library, technical reference
- et al.
A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations
J. Comput. Phys.
(1997) - et al.
A discontinuous hp finite element method for convection–diffusion problems
Comput. Methods Appl. Mech. Eng.
(Jul. 2, 1999) - et al.
A bounded artificial viscosity large eddy simulation model
SIAM J. Numer. Anal.
(2009)