# Gyrokinetic Continuum Simulation of Turbulence in Open-Field-Line Plasmas

## Abstract

The properties of the boundary plasma in a tokamak are now recognized to play a key role in determining the achievable fusion power and the lifetimes of plasma-facing components. Accurate quantitative modeling and improved qualitative understanding of the boundary plasma ultimately require five-dimensional gyrokinetic turbulence simulations, which have been successful in predicting turbulence and transport in the core. Gyrokinetic codes for the boundary plasma must be able to handle large-amplitude fluctuations, electromagnetic effects, open and closed magnetic field lines, magnetic X-points, and the dynamics of impurities and neutrals. The additional challenges of boundary-plasma simulation necessitate the development of new gyrokinetic codes or major modifications to existing core gyrokinetic codes.

In this thesis, we develop the first gyrokinetic continuum code capable of simulating plasma turbulence on open magnetic field lines, which is a key feature of a tokamak scrape-off layer. In contrast to prior attempts at this problem, we use an energy-conserving discontinuous Galerkin discretization in space. To model the interaction between the plasma and the wall, we design conducting-sheath boundary conditions that permit local currents into and out of the wall. We start by designing spatially one-dimensional kinetic models of parallel SOL dynamics and solve these systems using novel continuum algorithms. By generalizing these algorithms to higher dimensions and adding a model for collisions, we present results from the first gyrokinetic continuum simulations of turbulence on two types of open-field-line systems. The first simulation features uniform and straight field lines, such as found in some linear plasma devices. The second simulation is of a hypothetical model we developed of the NSTX scrape-off layer featuring helical field lines. In the context of discontinuous Galerkin methods, we also explore the use of exponentially weighted polynomials for a more efficient velocity-space discretization of the distribution function when compared to standard polynomials. We show that standard implementations do not conserve any moments of the distribution function, and we develop a modified algorithm that does. These developments comprise a major step towards a gyrokinetic continuum code for quantitative predictions of turbulence and transport in the boundary plasma of magnetic fusion devices.

September 2017 \copyrightyear2017
\adviserGregory W. Hammett \departmentAstrophysical Sciences

Program in Plasma Physics

###### Acknowledgements.

Many fall in the face of chaos, but not this one… not today.

Narrator, Darkest Dungeon

After more than six years at Princeton, I am finally ready to move on. Although I am eagerly anticipating what the future holds, I will also miss the tightly structured life that I have become accustomed to over the years. I would not have been able to complete my thesis without the direct support of many people. Of course, I must thank Greg Hammett for being my adviser. Although I did not decide to come to Princeton to work with him (or anyone else) in particular, I realized later on in my graduate career how lucky I was to find an adviser who was an expert in both turbulence and numerical methods, had a code-development project written in a modern programming language, and had the patience to teach a student who knew so little. It is fitting that we both have a tendency to send e-mails well-past midnight. I will be happy if I can one day attain even a small fraction of his level of knowledge, creativity, and intuition. I thank Tim Stoltzfus-Dueck for always making the time to share his incredible knowledge and powerful intuition about plasma physics with me. Tim has been a great mentor and teacher, often meeting with me several times a week, and I owe much of what I learned about SOL physics to him. Tim played an essential role in interpreting the physics of the simulations and suggesting systematic ways to test the code when the results appeared to be incorrect. Tim’s commitment to do physics in an objective and transparent way has greatly influenced the way I conduct and present my own work. I hope that Tim will one day use his remarkable aptitude for teaching to train future generations of plasma physicists in edge and SOL physics. I am very grateful to Ammar Hakim for creating the Gkeyll plasma-physics framework and for developing some of the additional capabilities required for my thesis. I owe a large part of the timely completion of my thesis (by Princeton standards, at least) to Ammar, who has designed the code with modern and sound software-engineering concepts. I would also like to thank Ammar for teaching me the basics about DG methods and for closely working with me for my second-year project. I thank John Krommes for being an excellent teacher, academic adviser, and even boss. Although he must have been disappointed in how little I learned from his ‘Irreversible Processes in Plasmas’ class when I took it as a second-year student, I hope he recognized that I finally started learning something in the two semesters during which I worked as a grader for his class. As my academic adviser, John always listened intently to whatever dilemmas I brought to him and drew upon his experience to provide profoundly wise advice. I also credit John with instilling a strong belief in adhering to established rules and conventions for technical writing. Unfortunately, I will not return your $31 if you catch five writing mistakes in my thesis. I thank my thesis readers, Stewart Zweben and Matt Kunz, for carefully reading my thesis and providing valuable feedback. Greg Hammett and Tim Stoltzfus-Dueck also provided valuable comments for my thesis. I would also like to thank Stewart Zweben (figures 2 and 3), the IAEA and Thomas Eich (figure 4), and Troy Carter and David Schaffner (figures 25, 26, and 38) for permission to use their figures in my thesis to illustrate important results. I thank my other research collaborators for taking the time to learn about my research or code: Qingjiang Pan, Tess Bernard, and Jimmy Juno. I credit Noah Mandell with the formidible task of building the Gkeyll code on various systems. I also thank Dan Ruiz and Jeff Parker for enlisting my help on a side project about DW–ZF interactions. I am very grateful to Curt Hillegas and Matt Kunz for their support for my use of the Perseus cluster at Princeton, which enormously sped up my time to graduate. I also thank Frank Jenko and Jason TenBarge for providing access to additional computational resources. I thank Troy Carter, Amitava Bhattacharjee, Paolo Ricci, Stewart Zweben, Frank Jenko, Bill Tang, Hantao Ji, and Bob Kaita for their encouragement and support. I also acknowledge useful discussions with John Krommes, Troy Carter, Paolo Ricci, Stewart Zweben, Denis St-Onge, Seung-Ho Ku, and Michael Churchill. I thank Dan Ruiz for his consistent friendship over the year and being a great office mate. He probably would have graduated even sooner if he did not have me as an office mate to provide constant distractions. I am in awe of his personal integrity and ability to conjure original research ideas. I am also grateful to the friends I made at Princeton: Mike Hay, Jeff Parker, Seth Davidovits, Lee Ellison, Jake Nichols, Yao Zhou, Kenan Qu, Chang Liu, and Yuan Shi. I thank Beth Leman, Barbara Sarfaty, and Jen Jones for shielding me from the host of complicated administrative details associated with being a graduate student at PPPL. Lastly, I thank my family for their unconditional support. Rita, I would not have been so determined in my goals if I did not have an older sibling to show me what could be achieved. Mom and Dad, thank you for all the sacrifices you have made to help your children pursue their dreams. You have given me so much and have never asked for a single thing in return. This work was funded by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466, through the Max-Planck/Princeton Center for Plasma Physics and the Princeton Plasma Physics Laboratory. Some simulations in Chapters \thechapter and \thechapter were performed on the Perseus cluster at the TIGRESS high-performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology’s Research Computing department. Initial development on the simulations in Chapter \thechapter used the Edison system at the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The simulations in Chapters \thechapter and \thechapter also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. \dedicationTo my parents.

## Chapter \thechapter Introduction

### 1 The Boundary Plasma

On a basic level, the plasma in a tokamak can be separated into core, edge,
and scrape-off-layer (SOL) regions.^{1}

While a significant fraction of the tokamak plasma lies in the core and edge, there exists a region outside the edge called the SOL where the magnetic field lines no longer trace out closed flux surfaces and instead intersect material walls after winding around toroidally a number of times (in practice it is not possible to create a magnetic field with field lines that are perfectly tangential to the wall at every point (Stangeby, 2000; Ricci, 2015)). Charged particles in the SOL are therefore rapidly lost to material surfaces in contact with the magnetic field lines, where recombination occurs. On a basic level, the SOL properties are set through a balance between plasma outflow from the edge, cross-field turbulent transport, and strong parallel losses at the divertor or limiter plates (Ricci et al., 2015; Mosetto, 2014), where a Debye sheath layer forms to keep electron and ion particle fluxes to the wall approximately equal. Plasma–surface interactions (PSIs) such as recycling and impurity influx are also important in setting the particle and power balances in the SOL. Due to plasma–wall interactions, the SOL plasma is much colder (– eV (Zweben et al., 2007; Stangeby, 2000; Stoltzfus-Dueck, 2009)) than the core and edge plasmas. The SOL and edge are separated by a boundary referred to as the last-closed flux surface (LCFS).

The properties of the boundary plasma constrain the performance and component lifetime of tokamak fusion reactor by affecting both the details of how heat is exhausted in the SOL as well as how much fusion power can be generated in the core, assuming that core profiles are stiff (Kotschenreuther et al., 1995; Doyle et al., 2007; Kardaun et al., 2008),. A thorough understanding of the physics of boundary plasma is therefore key to improving the overall viability of the tokamak concept, but there are still many gaps in our present understanding that need to be filled in before ITER begins the very long-awaited D–T operations in . Given the complexities of edge and SOL turbulence, numerical simulations have become important in furthering our theoretical understanding of the boundary plasma. This thesis attempts to add to this line of research by developing gyrokinetic simulations of turbulence in open-field-line (as in the SOL) plasmas using a class of grid-based (‘continuum’ or ‘Eulerian’) numerical methods. In this introduction, we briefly review some basic features of the SOL and motivate the usefulness of gyrokinetic simulations as a tool to address important boundary-plasma physics questions in present-day and future tokamaks.

#### 1.1 Basic Features of the SOL

The SOL plasma is defined as the plasma in the region extending from the LCFS to the material wall. Since the SOL plasma is in direct contact with solid surfaces, PSIs that lower plasma temperatures through radiative processes are inevitable (Stangeby, 2000). Just outside the LCFS, the SOL features steep exponential profiles and near-Gaussian probability distribution functions (PDFs). As one goes farther out radially in the SOL from the LCFS, the profiles become slowly decreasing or flat, and the PDFs become increasingly non-Gaussian (Zweben et al., 2007). Relative electron density fluctuation levels in the SOL are also fairly large, increasing from near the LCFS to near the first wall (Zweben et al., 2007). Figure 2 shows typical profiles of the D light emission measured in the boundary plasma on NSTX using a gas puff imaging (GPI) diagnostic (Zweben et al., 2017). The D light emission is influenced by the neutral D gas-puff density, the local electron density, and the local electron temperature (Zweben et al., 2017), but the relative fluctuation levels should mostly be set by plasma fluctuations.

On the open magnetic field lines in the SOL, charged particles rapidly flow along the field lines towards solid surfaces (e.g., walls, limiters, divertors), where they recombine and are lost (Chen, 1984). Since electrons are much more mobile than ions, they initially are lost to the surfaces in contact with the magnetic field lines faster than the ions, leaving the plasma with a thin layer of net positive charge, which is confined to a layer at the plasma–material boundary a few Debye lengths wide due to Debye shielding. This layer is called the electrostatic Debye sheath (Stangeby, 2000; Chen, 1984) and its primary purpose is to establish a potential barrier that repels incident electrons and accelerates incident ions into the surface, keeping the particle fluxes of electrons and ions lost to the surfaces approximately equal. The plasma itself maintains quasineutrality. Only electrons with sufficient parallel velocity can surmount the sheath potential drop to reach the surface. The sheath plays an important role in governing both the properties of the SOL plasma and how particles and energy are lost to solid surfaces.

As charged particles in the SOL rapidly flow along field lines to solid surfaces through parallel motion, they also undergo a much slower but non-negligible motion across the magnetic field. As in the core, this cross-field transport is dominated by turbulence, but unlike the core, a significant fraction of the cross-field transport in the SOL (Boedo, 2009; Boedo et al., 2003) is due to the radial convection of coherent structures called blobs or filaments (Zweben et al., 2004; Terry et al., 2007; Boedo et al., 2014; Zweben, 1985; Zweben & Gould, 1985). Plasma blobs are highly elongated along the field line, with typical parallel scales 1–10 m due to the rapid parallel motion of charged particles and much smaller cross-field scales 1–10 cm (Zweben et al., 2017). Blobs have elevated densities and/or temperatures compared to the background plasma. Although the blob-formation mechanisms are not as well-understood as their transport mechanisms (Krasheninnikov, 2016), they are commonly observed to form in the edge, from where they are ejected across the LCFS into the SOL (Boedo et al., 2003; Terry et al., 2003; Zweben et al., 2004). Cross-field transport in the SOL is highly intermittent due to blob propagation (Zweben et al., 2007) and is consequently poorly described in terms of effective diffusion coefficients and convective velocities (Naulin, 2007).

Figure 3 shows raw camera images from a GPI diagnostic viewed along the local magnetic field direction near the outer midplane of the NSTX tokamak. These images show the 2D cross-field blob structure of a 3D plasma filament as it is ejected across the LCFS into the SOL.

The standard picture of blob propagation involves the polarization of a blob by the curvature drift, which self-generates a vertical electric field across the blob that results in a radially outward drift (Krasheninnikov, 2001; Krasheninnikov et al., 2008; Grulke et al., 2006). The efficient and fast convective radial transport of blobs towards the main chamber walls (instead of along the field lines to the divertor) can result in damage to first-wall components and the contamination of the core plasma by wall impurities (Rudakov et al., 2005; Terry et al., 2007), which lowers the plasma temperature through radiative cooling. Both the generation of blobs in the edge and the structure and motion of blobs in the SOL (Zweben et al., 2016, 2006) are ongoing research topics of practical interest.

**
Particle and Heat Loss in the SOL** The loss of particles and heat in the SOL to plasma-facing components (PFCs) is a major concern for
future high-power devices like ITER (Lipschultz et al., 2012) and Demo because a significant fraction (20%) of the
heat produced in the core is transported across the LCFS into the SOL by convection and conduction
(Loarte et al., 2007), where the heat must be exhausted somehow.
The unmitigated steady-state parallel heat flux in the SOL is expected to
be GW m for ITER (Loarte et al., 2007) and GW m for Demo (Goldston, 2015),
while material limitations set the maximum tolerable heat flux normal to the divertor plates at MW m
for steady-state and MW m for transients (Loarte et al., 2007), such as from disruptions
and edge-localized modes.
An often-quoted comparison in the fusion community is the transient heat flux of approximately 6 MW m
experienced by some components of a space vehicle during atmospheric reentry.
The use of an extremely shallow incidence angle between the divertor plates and field lines (–)
significantly reduces the heat flux normal to the divertor plates,
but up to of the power (Goldston, 2015) may need to be dissipated in the SOL through various means
before reaching the divertor plates to bring normal heat fluxes down to tolerable levels.
Strategies to reduce the heat load on the divertor plates are still under development
and include radiative divertor detachment (Soukhanovskii, 2017), advanced divertor geometries (Kotschenreuther et al., 2016; Ryutov et al., 2008; Valanju et al., 2009),
and applied resonant magnetic perturbations (Ahn et al., 2010; Evans et al., 2005; Jakubowski et al., 2009).

The heat fluxes on the ITER divertor targets are believed to be enormous because experiments on present-day machines indicate that the heat exhausted into the SOL flows towards the divertor plates in an extremely narrow channel whose width at the outer midplane (quantified as an exponential decay length) is insensitive to the machine size (Loarte et al., 2007; Makowski et al., 2012; Eich et al., 2013). There exists major uncertainty about the validity of empirical extrapolations to ITER, however. The amount of power spreading along the ITER divertor legs is also not well understood (empirically or theoretically), and it is possible that power-spreading effects on ITER will be principal in setting the heat-flux width at the divertor plates, making such considerations of the midplane heat-flux width unimportant for predicting divertor-plate heat loads. Presumably, the location and strength of heat loads deposited on PFCs is set through a balance between confined plasma outflow across the LCFS, parallel losses at the sheaths, and cross-field turbulent or neoclassical transport in the SOL. Therefore, credible numerical investigations of the SOL heat-flux width require the use of sophisticated turbulence codes.

Figure 4 shows a plot from Eich et al. (2013) of the outer-midplane SOL heat-flux width versus the outer-midplane poloidal magnetic field computed from a multi-machine database of low-recycling H-mode tokamak discharges. In that study, the outer-midplane heat-flux width for each discharge was determined as a fit parameter to match infrared thermography measurements of the heat-flux profile at the outer divertor target (Eich et al., 2011). Eich et al. (2013) found that “the strongest and essentially only dependence amongst the regression variables tested, at least for the conventional aspect ratio tokamaks, is an inverse scaling with plasma current (or equivalently [an inverse] dependence on outboard midplane poloidal magnetic field)”. Eich et al. (2013) and Makowski et al. (2012) also extrapolated the scaling to the ITER H-mode diverted plasma, finding an outer-midplane SOL heat-flux width of mm, much smaller than the 5 mm value used in ITER design specifications and the 3.6 mm value from transport modeling (Kukushkin et al., 2013). Goldston (2012) developed a heuristic drift-based (HD) model of the SOL heat-flux width that was consistent with the results of Eich et al. (2013) and Makowski et al. (2012). The HD model notably does not attempt to include turbulence, which could spoil its extrapolation to new parameter regimes (such as the ITER SOL).

It is important to study the heat-flux-width problem using turbulence simulations to provide a first-principles-based check of the empirical predictions and to investigate ways to broaden the heat-flux width, such as by increasing turbulent heat transport along the divertor legs. A recent electrostatic gyrokinetic simulation predicted the outer-midplane heat-flux width on ITER to be approximately 5.6 mm (Chang et al., 2017), but no other gyrokinetic codes currently have the capability to cross-check this result. More research is required to understand the physics setting the SOL heat-flux width (at both the outer midplane and at the divertor plates), the validity of such empirical extrapolations to ITER, and the implications of such a narrow width on divertor-plate heat loads (Goldston, 2015).

#### 1.2 Basic Plasma Physics Experiments

While experimental measurements of the SOL plasma in tokamaks have been highly successful in improving our understanding of SOL turbulence (Zweben et al., 2007; Tynan et al., 2009), the complicated physics situation in a tokamak can make detailed comparisons of analytical theories and numerical simulations with experimental data extremely challenging. Fortunately, there are a number of basic plasma physics experiments that produce highly reproducible, well-diagnosed plasmas in simple, open-magnetic-field-line configurations. Measurements of turbulence and transport in these laboratory plasma devices can shed light on topics of relevance to tokamak plasmas, such as sheared-flow suppression of turbulence (Carter & Maggs, 2009; Schaffner et al., 2013; Gentle et al., 2010), turbulence intermittency (Carter, 2006; Windisch et al., 2011), and turbulence saturation. The plasmas are usually relatively cold ( eV) and low pressure (), so electrostatic 3D fluid simulations are commonly used to study the turbulence in basic plasma physics experiments (Friedman et al., 2013; Ricci & Rogers, 2009; Rogers & Ricci, 2010; Naulin et al., 2008). Although the plasmas created in these devices are not in the same parameter regime as fusion plasmas, the free energy sources, dissipation mechanisms, and nonlinear transfer processes identified in these devices may be generic and thus relevant to tokamaks (Tynan et al., 2009; Fasoli et al., 2006).

Linear devices create a column of plasma in a uniform, axial background magnetic field. Some examples of linear devices include LAPD (Gekelman et al., 1991, 2016), CSDX (Burin et al., 2005; Thakur et al., 2014), VINETA (Franck et al., 2002), HelCat (Gilmore et al., 2015), and Mirabelle (Pierre et al., 1987; Brandt et al., 2011). Turbulence in linear devices is typically driven by temperature or density profile gradients (drift-wave instabilities). Toroidal basic plasma physics devices, called simple magnetized tori (SMTs), create a magnetic geometry composed of helical magnetic field lines (a superposition of vertical and toroidal magnetic fields) that intersect the vessel on each end. Compared to linear devices, the turbulence and transport in SMTs are more relevant to a tokamak SOL because they can also be driven by magnetic-field-line-curvature effects, which are important in the generation of blobs (as we will see in Chapter \thechapter). Examples of SMTs include TORPEX (Fasoli et al., 2006) and Helimak (Gentle & He, 2008).

The relative simplicity, availability of detailed diagnostics, and ease of parameter scans in basic plasma physics experiments also facilitates comparisons between simulations and experiments. The Global Braginskii Code (GBS) has been used to perform rigorous validation studies with data from TORPEX (Ricci et al., 2015, 2011, 2009). In the development of complex numerical codes, comparisons with basic plasma physics experiments can be made relatively early on with a minimal set of features implemented in the code to benchmark and test the models and numerical algorithms in a physically relevant setting. We adopt this strategy in this thesis, where we first simulate turbulence in LAPD to test 5D gyrokinetic continuum simulations before adding additional levels of complexity.

### 2 Boundary-Plasma Modeling

Due to the complications of the edge and SOL, numerical simulations have become an important tool for elucidating the physics of the boundary plasma. The use of a particular simulation can be classified as interpretive, predictive, or generic (Mitchell et al., 2000). Since these classifications are seldom explained in plasma-physics literature, and sometimes incorrectly applied, we first provide a brief description of the different types of modeling.

Interpretive modeling. To aid in experimental analysis and physical understanding, interpretive modeling uses simulations with data collected from an experiment as inputs to estimate other variables described by reduced models. These estimated variables are often difficult or impossible to measure in the actual experiment due to diagnostic limitations. Some interpretive models have several free parameters that can be tuned automatically so that a subset of the output data matches the corresponding experimentally measured data (e.g., density and temperature profiles), making it straightforward (but somewhat questionable) to run additional simulations in a predictive manner by using hypothetical initial conditions and plasma parameters (Ricci, 2015). These tools have become standard in analyzing data from tokamak discharges to infer transport levels and energy confinement. Since interpretive simulations are limited by the fidelity of the underlying reduced models (which can be empirical rather than physics based) and by the availability and quality of experimental data, caution must be used when drawing physics conclusions based on interpretive simulations and when extrapolating the models to parameter regimes that are outside the range of demonstrated validity.

Predictive modeling. As their name suggests, the goal of predictive models is to make quantitative predictions about the physics given a collection of initial conditions and assumptions (e.g., plasma sources, boundary conditions, fluid vs. kinetic description). They are particularly valuable for testing and improving interpretive models and predicting the plasma behavior in situations that cannot be investigated by experiments (such as in future fusion reactors like ITER). The use of free parameters in predictive models should be minimized, and the insensitivity of results to specific values of the model free parameters must be verified. Ideally, predictive simulations should solve first-principles equations without ad hoc terms and empirical models (Barnes, 2008; Candy et al., 2009; Görler et al., 2014). However, such first-principles simulations for a given application (e.g., turbulence on transport time scales) can be computationally infeasible even on state-of-the-art high-performance computers, in which case reduced models must be used in some capacity (Kotschenreuther et al., 1995; Staebler et al., 2007; Kinsey et al., 2011). A code that can make quantitatively accurate predictions of the boundary plasma is likely to require a high degree of sophistication, having models of the numerous PSIs that occur in the SOL. An example of a boundary-plasma predictive simulation is the use of the XGC1 gyrokinetic PIC code to calculate the SOL heat-flux width in ITER (Chang et al., 2017), although this is a prediction that cannot be validated for many years. Predictive simulations for high-fusion-power ITER regimes are particularly valuable because mistakes in the design of ITER based on our understanding of existing experiments can result in irreparable damage to the machine and to the overall ITER project.

Generic modeling. Generic modeling concerns the modeling of hypothetical situations, often with a number of highly simplifying assumptions. This type of modeling is useful for helping assemble a theoretical picture to explain observed phenomena and for testing the theoretical models of others using first-principles simulations. Since the purpose is neither to analyze a particular experiment nor to make definite predictions, sophisticated physics models that are essential for accurate predictive modeling can be neglected, which reduces the computing requirements and code complexity. Some uses of generic modeling for the boundary plasma include understanding the low-to-high confinement-mode (L–H) transition in tokamaks (Connor & Wilson, 2000), which is a rapid bifurcation between two plasma states associated with the formation of an edge transport barrier. Many studies of the mechanisms driving turbulence in the edge and SOL also employ generic modeling (Scott, 2005a; Ribeiro & Scott, 2008; Halpern et al., 2016). The simulations that are presented in this thesis also fall into this category.

#### 2.1 Fluid Modeling

Fluid-based codes have been and still are commonly used to model the boundary plasma. One class of these codes solves simplified transport equations based on the Braginskii fluid equations in two dimensions assuming axisymmetry. Plasma turbulence is not captured in these models, so the use of ad hoc anomalous diffusion terms (Dekeyser et al., 2011) or coupling to a turbulence code (Schneider et al., 2006; Rognlien et al., 2004) is required to model the turbulent transport across the magnetic field. Perhaps the most widely used of these fluid plasma transport codes is the SOLPS package (Schneider et al., 1992; Rozhansky et al., 2009), which contains a 2D boundary-plasma transport component coupled to a kinetic Monte Carlo model for neutral transport (Reiter et al., 2005). SOLPS has become the standard tool used for ITER divertor and boundary-plasma modeling (Wiesen et al., 2015). Some physics models implemented in SOLPS include impurities, charge-exchange, ionization, radiation, and sputtering. Other boundary-plasma transport codes which also solve similar equations include UEDGE (Rognlien et al., 1999, 1994) and EDGE2D (Simonini et al., 1994). Although interpretive transport codes have become packed with features over decades of development, incorporating many sophisticated models for physics believed to be relevant in the boundary plasma (Schneider & Runov, 2007), the basic fact that they often do not capture plasma turbulence self-consistently makes it questionable to use them in a predictive capacity, especially as a principal tool for quantitative modeling of the challenging SOL in ITER, a machine which costs billion just to construct (Kramer, 2016). Some shortcomings of fluid plasma transport codes for boundary-plasma modeling are discussed by Boedo et al. (2009).

Boundary-plasma turbulence codes solving drift-ordered fluid equations
(Zeiler et al., 1997; Simakov & Catto, 2003; Scott, 1997; Xu & Cohen, 1998) in two or three dimensions
have also been developed (Ricci et al., 2008; Dudson et al., 2009; Naulin et al., 2008; Xu et al., 2000, 2008; Tamain et al., 2010),
although many codes employ approximations or omit physics that might have significant
effects on the results, such as neglecting adiabatic coupling^{2}^{3}

There are also boundary-plasma turbulence codes that solve 3D electromagnetic gyrofluid equations (Ribeiro & Scott, 2005, 2008; Xu et al., 2013; Kendl et al., 2010), which are more robust than drift-ordered fluid equations and model finite-Larmor-radius and Landau-damping effects (Dorland & Hammett, 1993; Snyder et al., 1997). Gyrofluid codes have models for the treatment of dynamics on the ion-gyroradius scale and smaller, a regime in which drift-ordered fluid equations break down (Scott, 2007b, a). Both fluid and gyrofluid turbulence codes have yielded many insights into edge and SOL turbulence, but attempts to make quantitative comparisons with experimental data from tokamaks are rare and have produced mixed results (Zweben et al., 2009; Halpern et al., 2017, 2015; Cohen et al., 2013). These codes keep just a few moments and cannot fully capture potentially important kinetic effects, such as trapped particles (Lackner et al., 2012), nonlinear wave–particle interactions, and suprathermal electrons, and their model assumptions can be violated in edge and SOL plasmas (Batishchev et al., 1997; Takizuka, 2017). While fluid and gyrofluid models have been useful in revealing the qualitative physics of the boundary plasma, satisfactory and reliable quantitative prediction of boundary-plasma properties are believed to require the use of kinetic simulations in some capacity (Cohen & Xu, 2008; Scott et al., 2010; Scott, 2003), including the direct or indirect coupling of a fluid transport code to kinetic turbulence code (Schneider & Runov, 2007).

For these reasons, there are efforts to develop first-principles gyrokinetic codes for boundary-plasma turbulence simulation (Chang et al., 2009; Shi et al., 2017; Dorf et al., 2016; Korpilo et al., 2016). Unlike drift-reduced Braginskii-fluid approaches, gyrokinetic approaches use equations that are valid across a wide range of collisionality regimes, even if the collisional mean free path is not small compared to the parallel scale length or if the ion drift-orbit excursions are not small compared to radial gradient length scales (Cohen & Xu, 2008). Gyrokinetic simulations, however, are much more computationally expensive than fluid simulations, so fluid-based transport and turbulence codes for the boundary plasma will remain useful for modeling the boundary plasma. The results from gyrokinetic simulations are also expected to aid in improving the fidelity of boundary-plasma fluid simulations (Ricci, 2015).

#### 2.2 Gyrokinetic Modeling

Full six-dimensional kinetic modeling of plasma turbulence in tokamaks on macroscopic time and length scales by solving the Vlasov–Maxwell or Vlasov–Poisson equations have memory and processing power requirements several orders beyond what is currently possible on present-day and near-term supercomputers. Fortunately, there is a way to reduce the often insurmountable full 6D problem to a tractable 5D one, given that certain assumptions are well satisfied. Gyrokinetic theory is a reduced five-dimensional description of low-frequency plasma dynamics constructed by systematically removing the details of the charged particles’ rapid gyromotion due to a magnetic field and other high-frequency phenomena (Krommes, 2012; Tronko et al., 2016; Krommes, 2010; Sugama, 2000; Brizard & Hahm, 2007; Brizard, 2000a). This time-scale separation is well justified for particles in a strong background magnetic field with weak spatial inhomogeneity, such as those present in tokamaks and stellarators, in which case the frequencies of turbulent fluctuations are much smaller than the ion gyrofrequency. The gyrokinetic system, which describes the evolution of a gyrocenter (the gyro-averaged particle position) distribution function over a 5D phase space, is much easier to simulate when compared to 6D kinetic descriptions of particle distribution functions because of the reduced dimensionality and from relaxing the restriction on the time step from the plasma period to turbulence time scales and the restriction on the grid spacing from the Debye length to the gyroradius (Garbet et al., 2010).

While 5D gyrokinetic simulations require much more computational resources than comparable 3D fluid simulations due to the high dimensionality of the gyrokinetic system, their use to study turbulent transport in the tokamak core has now become routine (Garbet et al., 2010). A number of important verification studies and cross-code benchmarks on core gyrokinetic codes have been performed (Dimits et al., 2000; Tronko et al., 2017; Lapillonne et al., 2010; McMillan et al., 2010; Rewoldt et al., 2007). More recently, gyrokinetic models have also been used to study astrophysical turbulence (Schekochihin et al., 2009; Numata et al., 2010). For many of the same reasons why quantitative modeling in the boundary plasma is difficult to approach analytically, gyrokinetic codes for boundary-plasma simulation are much less mature than than their core counterparts.

Some complications that must be faced by boundary-plasma codes include the need to handle large-amplitude fluctuations (invalidating conventional approaches (Hu & Krommes, 1994)), open and closed magnetic field lines with a LCFS and X-point (which can cause difficulties with coordinates), electromagnetic fluctuations (Scott, 2007b; Scott et al., 2010), a wide range of space and time scales, a wide range of collisionality regimes, sheath boundary conditions, plasma–wall interactions, atomic physics, and the existence of a high-frequency electrostatic shear Alfvén () mode in electrostatic simulations (Lee, 1987; Belli & Hammett, 2005) or sheath-interaction modes that one does not want to artificially excite. Major extensions to existing core gyrokinetic codes or new codes are required to handle the additional challenges of the edge and SOL regions.

**
Numerical Implementations of the Gyrokinetic Equations ** A variety of numerical methods have been developed for the computationally challenging
solution to the gyrokinetic equations.
The two main types of numerical methods for solving the gyrokinetic equations are continuum methods
(Jenko & Dorland, 2001) and the particle-in-cell (PIC) method (Lee, 1983; Bottino & Sonnendrücker, 2015; Birdsall & Langdon, 2004).
There is also a third approach called semi-Lagrangian methods, which are a hybrid
between continuum and PIC methods, but they have been seldom-used for gyrokinetic simulation
so far (Grandgirard et al., 2006, 2007).
While gyrokinetic simulations have been useful in elucidating the physical mechanisms behind
tokamak and stellarator microturbulence, the ultimate goal of gyrokinetic simulations
is to produce quantitatively reliable predictions of core and boundary plasma properties.
As the history of core gyrokinetic codes has demonstrated, it is important to explore both
PIC and continuum (and semi-Lagrangian) approaches as independent cross-checks against each other and to
continuously shore up the specific weaknesses of each method.

Continuum methods are Eulerian approaches to solve a kinetic equation (e.g., the 5D gyrokinetic equation) by discretizing the equation on a fixed phase-space mesh. Standard numerical methods developed for the solution of partial differential equations are used, including finite-difference, finite-volume, spectral, pseudospectral, finite-element, and discontinuous Galerkin (DG) methods. The PIC method is a Lagrangian approach that solves the kinetic equation using a finite set of particles called markers. Since it is often computationally infeasible to use a number of markers on the same order as the number of physical particles in a real plasma, a much lower number of markers are used in practice (Tskhakaya et al., 2007). Each marker in the simulation then represents a ‘macroparticle’ or ‘superparticle’ encapsulating many physical particles. Starting with a set of markers that sample the initial positions in phase space, the marker positions are advanced over a small time step according to the characteristics of the kinetic equation. The source terms for the 3D field equations are then computed on a fixed grid from the markers, and the resulting fields are interpolated back to the marker positions so that the markers can be advanced again for next time step. PIC methods can be considered as a kind of Monte Carlo method that uses a finite set of markers to approximate integrals involving the distribution function (Bottino & Sonnendrücker, 2015; Krommes, 2012). The 3D field equations in both continuum and PIC approaches are solved using standard grid-based algorithms.

Both classes of numerical methods are associated with a unique set of advantages, disadvantages, and challenges.
PIC methods automatically maintain the positivity of the distribution function and are relatively straightforward
to implement (Krommes, 2012).
They are also not subject to the Courant–Friedrichs–Lewy (CFL) condition when
explicit time stepping is used (Garbet et al., 2010),
which results in instability if violated in a continuum simulation.^{4}^{5}

Continuum methods are not Monte Carlo methods, so they avoid having to deal with the challenging
statistical-noise issues that PIC methods face.
These methods are however subject to a time-step size restriction from
the CFL condition if explicit time stepping is used, which can be highly restrictive
in electrostatic simulations (Lee, 1987; Belli & Hammett, 2005) and
for the treatment of diffusive and hyperdiffusive terms.
Therefore, the fastest dynamics of a system must always be resolved when explicit methods are used,
even if they do not affect the results.
Implicit or semi-implicit time stepping methods are required to avoid the CFL restriction, which can
require computationally expensive global matrix-inversion operations and processor communication.
Continuum methods usually do not automatically guarantee the positivity of the distribution function,
which can be a particularly difficult issue to address for high-order schemes (Garbet et al., 2010).
Continuum methods often use upwind techniques that introduce some numerical dissipation
that result in some smoothing of the particle distribution function at small scales, which
can be beneficial in addressing the ‘entropy paradox’ (Krommes & Hu, 1994) and numerical-recurrence issues
Garbet et al. (2010), which can affect the results if the grid is too coarse (Candy & Waltz, 2006).
The inclusion of electromagnetic effects in continuum codes has been a long-solved issue^{6}

Another difference between PIC and continuum codes is that PIC codes have increasingly good velocity space resolution for longer-wavelength modes, while continuum codes have a velocity-space resolution that is independent of spatial scale. The required velocity-space resolution is highly problem dependent, and the resonance broadening by nonlinear scattering (island overlap) and collisions should be accounted for. The effects described in Su & Oberman (1968) show that collisional diffusion sets a limit on velocity-space resolution requirements that scales as , where is the collision frequency, so the moderate collisionality of the boundary plasma often means that the velocity-space resolution does not need to be very high (see also Smith, 1997).

**
Gyrokinetic Codes for the Boundary Plasma** Currently, the most sophisticated gyrokinetic code for the boundary plasma is the XGC1 gyrokinetic
particle code (Chang et al., 2009),
which presently uses a ‘hybrid-Lagrangian scheme’ (Ku et al., 2016) and includes
a realistic diverted-plasma geometry, neutral particles, charge-exchange
and ionization interactions, and radiation cooling (Chang et al., 2017).
Recent electrostatic XGC1 simulations have predicted the midplane heat-flux width
on an attached ITER plasma to be 5.6 mm, which can be compared to the 5 mm
ITER design specification and the 1 mm
empirical scaling based on present experiments (Chang et al., 2017).
The authors were able to reproduce the measured midplane heat-flux widths in three major
tokamaks to build confidence in their ‘gyrokinetic projection’ to ITER.
The heat-flux width in XGC1 simulations of present-day tokamaks was found to be dominated
by ion-drift-orbit excursions, while the heat-flux width in the ITER simulation was
found to be dominated by turbulent electron-heat-flux spreading.
It remains to be seen what challenges the inclusion of electromagnetic effects in XGC1
simulations with kinetic electrons will present
and whether similar calculations can be obtained at much cheaper computational cost through
other codes or models (the ITER simulation
reported in Chang et al. (2017) used 300 billion markers and ran on 90% of the Titan
computer for a few days).
ELMFIRE (Korpilo et al., 2016; Heikkinen et al., 2008) is another 5D gyrokinetic PIC code being
extended to handle the boundary plasma, although it is still in an early stage of development.

On the other hand, the development of gyrokinetic continuum codes for the boundary plasma has lagged significantly behind gyrokinetic PIC codes, despite the promising review of three main efforts at the time by Cohen & Xu (2008). Unfortunately, negative results are seldom reported in science, and so it is unclear what issues have held up the development of these gyrokinetic continuum codes. Some hints are found in the literature, however. Recent papers reporting G5D simulations (Kawai et al., 2017; Idomura, 2014, 2012) focus only on core turbulence and no longer have indications that the code is being extended to the boundary plasma. The FEFI code (Scott et al., 2010; Scott, 2006a) tried to proceed directly to electromagnetic simulations in the SOL, but ran into difficulties arising from sheath-model stability and shear-Alfvén dynamics in low-density regions (Zweben et al., 2009). The development of TEMPEST (Xu et al., 2007) was apparently halted sometime after results from 4D gyrokinetic transport simulations were presented in Xu et al. (2010), perhaps due to issues stemming from the non-conservation properties of the underlying numerical scheme (Cohen & Xu, 2008). Some members of the TEMPEST team eventually began the development of the COGENT code (Dorf et al., 2012), which features a conservative fourth-order finite-volume discretization (Colella et al., 2011). COGENT was recently used to perform axisymmetric 4D gyrokinetic transport simulations in a realistic geometry with an anomalous radial-diffusion term to model radial transport due to turbulence (Dorf et al., 2016). A modified version of the GENE code (widely used for core turbulence simulation) has recently begun development for boundary-plasma applications (Pan et al., 2016). A gyrokinetic continuum code with similar capabilities of XGC1 would be extremely valuable to the plasma physics field, since XGC1 is the only gyrokinetic code capable of performing ITER boundary-plasma simulations at present. Having the same prediction made by two (or several) independent gyrokinetic codes for the boundary plasma using different numerical approaches would be highly reassuring.

In continuum codes, spectral techniques are commonly used in some directions, which can have problems with Gibbs phenomena that result in negative overshoots. Most algorithms used in magnetic fusion research are designed for cases in which viscous or dissipative scales are fully resolved and do not use limiters, and thus can have problems with small negative oscillations. Negative densities may result in unphysical behavior in the solution (for example, a negative density in the tail of the electron distribution function can reverse the slope of the sheath current versus sheath potential relation), and inaccuracies in the sheath boundary conditions can also lead qualitatively incorrect results. Some finite-difference algorithms make it easier to calculate derivatives across the LCFS with field-aligned coordinates, but may have problems with particle conservation, and small imbalances in electron and ion gyrocenter densities may drive large errors in the electric field.

For these reasons, we were motivated to investigate discontinuous Galerkin methods for gyrokinetic continuum simulation in boundary plasmas. DG methods are a class of finite-element methods that use discontinuous basis functions (typically piecewise polynomials) to represent the solution in each cell. The most popular version of DG is the Runge–Kutta discontinuous Galerkin method (Cockburn & Shu, 1998b, 2001; Shu, 2009), which uses a high-order DG method for space discretization and explicit, strong stability-preserving (SSP) high-order Runge–Kutta methods (Gottlieb et al., 2001) for time discretization. The RKDG method was originally developed for the solution of nonlinear, time-dependent hyperbolic systems and has found use in numerous applications (Cockburn et al., 2000), especially for Euler and Navier–Stokes equations. Since continuity in the solution is not required across cell interfaces, DG methods gain a number of important benefits that are not available to traditional finite-element methods. The RKDG method is attractive because it is highly local, highly parallelizable, able to handle complex geometries, allows high-order accuracy, and enforces local conservation laws. For an introduction to DG methods, the reader is referred to Shu (2009); Durran (2010); Hesthaven & Warburton (2008).

### 3 Thesis Overview

This thesis focuses on efforts towards the development of a gyrokinetic continuum code for the simulation of boundary-plasma turbulence. Specifically, we investigate the application of discontinuous Galerkin algorithms to handle the difficulties in open-field-line plasmas, which is the situation found in the SOL. The algorithms that we use have been implemented in the Gkeyll code, which is a framework for kinetic and fluid plasma simulations using a variety of grid-based numerical algorithms. The Gkeyll code is primarily developed at the Princeton Plasma Physics Laboratory, with contributors from a variety of institutions around the United States. We note that we have not yet performed a ‘gyrokinetic continuum simulation of a tokamak SOL’. At a minimum, a tokamak SOL simulation needs to include a confined edge region where plasma is sourced and a realistic diverted geometry including a LCFS and X-point. Additionally, kinetic modeling of wall-recycled neutrals and models of radiative power losses, charge-exchange interactions, and ionization are required for quantitative prediction, since these processes play important energy-dissipation roles in the SOL. Nevertheless, major steps towards this goal have been completed in the course of this thesis.

In Chapter \thechapter, we discuss the important models and algorithms used in other chapters of the thesis. We first describe the gyrokinetic model that has been implemented in Gkeyll, which at present employs a number of simplifications (electrostatics, long-wavelength, linear polarization) to make the problem tractable in the scope of a PhD thesis. Next, we discuss discontinuous Galerkin algorithms, starting from a 1D example before covering the specific energy-conserving version of DG that we use in our simulations. We also discuss important aspects of a simplified Lenard–Bernstein collision operator and how positivity issues in the distribution function are addressed. A key component of simulating plasma dynamics on open magnetic field lines is the sheath-model boundary condition applied at the material interfaces, and we describe two kinds of sheath models that can be used and note an important shortcoming of these models concerning the Bohm sheath criterion that is not usually acknowledged.

In Chapter \thechapter, we present results from our initial efforts to investigate the feasibility of using DG methods for gyrokinetic continuum simulation in the boundary plasma in spatially 1D kinetic simulations. We describe the construction of a simplified 1D1V (one position dimension, one velocity dimension) gyrokinetic model that incorporates ion polarization effects through a specified perpendicular wavenumber in a modified gyrokinetic Poisson equation. Combined with logical-sheath boundary conditions, this model is then used to simulate the parallel propagation of an ELM heat pulse in the SOL, which is a problem that has been studied before using kinetic simulations that fully resolved the sheath and fluid simulations that used sheath boundary conditions. This model is then extended to 1D2V and some collisional effects are added through a Lenard–Bernstein collision operator. Despite not directly resolving the sheath, our 1D1V and 1D2V gyrokinetic simulations agree quantitatively well with comparable fully kinetic simulations that are much more computationally expensive due to restrictive spatial and temporal resolution requirements.

Chapter \thechapter presents the key accomplishment of the thesis, which are the first 5D gyrokinetic continuum simulations of turbulence in a straight-magnetic-field-line geometry. Specifically, we present simulations of the Large Plasma Device (LAPD), which is a basic plasma physics experiment at the University of California, Los Angeles. Compared to a realistic scrape-off layer, the LAPD plasma is at a much colder temperature and is not subject to magnetic-curvature effects. The LAPD plasma is very well diagnosed, so we compare turbulence characteristics from our simulations to previous LAPD measurements published elsewhere. We also describe a simple modification to the sheath-model boundary conditions that allows us to simulate a set of LAPD experiments to investigate sheared-flow-suppression of turbulence using bias-induced flows.

In Chapter \thechapter, we add some additional complexity to the open-field-line simulations by adding magnetic-curvature effects to simulate turbulence on helical field lines. While the new magnetic geometry makes the simulations particularly suitable for simulating the plasma turbulence in SMTs, we design a test case for a helical SOL using parameters relevant for the NSTX SOL. The helical-SOL simulations are qualitatively different from the LAPD simulations, with the generation and radial propagation of blobs playing an important role in transporting plasma across the magnetic field lines. We also show how the magnetic-field-line incidence angle affects plasma profiles and turbulence characteristics in these simulations. In this simple model, we show that the heat-flux width is strongly affected by the strength of the vertical (poloidal) magnetic field.

In Chapter \thechapter, we present a numerical method that uses exponentially weighted polynomials to represent the solution while maintaining important conservation properties. The use of non-polynomial basis functions is motivated by the need to use as few pieces of data to represent the solution as possible, and exponentially weighted polynomials appear to be a reasonable choice to represent the distribution function in problems in which collisions are strong. Previous work (Yuan & Shu, 2006) in using non-polynomial basis functions in standard DG methods does not conserve number, momentum, and energy. Using 1D numerical tests of collisional relaxation, we show how the new method conserves important quantities to machine precision. Results from a non-trivial calculation of the parallel heat flux in a simplified Spitzer–Härm test problem are then presented to compare the accuracy and efficiency of the new method with standard DG methods using polynomials. The generalization of this method to high dimensions and the implementation of this method in Gkeyll is left for future work.

In Chapter \thechapter, we summarize the main results of this thesis and discuss what are high-priority directions for near-term future work.

#### 3.1 A Note on Color Maps

Due to the continued popularity of the rainbow color map for representing interval data in plasma-physics papers, we briefly explain why authors should avoid using the rainbow color map for such purposes. We decided to place this discussion here in the introduction chapter instead of in the appendices to help spread general awareness of this issue. In fact, even we are guilty of using a rainbow color map for the 2D plots in a recent publication (Shi et al., 2017) before this issue was brought to our attention.

As done by most authors, we present 2D data using pseudocoloring, in which data is displayed by mapping scalar data values to colors according to a color map. While many authors still use a color map that is ordered according to the visible light spectrum, known as the rainbow color map, data visualization experts have long recognized that the rainbow color map is confusing and misleading (Eddins, 2014; Borland & Taylor II, 2007; Ware, 1988; Rogowitz & Treinish, 1998). Many of these issues arise from the lack of perceptual ordering and perceptual uniformity in the rainbow color map. Perceptual ordering refers to a color map that uses a sequence of colors with a consistent, inherent ordering: Given a set of distinct colors, will most people order the color in the same sequence based on their color perception (and not a mnemonic)? Many issues in the rainbow color map come from our perception of yellow as the brightest color, while the rainbow color map assigns the largest data values to red. Perceptual uniformity refers how the same difference between two data values corresponds to the same perceived difference in color on the entire color scale. The perceived difference between the colors representing the values 1 and 2 should also be the same between the colors representing the values 9 and 10.

The color map used for 2D data visualization in this thesis is the ‘inferno’ color map (Smith & van der Walt, 2016), which is available Matlab, Matplotlib, and R. The inferno color map is perceptually uniform and perceptually ordered. Figure 5 shows a comparison of electron-density snapshots from a gyrokinetic simulation of LAPD (discussed in Chapter \thechapter) using the rainbow color map (called the jet color map in Matlab) and the inferno color map. By comparing the two sets of plots, several visual artifacts can be identified in the plots with the rainbow color map. Additionally, the detail in green and cyan regions is wiped out due to the perceptual similarity of these two colors.

## Chapter \thechapter Models and Numerical Methods

In this chapter, we discuss some of the common models and numerical algorithms used elsewhere in this thesis. Although it is well known that scientific results must be reproducible, authors sometimes do not describe their simulations in published papers in enough detail for other researchers to replicate their results. One should strive to publicly document the model equations, initial conditions, boundary conditions, grid resolution, and any special techniques that were necessary to obtain stable simulations. We hope to facilitate any future attempts to reproduce our results and to break away from the tendency to be vague about the ‘ugly parts’ of a simulation by providing a clear description of our approach for the open-field-line simulations discussed in Chapters \thechapter and \thechapter, including the difficulties we encountered and what was done to address them.

In Section 4, we describe the gyrokinetic model we solve in Chapters \thechapter and \thechapter. We then provide an overview of discontinuous Galerkin (DG) methods and discuss how we apply an energy-conserving DG method (Liu & Shu, 2000) to solve the gyrokinetic system in Section 5. We provide details about the numerical implementation of the Lenard–Bernstein collision operator in Section 6, which also motivates the need for a positivity-adjustment procedure described in Section 7. Lastly, we specify the sheath-model boundary conditions that are used in our simulations in Section 8.

These models and algorithms are implemented in the Gkeyll code, which is a framework for kinetic and fluid plasma-physics problems. Recently, Gkeyll has been used for fluid studies of magnetic reconnection (Wang et al., 2015; Ng et al., 2015), kinetic simulations of the Vlasov–Maxwell system (Juno et al., 2017), and kinetic and multi-fluid sheath modeling (Cagas et al., 2017). Gkeyll is composed of a core component written in C++ and problem-specific configuration components written in the Lua scripting language (Ierusalimschy et al., 1996).

Gkeyll was originally designed and developed at the Princeton
Plasma Physics Laboratory (PPPL) by A. Hakim, who created the core framework
components (e.g., Lua integration, data structures, and grid classes)
used in the various applications.
For the simulations presented in this thesis, A. Hakim added additional
domain-decomposition capabilities, developed the gyrokinetic-Poisson-equation solver,
and contributed to the extension of the energy-conserving Liu & Shu algorithm (Liu & Shu, 2000)
to generic Hamiltonian systems (Section 5.2).
There are currently several code contributors from various academic institutions
across the United States.
At the time of writing, the Gkeyll repository is stored on Bitbucket and
access must be requested.^{7}

The core C++ component contains implementations of the basic self-contained classes, such as data structures to store solutions, rectangular finite-element meshes, and solvers, which generically perform operations on input data to produce output data. Currently, Gkeyll has classes that implement discontinuous Galerkin, finite-element, finite-volume, and finite-difference algorithms. Most of the work in this thesis required the addition of new capabilities to the Gkeyll code in the form of solvers.

The configuration component uses Lua for problem-specific applications of the Gkeyll code. The core Gkeyll code only needs to be compiled by the user once, while a Lua script is provided as an input argument to the executable and is automatically compiled at run time. Lua is used for much more than supplying parameter values for a simulation. In a Lua script, a user defines which objects (implemented in the core code) are to be created for a simulation and how these objects interact with each other. To elaborate, there is no implementation of an ‘open-field-line simulation’ in the core C++ code. There are several classes relevant to this problem, such as implementations of sheath boundary conditions, basis functions, a Poisson-equation solver, a kinetic-equation solver, and more, but a Lua script is required to connect all these pieces together and determine the sequence in which various tasks are performed. Lua scripts for the 5D gyrokinetic continuum simulations have approximately 3000 lines of code.

### 4 Gyrokinetic Model

The gyrokinetic model discussed in this section is used in simulations of the Large Plasma Device (LAPD) (Chapter \thechapter) and a helical scrape-off-layer (SOL) model (Chapter \thechapter). Several versions of full- gyrokinetic equations have been derived with various formulations, ordering assumptions, and levels of accuracy (Brizard & Hahm, 2007; Sugama, 2000; Hahm et al., 2009; Parra & Calvo, 2011; Parra et al., 2014; Dimits, 2012; McMillan & Sharma, 2016), and can generically be written in the form of an evolution equation for the gyrocenter distribution function, with a Poisson bracket for the phase-space velocities and expressions for the Lagrangian/Hamiltonian, coupled to field equations to determine the potentials. Here, we solve a long-wavelength (drift-kinetic) limit of electrostatic full- gyrokinetic equations with a linearized polarization term for simplicity, as summarized by Idomura et al. (2009). As the code is further developed, it can be extended to more accurate and more general equations, though the equations will still have this generic structure.

The fundamental assumption of standard gyrokinetics is that there is a coordinate system in which things change slowly compared to the gyrofrequency. In some gyrokinetic derivations, the ordering assumptions are written in a more restrictive form requiring that the fluctuation amplitudes must be small. But as discussed in various places (such as Hahm et al., 2009; Dimits, 2012; McMillan & Sharma, 2016), more general derivations that are appropriate for the edge region of fusion devices are possible, such as using a small vorticity ordering (McMillan & Sharma, 2016), which allows large flows and large-amplitude fluctuations at long wavelengths.

We solve a full- gyrokinetic equation written in the conservative form (Brizard & Hahm, 2007; Sugama, 2000; Idomura et al., 2009)

(1) |

where is the gyrocenter distribution function for species , is the Jacobian of the gyrocenter coordinates, , , represents the effects of collisions, , and represents plasma sources (e.g., neutral ionization or core plasma outflow). In a straight-field-line geometry, simplifies to . The phase-space advection velocities are defined as and , where the gyrokinetic Poisson bracket is

(2) |

The gyrocenter Hamiltonian is

(3) |

where is the gyro-averaged potential (with the gyro-angle denoted by ). In the simulations discussed in Chapters \thechapter and \thechapter, we consider a long-wavelength limit of the gyrokinetic system and neglect gyroaveraging in the Hamiltonian to take . This system has similarities to some versions of drift kinetics (and is sometimes referred to as the drift-kinetic limit of gyrokinetics (Dorf et al., 2016, 2013; Cohen & Xu, 2008)), but is unlike versions that include the polarization drift in the kinetic equation or determine the potential from some other equation. In a straight-magnetic-field geometry, (1)–(3) reduce to the description of parallel streaming, an drift, and acceleration along the field line due to (see (126), which is solved in our LAPD simulations).

The potential is solved for using the long-wavelength gyrokinetic Poisson equation with a linearized ion polarization density

(4) |

where , , and is the background ion gyrocenter density that we will take to be a constant in space and in time. Gyroaveraging in the gyrocenter densities is neglected in this long-wavelength limit. The replacement of by on the left-hand side of (4) is analogous to the Boussinesq approximation employed in some Braginskii fluid codes (Dudson et al., 2015; Halpern et al., 2016; Angus & Umansky, 2014). We note that the use of a linearized ion polarization charge density is formally valid when ion density fluctuations are small.

Note that (4) is a statement of quasineutrality, where the right-hand side is the gyrocenter component of the charge density , and the left-hand side is the negative of the ion polarization charge density, (due to the plasma response to a cross-field electric field), so this equation is equivalent to . The simulations are done in a Cartesian geometry with and being used as coordinates perpendicular to the magnetic field, which lies solely in the direction. Therefore, .

### 5 Discontinuous Galerkin Algorithms

Consider a numerical method to solve the following time-dependent partial differential equation over a domain :

(5) |

where is an operator that involves spatial derivatives in and an initial condition and boundary conditions are prescribed. One class of numerical methods, called series-expansion methods, solves (5) by approximating as a linear combination of a finite number of predetermined basis functions:

(6) |

In almost all practical cases of interest, it is not possible to find a solution that satisfies (5) by this expansion because the ’s are generally not eigenfunctions of (Durran, 2010), so the sign is used to relate to the numerical solution . Therefore, we settle for finding the degrees of freedom that minimize some kind of error that quantifies the degree to which the numerical solution fails to satisfy (5). For this purpose, it is convenient to define the residual

(7) |

and try to minimize a function involving the residual in an integral sense over or at a set of points in .

In Galerkin methods, a system of equations for the time evolution of the degrees of freedom, , is obtained by requiring that the residual be orthogonal to each basis function:

(8) |

The solution to (8) determines the time evolution of the degrees of freedom such that the squared--norm error

(9) |

is minimized. To see this, we substitute the expansion (6) into (9), take a derivative with respect to , and look for critical points:

(10) | ||||

(11) | ||||

(12) |

This choice of minimizes the squared--norm error because the second derivative of with respect to is positive:

(13) |

The Galerkin approximation (8) is used in many series-expansion methods, including the spectral method, some finite-element methods, and the discontinuous Galerkin method. The Runge–Kutta discontinuous Galerkin (RKDG) method (Cockburn & Shu, 1998b, 2001; Shu, 2009) is a semi-discrete numerical method that uses a discontinuous Galerkin discretization for the spatial variables and explicit high-order-accurate Runge–Kutta methods (made of convex combinations of first-order Euler steps) for time discretization. The method is particularly well-suited for the solution of nonlinear, time-dependent hyperbolic conservation laws and has found use in numerous applications (Cockburn et al., 2000). As its name implies, the discontinuous Galerkin method is a Galerkin method that uses discontinuous basis functions. In contrast to other Galerkin methods, the condition (8) is enforced element-by-element instead of globally over the entire domain.

The use of discontinuous basis functions that are smoothly varying within a cell but zero everywhere outside
of it enables several benefits for DG methods.
Computations are highly localized in the sense that data only needs to be shared with immediate neighbors
regardless of the basis function degree,^{8}

#### 5.1 DG for 1D Conservation Laws

We review the RKDG method for a 1D nonlinear conservation law

(14) |

As in finite-element methods, we first partition the domain into a number of cells , for and define an approximation space for the discrete solution

(15) |

where the local space is usually taken to be , the space of polynomials up to degree for . Local basis functions that span are required for numerical implementation, and the solution in is expressed as

(16) |

so the RKDG method prescribes a way to solve for the evolution of the degrees of freedom

(17) |

Legendre and Lagrange polynomials are typically used as basis functions in modal and nodal DG representations, respectively. We have also implemented Serendipity basis functions (Arnold & Awanou, 2011), which are an attractive alternative to Lagrange polynomials. The Serendipity finite-element space has fewer basis functions (smaller dimension) than the Lagrange finite-element space, while achieving the same convergence rate. The Serendipity finite-element space can also represent solutions that are continuous across elements. The disparity between the size of these two finite element spaces grows with the degree of the finite element space and the space dimension. For example, degrees of freedom are required to represent the solution in 5D for a second-degree Lagrange finite element, while only degrees of freedom are needed for a second-degree Serendipity finite element, and one could expect a factor of five speedup when using Serendipity finite elements in this case due to prevalence of matrix operations of size in DG methods. We have not yet explored the use of higher-order Serendipity elements for 5D gyrokinetic simulations and use first-order elements (where the two finite-element spaces are identical) in those simulations for simplicity. DG methods also allow for the use of non-polynomial basis functions (Yuan & Shu, 2006), which can also result in significant savings over polynomial basis functions. These ideas are explored in Chapter \thechapter.

Recalling that Galerkin methods minimize the squared--norm error by requiring that the residual be orthogonal to the basis functions (8), we multiply (14) by an arbitrary test function and integrate over a cell

(18) |

where an integration by parts was performed to move the spatial derivative on onto the test function . Next, the exact solution is replaced by the numerical solution , the flux evaluated at the element interfaces is replaced by the numerical flux , and the test function is replaced by :

(19) |

The and superscripts indicate that a discontinuous function is evaluated at the interface using the left and right limits of the discontinuous numerical solution, respectively. In (19), the test functions are evaluated inside .

The numerical flux is an approximation to the exact flux and depends on the value of on each side of a boundary:

(20) |

The single-valued numerical flux is only defined at cell interfaces, where is discontinuous. The numerical flux does not appear in the volume integrals because there is no ambiguity in the value of to use in the interior of a cell. The choice of numerical flux must be consistent with the physical flux that it approximates: . The numerical flux must also be a non-decreasing function of its first argument and a non-increasing function of its second argument, which are required for the numerical scheme to reduce to a monotone finite-volume scheme when the solution is approximated by a constant in each cell (Cockburn & Shu, 2001). A standard choice for the numerical flux is to use the upwind numerical flux. If , then the upwind flux is

(21) |

It is important to use a numerical flux suitable for the problem at hand because the choice affects the approximation quality (Cockburn & Shu, 2001).

We can now obtain a system of coupled equations that can be solved for the unknowns in each cell by substituting the basis function expansion for (16) and taking for in (19):

(22) |

where are the components of the mass matrix for . Equation (22) can be written as a matrix equation , where the element of the vector is the right hand side of (22) evaluated for . In this example, is independent of time and is the same for all , and so it is advantageous to compute once and store its inverse at the beginning of a simulation.

Integrals are numerically evaluated using Gaussian quadrature methods. In the Gkeyll code, Gauss–Legendre quadrature is typically used to approximate the definite integral of a function from to by taking a weighted sum of the evaluation of the function to be integrated at a set of points:

(23) |

where the approximation is exact when is a polynomial of degree or less. On the interval , the Gauss–Legendre quadrature nodes are the zeros of , the th Legendre polynomial for , and the associated weights are computed as (Press et al., 2007)

(24) |

The quadrature rule on the interval can be scaled to an arbitrary interval using the relation

(25) |

so the quadrature rule becomes (Press et al., 2007)

(26) |

We note that other quadrature rules with different weights and nodes can be constructed for general integrals of the form , where is an arbitrary (possibly non-polynomial) weight function (for an investigation of several choices, see Landreman & Ernst, 2013). Gaussian quadrature rules for integrals involving multiple directions are simply computed by taking tensor products of the 1D quadrature rule.

In order to apply Gaussian quadrature methods to evaluate integrals as they appear in (22),
auxiliary matrices^{9}

(27) |

from the basis functions whose analytical form is known.

We also compute and store a matrix that calculates in terms of the same basis functions used to represent the spatial dependence of , which is used when the evaluation of spatial derivatives is required (such as in the volume integral on the right-hand side of (22)). That is, we seek the weights in element such that

(28) |

We multiply by each basis function , substitute the basis function expansion for , and integrate over to get the system of equations that can be solved for the :

(29) |

By defining an additional auxiliary matrix whose elements are

(30) |

we see that .

Finally, once (22) has been solved for the weights , we use an explicit high-order strong-stability-preserving (SSP) time stepping method to advance the solution in time. These methods achieve high-order accuracy by taking convex combinations of first-order forward Euler steps. Specifically, the three-stage, third-order SSP Runge–Kutta method (Gottlieb et al., 2001) is used in the simulations discussed in this thesis and is the most popular choice (Shu, 2009) for solving equations of the form

(31) |

where is the space discretization of the operator in (14). If the first-order Euler time discretization of (31) is stable under a certain norm for a sufficiently small time step (e.g., satisfying the Courant–Friedrichs–Levy condition), high-order SSP methods are designed to automatically maintain the strong stability property for certain higher-order time discretizations (Gottlieb et al., 2001) under a possibly more restrictive time step. To advance to , the third-order SSP Runge–Kutta method is

(32) | ||||

(33) | ||||

(34) |

In this method, the time step restriction is the same as the time step needed for the forward-Euler method to be strongly stable. One downside to this method is the considerable storage requirement. At any given instant of the algorithm, three versions of need to be stored. For two plasma species (electrons and one ion species), the storage requirement is doubled. It might be useful to consider the use of the low-storage, third-order SSP Runge–Kutta method, which only needs to store two versions of , but requires a 3.125 times smaller time step (Gottlieb et al., 2001).

#### 5.2 DG for the Gyrokinetic System

An energy-conserving (in the continuous-time limit) discontinuous Galerkin algorithm (Liu & Shu, 2000) is used to discretize the equations in space. Although Liu & Shu (2000) presented their algorithm for the 2D incompressible Euler and Navier–Stokes equations, G. Hammett recognized the general applicability of their algorithm for Hamiltonian systems, and A. Hakim contributed to the generalization. Upwind interface fluxes (21) are used in the discretization of the gyrokinetic equation (1). This algorithm requires that the Hamiltonian be represented on a continuous subset of the basis set used to represent the distribution function. Therefore, the distribution function is represented using discontinuous () polynomials, while the electrostatic potential is represented using continuous () polynomials (equivalent to continuous finite elements). In addition to contributing to the formulation of the generalized algorithm, A. Hakim also implemented and performed a number of two-dimensional tests (e.g. incompressible Euler equations, 1D1V Vlasov–Poisson system) to benchmark the convergence and conservation properties of this algorithm. Additional details regarding the generalization of the Liu & Shu (2000) algorithm will be provided in future paper by A. Hakim et al.

Liu & Shu (2000) presented their algorithm for 2D incompressible Euler equations in a
vorticity–stream function formulation: ,
with the stream function given by ,
and the velocity .
They showed analytically that the DG space discretization of these equations conserves
energy if the basis functions for are in a continuous subspace of the basis
functions used for the vorticity .^{10}

Here, we consider how the algorithm of Liu & Shu (2000) is applied to the gyrokinetic equation system (1)–(4), although we will neglect the collisions and source terms in the gyrokinetic equation (1) to focus on the Hamiltonian part of the system. First, we discretize the -dimensional phase-space domain by dividing it into a number of elements . We also define the configuration-space domain and its respective partition .

Next, we define the approximation spaces (following the notation of Liu & Shu, 2000)

(35) | ||||

(36) |

where is the space of polynomials in variables with each variable degree at most for .

We note that the gyrokinetic equation (1) can be written in the general form

(37) |

where the coordinates and , The Poisson bracket can be written as

(38) |

so . The Poisson matrix is assumed to be antisymmetric. As a consequence of this antisymmetry, the component of the characteristic velocities normal to a surface is continuous on that surface. This means that evaluated at an interface between two cells has the same value when approached from either side of the interface. Equivalently, this means that evaluated on the interface between two cells can be correctly evaluated using the solution from only one of the cells. To see this, we compute

(39) |

which is a term that will appear in surface integrals when we perform an integration by parts on the weak form of the advection equation. We also see that is orthogonal to :

(40) |

Since we have required that the Hamiltonian be continuous, we see from (39) and (40) that is continuous on cell surfaces. While can be discontinuous at cell surfaces because is only required to be continuous, will be continuous on a surface if is also continuous on that surface.

As in standard DG methods, we multiply the gyrokinetic equation (1) by a test function and integrate over each cell :

(41) |

By taking and summing over all cells, we see that particle number is conserved by the discrete scheme (assuming that the constant 1 is in the approximation space):

(42) |

For the 5D gyrokinetic system, has the following form (see (2)):

(43) |

Following a standard finite-element approach, in order to satisfy continuity constraints on the potential, we project the gyrokinetic Poisson equation onto the space of basis functions for the potential, yielding a single non-local 3D solve to find the potential. Multiplying (4) by a test function and integrating over the entire configuration-space domain (not just over an individual cell, since the basis functions for couple the cells together):

(44) |

where . There are ways to solve (44) by performing a set of independent 2D solves, followed by a local self-adjoint smoothing/interpolation operation, but we have not yet implemented this potentially more efficient approach. Although a non-local solve in three dimensions is required for the potential, the 5D gyrokinetic equation itself can be solved in a highly local manner.

We can see that there will be a conserved energy by taking to be the discrete Hamiltonian:

(45) |

The third term is zero because . By summing (45) over all elements, the second term vanishes because the surface integrals not on the domain boundaries appear in equal and opposite pairs due to the continuity of , while the surface integrals on the domain boundaries are zero by boundary conditions (e.g., zero flux or ). Therefore, we have

(46) |

To identify the conserved energy for the gyrokinetic system, we insert the discrete Hamiltonian into (46) and sum over all species:

(47) |

Carrying out the velocity-space integrals and identifying the first two terms as a thermal energy, we have

(48) |

where and the species summation on the last term was carried out to write it in terms of the gyrocenter charge density . Next, we take a time derivative of (44) and set to get

(49) |

where the surface terms are zero from the boundary conditions (e.g., periodic or on the side walls). Therefore, we have the energy conservation law for the discrete system

(50) |

where the energy .

For simplicity, we use nodal, linear basis functions to approximate the solution in each element for the 5D gyrokinetic simulations in Chapters \thechapter and \thechapter. This choice leads to degrees of freedom per cell in the 5D phase-space mesh ( degrees of freedom in the 3D configuration-space mesh). With the degrees of freedom specified in a cell, can be evaluated anywhere within the cell without additional approximation. With knowledge of the basis functions, we can generate data for plotting without having to rely on approximate interpolation methods. This means that the data we show in various plots in this thesis have not been filtered or smoothed in post processing to make a more visually attractive image. Some more details about the plotting procedures used to create figures for this thesis are provided in Appendix \thechapter.

The choice of a linear-polynomial basis set means that , which appears in the Hamiltonian (3), cannot be exactly represented. We therefore approximate in the linear basis set by requiring that the piecewise-linear approximation equal at the DG nodes. These nodes are located on the vertices of uniform rectangular cells, so the nodes of cell either have or , where is the coordinate of the center of cell . By approximating in such a manner, will vary linearly in cell from to . Note that the approximated is continuous across elements, as required for numerical energy conservation, and its first derivative is discontinuous across elements, i.e. . We use rectangular meshes with uniform cell spacing, but we note that most of the DG algorithms discussed in this thesis are trivially generalizable to non-uniform and non-rectangular meshes.

### 6 Collision Operator

Electron–electron and ion–ion collisions are implemented using a Lenard–Bernstein model collision operator (Lenard & Bernstein, 1958)

(51) |

where standard expressions are used for collision frequency (Huba, 2013, p. 37), , and . The exact expressions used in our simulations for and are (Fitzpatrick, 2011)

(52) | ||||

(53) |

where is Coulomb logarithm for expressed in and