Electromagnetic full-$f$ gyrokinetics in the tokamak edge with discontinuous Galerkin methods
N. R. Mandell, A. Hakim, G. W. Hammett, M. Francisquez
TL;DR
This paper introduces an energy-conserving discontinuous Galerkin scheme for full-$f$ electromagnetic gyrokinetics in tokamak edge plasmas, successfully performing nonlinear simulations on open field lines without Ampère cancellation issues.
Contribution
It develops the first nonlinear electromagnetic gyrokinetic simulation on open field lines using a novel discontinuous Galerkin method that conserves energy and avoids Ampère cancellation.
Findings
Scheme successfully avoids Ampère cancellation problem.
First nonlinear electromagnetic gyrokinetic simulation on open field lines.
Comparison shows differences between electromagnetic and electrostatic results.
Abstract
We present an energy-conserving discontinuous Galerkin scheme for the full- electromagnetic gyrokinetic system in the long-wavelength limit. We use the symplectic formulation and solve directly for , the inductive component of the parallel electric field, using a generalized Ohm's law derived directly from the gyrokinetic equation. Linear benchmarks are performed to verify the implementation and show that the scheme avoids the Amp\`ere cancellation problem. We perform a nonlinear electromagnetic simulation in a helical open-field-line system as a rough model of the tokamak scrape-off layer using parameters from the National Spherical Torus Experiment (NSTX). This is the first published nonlinear electromagnetic gyrokinetic simulation on open field lines. Comparisons are made to a corresponding electrostatic simulation.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Electromagnetic full- gyrokinetics in the tokamak edge with discontinuous Galerkin methods
N. R. Mandell\aff1 \corresp
A. Hakim\aff2
G. W. Hammett\aff2
M. Francisquez\aff3
\aff1Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA \aff2Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USA \aff3MIT Plasma Science and Fusion Center, Cambridge, MA, 02139, USA
Abstract
We present an energy-conserving discontinuous Galerkin scheme for the full- electromagnetic gyrokinetic system in the long-wavelength limit. We use the symplectic formulation and solve directly for , the inductive component of the parallel electric field, using a generalized Ohm’s law derived directly from the gyrokinetic equation. Linear benchmarks are performed to verify the implementation and show that the scheme avoids the Ampère cancellation problem. We perform a nonlinear electromagnetic simulation in a helical open-field-line system as a rough model of the tokamak scrape-off layer using parameters from the National Spherical Torus Experiment (NSTX). This is the first published nonlinear electromagnetic gyrokinetic simulation on open field lines. Comparisons are made to a corresponding electrostatic simulation.
††journal: Journal of Plasma Physics
1 Introduction
Understanding turbulent transport physics in the tokamak edge and scrape-off layer (SOL) is critical to developing a successful fusion reactor. The dynamics in these regions plays a key role in determining the L–H transition, the pedestal height, and the heat load to the vessel walls. While the edge is often modelled by Braginskii-type fluid models that have provided valuable results and insights (Xu et al., 2008; Tamain et al., 2010; Ricci et al., 2012; Zhu et al., 2017; Francisquez et al., 2017), a kinetic treatment will inevitably be necessary for reliable quantitative predictions in some cases (Jenko & Dorland, 2001; Cohen & Xu, 2008). Gyrokinetic theory and direct numerical simulation have become important tools for studying turbulence and transport in fusion plasmas, especially in the core region (Parker et al., 1993a; Kotschenreuther et al., 1995; Lin et al., 2000; Dimits et al., 2000; Dorland et al., 2000; Jenko, 2000; Jost et al., 2001; Candy & Waltz, 2003; Idomura et al., 2003; Watanabe & Sugama, 2005; Jolliet et al., 2007; Idomura et al., 2008; Peeters et al., 2009; Lanti et al., 2019). In the edge and SOL, gyrokinetic simulations are particularly challenging because the large, intermittent fluctuations in the SOL make assumptions of scale separation between equilibrium and fluctuations not strongly valid. This necessitates a full- approach that self-consistently evolves the full distribution function, (as opposed to the approach commonly used in the core, where one assumes with a fixed background so that only perturbations must be evolved, and the parallel electric field nonlinearity is frequently neglected). Steady progress in gyrokinetic edge/SOL modelling has been made with both particle-in-cell (PIC) (Ku et al., 2009; Korpilo et al., 2016; Ku et al., 2016) and continuum (Shi et al., 2019; Pan et al., 2018; Dorf et al., 2016) methods. Another challenge is the magnetic geometry of the edge/SOL region, which requires treatment of open and closed magnetic field-line regions and the resulting plasma interactions with material walls on open field lines. The X-point in a diverted geometry is an additional complication which makes the use of field-aligned coordinates challenging. Currently only the XGC1 hybrid-Lagrangian PIC code (Ku et al., 2016) can simulate gyrokinetic turbulence in a three-dimensional diverted geometry with an X-point.
The edge/SOL region also features steep pressure gradients, especially in the H-mode transport barrier and SOL regions, which contributes to the importance of electromagnetic effects. In this regime, the parallel electron dynamics is no longer fast relative to the drift turbulence, so electrons can no longer be treated adiabatically (Scott, 1997). This leads to coupling of the perpendicular vortex motions and kinetic shear Alfvén waves, which results in field-line bending (Xu et al., 2010). Including electromagnetic effects in gyrokinetic simulations has proved numerically and computationally challenging, both in the core and in the edge. The so-called Ampère cancellation problem is one of the main numerical issues that has troubled primarily PIC codes (Reynders, 1993; Cummings, 1994). Various PIC schemes to address the cancellation problem have been developed and there are interesting recent advances in this area (Chen & Parker, 2003; Mishchenko et al., 2004; Hatzky et al., 2007; Mishchenko et al., 2014; Startsev & Lee, 2014; Bao et al., 2018). Meanwhile, some continuum core codes avoided the cancellation problem completely (Rewoldt et al., 1987; Kotschenreuther et al., 1995), while others had to address somewhat minor issues resulting from it (Jenko, 2000; Candy & Waltz, 2003). Meanwhile, some continuum core codes avoided the cancellation problem completely (Rewoldt et al., 1987; Kotschenreuther et al., 1995), while others had to address somewhat minor issues resulting from it (Jenko, 2000; Candy & Waltz, 2003). With respect to the cancellation problem, one possible reason for the differences might be that in continuum codes the fields and particles are discretized on the same grid, whereas in PIC codes the particle positions do not coincide with the field grid. Because particle positions are randomly located relative to the field grid, one might need to be more careful in some way when treating the interaction of the particles and electromagnetic fields.
To this point, all published nonlinear electromagnetic gyrokinetic results have focused on the core region, mostly within the formulation neglecting the nonlinearity, although the ORB5 PIC code includes the nonlinearity and is effectively full- (Lanti et al., 2019). The XGC1 code is also full- and is focused on both the core and the edge/SOL; it has an option for a gyrokinetic ion/drift-fluid massless electron hybrid model (Hager et al., 2017), with a fully kinetic implicit electromagnetic scheme based on Chen & Chacon (2015) recently implemented and under further development (Ku et al., 2018b). Other gyrokinetic codes working on the SOL are not yet electromagnetic. Thus, to our knowledge, the results presented here are the first nonlinear electromagnetic full- gyrokinetic turbulence simulations on open field lines.
In this paper we present a numerical scheme for simulating the full- electromagnetic gyrokinetic system using a continuum approach. We use an energy-conserving discontinuous Galerkin (DG) scheme for the discretization of the gyrokinetic system in phase space, building on the work of Hakim et al. (2019); Shi et al. (2017); Shi (2017); Liu & Shu (2000). DG methods are attractive because they are highly local (enabling fairly straightforward parallelization schemes), allow high-order accuracy, and enforce local conservation laws (Durran, 2010). The present target of the scheme is simulating the edge and SOL of tokamaks, although the scheme could in principle be used for whole-device modelling including the core. Our scheme has been implemented as part of the gyrokinetics solver (Shi et al., 2017, 2019; Bernard et al., 2019) of the Gkeyll computational plasma framework, which also includes solvers for the Vlasov–Maxwell system (Cagas et al., 2017; Juno et al., 2018) and multi-moment fluid equations (Wang et al., 2015).
The paper is organized as follows. In Section 2, we describe the electromagnetic gyrokinetic system and some of its conservation properties. Section 3 describes the discontinuous Galerkin phase-space discretization of the system, and also presents proofs that the scheme preserves particle and energy conservation. The time-discretization scheme is handled in Section 4. In Section 5 we present some linear electromagnetic benchmarks that validate the scheme and also demonstrate the avoidance of the cancellation problem. We present nonlinear results showing the first electromagnetic gyrokinetic turbulence simulation on open field lines in Section 6, along with comparisons to a corresponding electrostatic simulation. We summarize and address future work in Section 7.
2 The electromagnetic gyrokinetic system
2.1 Basic equations
We solve the full- electromagnetic gyrokinetic (EMGK) equation in the symplectic formulation (Brizard & Hahm, 2007), which describes the evolution of the gyrocentre distribution function f_{s}(\mbox{\boldmathZ},t)=f_{s}(\mbox{\boldmathR},v_{\parallel},\mu,t) for each species , where is a phase-space coordinate composed of the guiding centre position \mbox{\boldmathR}=(x,y,z), the parallel velocity and the magnetic moment . In terms of the gyrocentre Hamiltonian and the Poisson bracket in gyrocentre coordinates, and also including collisions and sources (which do not derive from the bracket), the gyrokinetic equation is given by111One can use extended gyrocentre phase-space coordinates, which include time and the canonically conjugate energy , to include the time derivative terms in Eq. (1) inside an extended Poisson bracket (Brizard & Hahm, 2007). For ease of presentation we do not take this approach.
[TABLE]
or equivalently,
[TABLE]
where the gyrokinetic Poisson bracket is given by
[TABLE]
and we take the gyrocentre Hamiltonian to be
[TABLE]
Here we have taken the long-wavelength (drift-kinetic) limit to neglect gyroaveraging of the electrostatic potential , and we have also dropped higher- order terms in the Hamiltonian that appear in e.g. Brizard & Hahm (2007); extensions to include gyroaveraging will be included in later work, but these additions will not change the overall scheme presented here. The nonlinear phase-space characteristics are given by
[TABLE]
Here, B_{\parallel}^{*}=\mathbf{\hat{b}}\mbox{\boldmath\cdot}\mbox{\boldmathB^{}} is the parallel component of the effective magnetic field \mbox{\boldmathB^{}}=\mbox{\boldmathB}+(m_{s}v_{\parallel}/q_{s})\nabla\times\mathbf{\hat{b}}+\delta\mbox{\boldmathB}, where \mbox{\boldmathB}=B\mathbf{\hat{b}} is the equilibrium magnetic field and \delta\mbox{\boldmathB}=\nabla\times(A_{\parallel}\mathbf{\hat{b}})\approx\nabla A_{\parallel}\times\mathbf{\hat{b}} is the perturbed magnetic field (assuming that the equilibrium magnetic field varies on spatial scales longer than perturbations so that can be neglected). We neglect higher-order parallel compressional fluctuations of the magnetic field, so that \delta\mbox{\boldmathB}=\delta\mbox{\boldmathB}_{\perp}. The species charge and mass are and , respectively. In Eq. (6), note that we have separated into a term that comes from the Hamiltonian, , and another term proportional to the inductive component of the parallel electric field, . We use this notation for convenience, and so that the time derivative of the parallel vector potential appears explicitly. Further, we will assume a field-aligned coordinate system (e.g. Beer et al., 1995), and we will take the perpendicular directions to be and , and the parallel direction to be .
In the absence of collisions and sources , Eq. (1) can be recognized as a Liouville equation, which shows that the distribution function is conserved along the nonlinear characteristics. Liouville’s theorem also shows that phase-space volume is conserved,
[TABLE]
where is the Jacobian of the gyrocentre coordinates, and we will make the approximation \mathbf{\hat{b}}\mbox{\boldmath\cdot}\nabla\times\mathbf{\hat{b}}\approx 0 so that .
We can now write the gyrokinetic equation in conservative form,
[TABLE]
Here, we have used the symplectic formulation of electromagnetic gyrokinetics, where the parallel velocity is used as an independent variable (as opposed to the Hamiltonian formulation which uses the parallel canonical momentum as an independent variable) (Brizard & Hahm, 2007; Hahm et al., 1988). Notably, in the symplectic formulation, the time derivative of appears explicitly in the gyrokinetic equation, Eq. (8), and appears in but not in the Hamiltonian.
The electrostatic potential is determined by the quasi-neutrality condition in the long-wavelength limit, given by
[TABLE]
with the guiding centre charge density (neglecting gyroaveraging in the long-wavelength limit)
[TABLE]
Here, we have defined d\mbox{\boldmathw}=2\pi\,m_{s}^{-1}dv_{\parallel}\,d\mu=m_{s}^{-1}dv_{\parallel}\,d\mu\int d\alpha as the gyrocentre velocity-space volume element with the gyroangle integrated away and the Jacobian factored out. The polarization vector is
[TABLE]
where \nabla_{\perp}=\nabla-\mathbf{\hat{b}}(\mathbf{\hat{b}}\mbox{\boldmath\cdot}\nabla) is the gradient perpendicular to the background magnetic field. We use a linearized polarization density that we take to be a constant in time, which is consistent with neglecting a second-order energy term in the Hamiltonian. While the validity of this approximation in the SOL can be questioned due to large density fluctuations, a linearized polarization density is commonly used for computational efficiency (Ku et al., 2018a; Shi et al., 2019). Future work will include the nonlinear polarization density along with the second-order energy term in the Hamiltonian. The quasi-neutrality condition can then be rewritten as the long-wavelength gyrokinetic Poisson equation,
[TABLE]
Even in the long-wavelength limit with no gyroaveraging, the first-order polarization charge density on the left-hand side of Eq. (12) incorporates some finite Larmor radius (FLR) effects.
The parallel vector potential is determined by the parallel Ampère equation,
[TABLE]
Note that we can also take the time derivative of this equation to get a generalized Ohm’s law which can be solved directly for , the inductive component of the parallel electric field (Reynders, 1993; Cummings, 1994; Chen & Parker, 2001)
[TABLE]
Writing the gyrokinetic equation as
[TABLE]
where denotes all the terms in the gyrokinetic equation (including sources and collisions) except the term, Ohm’s law can be rewritten (after an integration by parts) as
[TABLE]
As we will show in Section 4, this form allows for the use of an explicit time-stepping scheme in which one can first compute (which does not involve , then compute , and finally compute . Note, however, that in some PIC approaches (Reynders, 1993; Chen & Parker, 2001), one must expand the right-hand side of Eq. (16) by inserting the gyrokinetic equation so that the right-hand side involves only moments of without time derivatives. In our continuum scheme we can compute directly and then perform the integration. Further, note that, although we are using the symplectic () formulation of EMGK, our Ohm’s law from Eq. (16) contains two integral terms which must cancel exactly. This is the root of the cancellation problem that appears in Ampère’s law in the Hamiltonian () formulation, and in Appendix A we show that the same cancellation problem could arise from Eq. (16) if the integrals are not treated consistently.
To model the effect of collisions we use a conservative Lenard–Bernstein (or Dougherty) collision operator (Lenard & Bernstein, 1958; Dougherty, 1964),
[TABLE]
where
[TABLE]
with n=\int d\mbox{\boldmathw}\ \mathcal{J}f. This collision operator contains the effect of drag and pitch-angle scattering, and it conserves number, momentum and energy density. Consistent with our present long-wavelength treatment of the gyrokinetic system, finite Larmor radius effects are ignored. For simplicity we restrict ourselves to the case in which the collision frequency is velocity independent, i.e. . Further details about this collision operator, including its conservation properties and its discretization, are left to a separate paper (Francisquez et al., 2020). In this work, we include only the effects of like-species collisions, which neglects electron–ion collisions and the resulting resistivity. A conservative scheme for cross-species collisions has also been implemented and will be included in later work. Extensions to a more complete collision operator are in progress.
2.2 Conservation properties
In the absence of collisions and sources, the Hamiltonian structure of the gyrokinetic system guarantees conservation of arbitrary functions of along the characteristics,
[TABLE]
along with corresponding Casimir invariants \int d\mbox{\boldmathR}\ d\mbox{\boldmathw}\ \mathcal{J}G(f), where d\mbox{\boldmathR}=dx\,dy\,dz. Thus, the system has an infinite number of conserved quantities such as the total particle number (or norm) N=\int d\mbox{\boldmathR}\ d\mbox{\boldmathw}\ \mathcal{J}f, the norm M=\int d\mbox{\boldmathR}\ d\mbox{\boldmathw}\ \mathcal{J}f^{2} and the kinetic entropy S=-\int d\mbox{\boldmathR}\ d\mbox{\boldmathw}\ \mathcal{J}f\ln f (Idomura et al., 2008).
The system also conserves total energy, , where the kinetic particle energy (neglecting the kinetic energy of the flow) is
[TABLE]
the (non-vacuum) electrostatic field energy (equivalent to the kinetic energy associated with the flow of particles) is
[TABLE]
the (perturbed) electromagnetic field energy is
[TABLE]
and
[TABLE]
(Note that is the sum of the particle kinetic energy and twice the potential energy, because every pair of particle interactions is double counted in the raw integral of .)
Assuming the boundary conditions are periodic or that the distribution function vanishes at the boundary so that surface terms vanish, the evolution of these quantities can be calculated as
[TABLE]
so that total energy is indeed conserved:
[TABLE]
3 The discrete EMGK system
In this section we describe the phase-space discretization of the electromagnetic gyrokinetic system used in Gkeyll.
3.1 Discrete equations
We use an energy-conserving discontinuous Galerkin scheme to discretize the gyrokinetic system in phase space. The scheme generalizes the algorithm of Liu & Shu (2000) (originally for the two-dimensional incompressible Euler and Navier-Stokes equations) to arbitrary Hamiltonian systems (Hakim et al., 2019; Shi et al., 2017; Shi, 2017). However, unlike the nodal approach used in Shi et al. (2017); Shi (2017), we use a modal DG scheme.
We start by decomposing the global phase-space domain into a structured phase-space mesh with cells . We then introduce a piecewise-polynomial approximation space for the distribution function f(\mbox{\boldmathR},v_{\parallel},\mu),
[TABLE]
where \mbox{\boldmathP}^{p} is some space of polynomials with maximum degree (by some measure). That is, are polynomial functions of in each cell, and \mbox{\boldmathP}^{p} is the space of the linear combination of some set of multi-variate polynomials. In this work, we choose \mbox{\boldmathP}^{p} to be an orthonormalized serendipity polynomial element space (Arnold & Awanou, 2011). The serendipity basis set has the advantage of using fewer basis functions while giving the same formal convergence order (although it is less accurate) as the Lagrange tensor basis, although note that, for , the serendipity basis is equivalent to the Lagrange tensor basis. We can then obtain the discrete weak form of the gyrokinetic equation by multiplying Eq. (8) by any test function and integrating (by parts) in each cell
[TABLE]
Solving this equation for all test functions in all cells yields the discretized distribution function , where the subscript denotes a discrete quantity in . In the surface terms, d\mbox{\boldmaths}_{R} is the differential element on a configuration-space surface (pointing outward normal to the surface), ds_{w}=2\pi\,m_{s}^{-1}d\mu\,\mbox{\boldmathn}\mbox{\boldmath\cdot}(\partial{\mbox{\boldmathZ}}/\partial{v_{\parallel}}) is the differential element on a surface, and the notation () indicates that the function is evaluated just inside (outside) the surface . The notation indicates a ‘numerical flux’, which takes a single value at the cell surface and in general can depend on the solution on both sides of the surface since the solution is discontinuous at the surface. Here, we choose to use standard upwind fluxes, which depend on the local value of the phase-space characteristic flow normal to the surface evaluated at each Gaussian quadrature point on the surface. Denoting the flow as \mbox{\boldmath\alpha}_{h}, the upwind flux can be expressed as
[TABLE]
where \mbox{\boldmathn}=d\mbox{\boldmaths}/|d\mbox{\boldmaths}| is the unit normal pointing out of the surface.
We will introduce a subset of where the piecewise polynomials are continuous across cell interfaces, denoted by \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{V}\cr} }}}^{p}_{h}. As we will show later, in order to preserve energy conservation in our discrete scheme, we will require that the discrete Hamiltonian be continuous across cell interfaces, i.e. H_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{V}\cr} }}}^{p}_{h} (Hakim et al., 2019; Liu & Shu, 2000; Shi et al., 2017; Shi, 2017). Note that one can show that this ensures that the discrete phase-space characteristics, \dot{\mbox{\boldmathR}}_{h}=\{\mbox{\boldmathR},H_{h}\} and , are also continuous across cell interfaces.222In a general non-orthogonal field-aligned geometry this is not necessarily true. This is because \mbox{\boldmathB}^{*}\mbox{\boldmath\cdot}\nabla z contains , which can be discontinuous in the direction. This makes the characteristic speed \dot{\mbox{\boldmathR}}_{h}\mbox{\boldmath\cdot}\nabla z discontinuous across cell interfaces. This will be addressed in a separate paper dealing with non-orthogonal field-aligned geometry.
We must also discretize the field equations. We introduce the restriction of the phase-space mesh to configuration space, , and we denote the configuration-space cells by for , where is the number of configuration-space cells. We also restrict to configuration space as
[TABLE]
Further, we introduce the subset of polynomials that are piecewise continuous across configuration-space cell interfaces \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}_{h}^{p}\subset\mathcal{X}_{h}^{p}, along with an additional subset \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}_{h}^{p}\subset\mathcal{X}_{h}^{p} where continuity is required in the directions perpendicular to the magnetic field, but not in the direction parallel to the field. Assuming a field-aligned coordinate system (e.g. Beer et al., 1995), we will take the perpendicular directions to be and , and the parallel direction to be .
Since we require to be continuous across all cell interfaces, this means that we require to be continuous, i.e. \phi_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}_{h}^{p}. Thus to solve the Poisson equation we use the (continuous) finite-element method (FEM). While one could ensure is continuous in all directions by using a three-dimensional FEM solve, we instead use a two-dimensional FEM solve in the and directions, followed by a one-dimensional smoothing operation in the direction. That is, we first solve for \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} using a two-dimensional FEM solve, and then we use a smoothing/projection operation to ensure continuity in the direction. We will denote this operation as \phi_{h}=\mathcal{P}_{z}[\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}] and define it below. We can make this splitting because only produces coupling in the and (perpendicular) directions.
For the two-dimensional solve, we solve for \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} by multiplying Eq. (12) by a test function \xi\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} and integrating (by parts) in each configuration-space cell to obtain the discrete local weak form
[TABLE]
where denotes the restriction of to cell and
[TABLE]
with the restriction of to velocity space. The global weak form is then obtained by summing Eq. (32) over cells in and (but not in ), which results in cancellation of the surface terms at cell interfaces and leaves only a global boundary term. Note that in order to maintain energetic consistency (as we will see below), the introduction of necessitates the modification of the right-hand side of Eq. (32) with , the adjoint of , defined as
[TABLE]
For the smoothing operation \phi_{h}=\mathcal{P}_{z}[\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}], we use a one-dimensional FEM solve in the direction. This can be written as the solution of the global (in ) weak equality
[TABLE]
where , with a subset of the configuration-space basis where continuity is required only in the direction. Here, denotes a restriction of the domain that is global in but cell-wise local in and . We remark that using an FEM solve for this operation makes self-adjoint, so that . Note, however, that one could instead use a different, local smoothing operation that is not self-adjoint, so we will keep the distinction between and . Also note that is a projection operator, in that \mathcal{P}_{z}[\mathcal{P}_{z}[\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}]]=\mathcal{P}_{z}[\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}].
The continuous discrete Hamiltonian H_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{V}\cr} }}}_{h}^{p} is then given by
[TABLE]
where is the projection of onto \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{V}\cr} }}}^{p}_{h}. Note that this is only necessary when is not in the basis, i.e. when , where is the maximum degree of the monomials in the basis set.
For the parallel Ampère equation we will take A_{\parallel h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} so that is continuous in and but discontinuous in . Multiplying Eq. (13) by a test function \varphi\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} and integrating, we can obtain the discrete weak form of this equation. The local weak form in cell is
[TABLE]
where again the surface terms will cancel on summing over cells except at the global boundary, and
[TABLE]
Here, note that we have replaced the in the definition from Eq. (13) with ; this will be required for energy conservation in the case, since when is not in the basis. Instead, for , , the piecewise-constant projection of . As before, we solve Eq. (37) using a two-dimensional FEM solve in the and directions. Note, however, that we do not require the smoothing operation in here because is allowed to be discontinuous in the direction, since it does not appear in the Hamiltonian in the symplectic formulation of EMGK.
The discrete weak form of Ohm’s law can be obtained by taking the time derivative of Eq. (37); after some manipulation, which we leave to Appendix B, the local weak form becomes
[TABLE]
where \partial A_{\parallel h}/\partial t\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h}, and
[TABLE]
so that the gyrokinetic equation can be written as
[TABLE]
Note that some special attention is required to ensure that upwinding of the numerical fluxes is handled consistently in Eqs. (39) and (42) in the case. The upwind flow for the surface terms is ; this is somewhat problematic because we cannot readily solve for from Eq. (39) without first knowing the upwind direction, which depends on . Thus for only, we use an approximate , calculated using Eq. (40) (which contains no surface term contributions), to compute the upwind direction for the surface terms in Eqs. (39) and (42). (One could extend this algorithm by iterating with a new estimate of the upwind direction based on the previous estimate of , but we leave that for future work. The present algorithm seems to work well for the cases tested so far, and we expect that results in the correct upwind direction most of the time.)
In our modal DG scheme, integrals in the above weak forms are computed analytically using a quadrature-free scheme that results in exact integrations (of the discrete integrands). (This means there are no aliasing errors, and that integration by parts operations that led to these integrals are treated exactly, for the specified discrete representation of and other factors in the integrand.) This is important for ensuring the conservation properties of the scheme, since the conservation laws in the EMGK system are indirect, involving integrals of the gyrokinetic equation. The fact that integrations are exact also has important implications for the cancellation problem. Since integrals in the discrete Ohm’s law are computed exactly, the discretization errors (which are solely embedded in the discrete integrands) cancel exactly, avoiding the cancellation problem. For more details about the modal scheme, the analytical integrations and the avoidance of the cancellation problem, we have included in Appendix C a derivation of a semi-discrete Alfvén wave dispersion relation that results from our scheme.
3.2 Discrete conservation properties
Now we would like to show that the discrete system (in the continuous-time limit) preserves various conservation laws of the continuous system. As with the continuous system, we will consider the conservation properties in the absence of collisions, sources and sinks, and we will assume that the boundary conditions are either periodic or that the distribution function vanishes at the boundary.
Proposition 1
The discrete system conserves total number of particles (the norm).
Proof 3.1**.**
Taking in the discrete weak form of the gyrokinetic equation, Eq. (29), and summing over all cells, we have
[TABLE]
where the surface terms cancel exactly at cell interfaces because the integrands (both the phase-space characteristics and the numerical fluxes) are continuous across the interfaces.
Proposition 3.2**.**
The discrete system conserves a discrete total energy, , where
[TABLE]
and
[TABLE]
Proof 3.3**.**
The proof follows from Proposition 3.2 in Hakim et al. (2019). We start by calculating
[TABLE]
The first term can be calculated by taking in Eq. (29) and summing over cells and species, since and H_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\displaystyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\textstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptstyle\mathcal{V}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\leaders\hrule\hfill\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{V}\cr} }}}_{h}^{p}\subset\mathcal{V}_{h}^{p}:
[TABLE]
Here, we see why we must require to be continuous; we want the surface terms to vanish, which means the integrands must be continuous across cell interfaces so that the contributions from either side of the interface cancel exactly when we sum over cells. The numerical flux is by definition continuous across the interface, and we have already noted above that the phase-space characteristics \dot{\mbox{\boldmathR}}_{h} and are also continuous across cell interfaces. This leaves the Hamiltonian, which we require to be continuous so that the surface terms do indeed vanish. Further, the first volume term vanishes exactly because \dot{\mbox{\boldmathR}}_{h}\mbox{\boldmath\cdot}\nabla H_{h}+\dot{v}^{H}_{\parallel h}\partial H_{h}/\partial v_{\parallel}=\{H_{h},H_{h}\}=0 by definition of the Poisson bracket. This leaves
[TABLE]
where here we see why we have defined using the derivative of instead of , as noted after Eq. (38). We now have the desired result for this term. For the second term in Eq. (55), we have
[TABLE]
Thus we have
[TABLE]
which is consistent with Eq. (24).
Next, we calculate
[TABLE]
where we have used \xi^{(j)}=\partial\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}\mskip 0.02998mu_{h}/\partial t in Eq. (32) to make the second equality, noting that the surface term vanishes upon summing over cells because \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} is continuous in the perpendicular directions. Here, we see why we modified the right-hand side of Eq. (32) with , so that the resulting term in Eq. (108) matches the one in Eq. (74).
Finally, we calculate
[TABLE]
where we have used in Eq. (37) to make the second equality, again noting that the surface term vanishes upon summing over cells because \partial A_{\parallel h}/\partial t\in\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\mathcal{X}\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\mathcal{X}\cr} }}}^{p}_{h} is continuous in the perpendicular directions.
We now have conservation of discrete total energy:
[TABLE]
We note that this proof did not rely on the particular choice of numerical flux function.
Proposition 3.4**.**
The discrete system exactly conserves the norm of the distribution function when using a central flux, while the distribution function norm monotonically decays when using an upwind flux.
Proof 3.5**.**
The proof is given as Proposition 3.3 in Hakim et al. (2019).
Proposition 3.6**.**
If the discrete distribution function remains positive definite, then the discrete scheme grows the discrete entropy monotonically,
[TABLE]
Proof 3.7**.**
The proof is given as Proposition 3.4 in Hakim et al. (2019).
4 Time-discretization scheme
So far we have considered only the discretization of the phase space for the system, and we have considered the conservation properties of the scheme in the continuous-time limit. Indeed, in the discrete-time system the conservation properties are no longer exact due to truncation error in the non-reversible time-stepping methods that we consider. However the errors will be independent of the phase-space discretization, and errors can be reduced by taking a smaller time step or by using a high-order time-stepping scheme to improve convergence. Following the approach of the Runge-Kutta discontinuous Galerkin method (Cockburn & Shu, 1998, 2001; Shu, 2009), we have implemented several explicit multi-stage strong stability-preserving Runge-Kutta high-order schemes (Gottlieb et al., 2001; Shu, 2002); the results in this paper use a three-stage, third-order scheme (SSP-RK3). These schemes have the property that a high-order scheme can be composed of several forward-Euler stages. Thus we will detail our time-stepping scheme for a single forward-Euler stage, which can then be combined into a multi-stage high-order scheme. Note that although we present the time-discretization scheme in this section in terms of our DG phase-space discretization, the scheme could be generalized to any spatial discretization.
Given and at time , the steps of the forward-Euler scheme to advance to time are as follows:
Calculate \mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}^{n} using Eq. (32), and then \phi_{h}^{n}=\mathcal{P}_{z}[\mathop{\mathchoice{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\displaystyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\textstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptstyle\phi\cr} }}{\vbox{ \offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern 1.0pt\cr\scriptscriptstyle\phi\cr} }}}_{h}^{n}] using Eq. (35).
[TABLE] 2. 2.
Calculate the partial EMGK update using Eq. (41).
[TABLE] 3. 3.
Calculate from Eq. (40) [for , this is only a provisional value, which we will denote as ].
[TABLE] 4. 4.
() Use the provisional from step 3 to calculate the upwinding direction in the surface terms in Eq. (39), and then calculate .
[TABLE] 5. 5.
Calculate the full EMGK update, , using Eq. (42). For , the provisional from step 3 should again be used to calculate the upwinding direction in the surface terms for consistency.
[TABLE] 6. 6.
Advance and to time .
[TABLE]
Note that the parallel Ampère equation, Eq. (37), is only used to solve for the initial condition of . For all other times, Eq. (143) is used to advance . This prevents the system from being over-determined and ensures consistency between and .
5 Linear benchmarks
5.1 Kinetic Alfvén wave
As a first benchmark of our electromagnetic scheme, we consider the kinetic Alfvén wave. In a slab (straight background magnetic field) geometry, with stationary ions (assuming ), the gyrokinetic equation for electrons reduces to
[TABLE]
Taking a single Fourier mode with perpendicular wavenumber and parallel wavenumber , the field equations become
[TABLE]
After linearizing the gyrokinetic equation by assuming a uniform Maxwellian background with density and temperature , so that , the dispersion relation becomes
[TABLE]
where , with , is the electron thermal speed, is the ion sound gyroradius, and is the plasma dispersion function (Fried & Conte, 1961). In the limit the wave becomes the standard shear Alfvén wave from magnetohydrodynamics (MHD), which is an undamped wave with frequency , where is the Alfvén velocity. For larger values of , the mode is damped by kinetic effects.
In Figure 1, we show the real frequencies () and damping rates () obtained by solving Eq. (148) for a few values of . We also show numerical results from Gkeyll, which match the analytic results very well. These results are a good indication that our scheme avoids the Ampère cancellation problem, which can cause large errors for modes with (see Appendix A); we see no such errors, even for the case with . Each Gkeyll simulation was run using piecewise-linear basis functions () in a reduced dimensionality mode with one configuration space dimension and one velocity space dimension, with the number of cells in each dimension. The perpendicular dimensions ( and ), which appear only in the field equations in this simple system, were handled by replacing by , as in Eqs. (145) and (146). We use periodic boundary conditions in and zero-flux boundary conditions in .
We also show in Figure 2 the fields and for the case with and , which gives . For these parameters the system is near the MHD limit, which means we should expect . While this condition is never enforced, getting the physics correct requires the scheme to allow . The fact that our scheme allows discontinuities in in the parallel direction is an advantage in this case. Because is piecewise-linear here, is piecewise-constant; this is necessarily discontinuous for non-trivial solutions. Thus the scheme produces a piecewise-constant in this MHD-limit case, as shown in Figure 2, resulting in . If our scheme did not allow discontinuities in , a continuous would never be able to exactly cancel a discontinuous , and the resulting would make the solution inaccurate. Notably, this would be the case had we chosen the Hamiltonian () formulation of the gyrokinetic system, which uses as the parallel velocity coordinate. This is because is included in the Hamiltonian in the formulation, which would require continuity of (and thereby ) to conserve energy in our discretization scheme.
5.2 Kinetic ballooning mode (KBM)
We use the kinetic ballooning mode (KBM) instability in the local limit as a second linear benchmark of our electromagnetic scheme. The dispersion relation is given by solving (Kim et al., 1993)
[TABLE]
where
[TABLE]
with , , , and with the modified Bessel function. Here, the wavenumbers and are normalized to and , respectively, and the frequencies and are normalized to . In the local limit, and do not vary along the field line. Note that in Eq. (149) we have modified the FLR terms from (Kim et al., 1993) so that we can take while keeping to neglect all FLR effects except for the first-order polarization term, which is consistent with our long-wavelength Poisson equation.
The local limit can be achieved by simulating a helical flux tube with no magnetic shear, which gives a system with constant magnetic curvature that corresponds to . This geometry has been previously used for SOL turbulence studies with Gkeyll (Shi et al., 2019; Bernard et al., 2019), except in this section we take the boundary condition along the field lines to be periodic. We will provide further details about the helical geometry and the coordinates in the next section.
We show the results of Gkeyll simulations of the KBM instability in the local-limit helical geometry for several values of in Figure 3. The results agree well with the analytic result obtained by numerically solving Eqs. (149-150). The parameters are chosen to match those used in figure 1 of Kim et al. (1993), although the differences in FLR terms () cause our growth rates to be larger than those in Kim et al. (1993). These are fully five-dimensional simulations with the real deuterium–electron mass ratio using piecewise-linear () basis functions, with . The boundary conditions are periodic in the three configuration-space dimensions and zero flux in the velocity dimensions. The initial distribution function of each species is composed of a background Maxwellian with gradients in the density and temperature corresponding to the desired and , plus a perturbed Maxwellian (for the electrons only) with small sinusoidal variations in the density corresponding to the desired and . Note that since we are using a full- representation, the presence of a background gradient in the distribution function means that we must apply the periodic boundary conditions by first subtracting off the initial background distribution function, then applying periodicity to the perturbations only and then adding back the background distribution. Finally, we note that since Gkeyll is designed primarily for nonlinear calculations, the fact that Fourier modes are not eigenfunctions of the DG discretization of the system makes these linear tests somewhat difficult for Gkeyll. This may play a role in the small deviation of the results from the analytical theory. Because of this, Fourier modes other than the one initialized can grow and pollute the results. In particular, we have not included results from the ion temperature gradient (ITG) branch because we find that a mode with grows and overcomes the initialized finite mode before its growth rate has converged.
6 Nonlinear results
We now present preliminary nonlinear electromagnetic results from Gkeyll. We simulate turbulence on helical, open field lines as a rough model of the tokamak scrape-off layer. These simulations are a direct extension of the work of Shi et al. (2019) to include electromagnetic fluctuations. As such, we use the same simulation geometry and similar NSTX-like parameters. In the non-orthogonal field-aligned geometry, is the radial coordinate, is the coordinate along the field lines, and is the binormal coordinate which labels field lines at constant and . These coordinates map to physical cylindrical coordinates ( via , , . (Note that this fixes an error in the mapping in (Shi et al., 2019).) The field-line pitch is taken to be constant, with the vertical component of the magnetic field (analogous to the poloidal field in typical tokamak geometry), and the total magnitude of the background magnetic field. Further, is the radius of curvature at the centre of the simulation domain, with the device major radius and the minor radius. As in Shi et al. (2019), we neglect all geometrical factors arising from the non-orthogonal coordinate system, except for the assumption that perpendicular gradients are much stronger than parallel gradients so that we can approximate
[TABLE]
where we have used \mbox{\boldmathB}=B_{\text{axis}}(R_{0}/R)\mathbf{\hat{e}}_{z} in the last step, with the magnetic field strength at the magnetic axis.
We use this geometry to simulate a flux-tube-like domain on the outboard side that wraps helically around the torus and terminates on conducting plates at each end in . The simulation box is centreed at with dimensions cm, cm, and m, where m and . Periodic boundary conditions are used in the direction, and a Dirichlet boundary condition is applied in x, which effectively prevents flows into the boundaries in . Conducting-sheath boundary conditions are applied in the direction (Shi et al., 2017, 2019), which partially reflect one species (typically electrons) and fully absorbs the other species depending on the sign of the sheath potential. This involves solving the gyrokinetic Poisson equation to evaluate the potential at the boundary, corresponding to the sheath entrance, and using the resulting sheath potential to determine a cutoff velocity below which particles are reflected by the sheath. Notably, our sheath boundary condition allows current fluctuations in and out of the sheath, which we will find to important later in this section. This is different from the standard logical sheath boundary condition (Parker et al., 1993b) that imposes that there is no net current to the sheath by assuming that the ion and electron currents at the sheath entrance are equal at all times. The velocity-space grid has extents and , where and . We use piecewise-linear () basis functions, with . Note that although the domain that we simulate is a flux tube, the simulations are not performed in the local limit; the simulations include radial variation of the magnetic field and the profiles, and are thus effectively global.
The simulation parameters are similar to those used in Shi et al. (2019), roughly approximating an H-mode deuterium plasma in the NSTX SOL: T, m, m, eV. For the particle source, we use the same form as in Shi et al. (2019) but we increase the source particle rate by a factor of 10 to access a higher regime where electromagnetic effects will be more important. The source is localized in the region , with m and m. The location , which separates the source region from the SOL region, can be thought of as the separatrix. A floor of one tenth the peak particle source rate is used near the midplane to prevent regions of from developing at large . Unlike in Shi et al. (2019) we do not use numerical heating to keep despite the fact that our DG algorithm does not guarantee positivity. While the simulations appear to be robust to negative in some isolated regions, lowering the source floor in the SOL region leads to simulation failures due to positivity issues at large . A more sophisticated algorithm for ensuring positivity is left to future work. We also artificially lower the collision frequency to one tenth the physical value to offset the increased particle source rate so that the time-step limit from collisions does not become too restrictive. Further, we model only ion–ion and electron–electron collisions, leaving cross-species collisions to future work.
Using the novel electromagnetic scheme described in this paper, we ran a simulation in this configuration to ms, with a quasi-steady state being reached around when the sources balance losses to the end plates. For reference, the ion transit time is . In Figure 4 we show snapshots of the density, temperature and of electrons (top row) and ions (bottom row). The snapshots are taken at the midplane () at s. We can see a blob with a mushroom structure being ejected from the source region. We also show in Figure 6 snapshots of the electromagnetic fields taken at the same time and location. We show the electrostatic potential , the parallel magnetic vector potential and the normalized magnetic fluctuation amplitude (top row), along with the components of the parallel electric field (bottom row). Note that only , and are evolved quantities in the simulation, with the other quantities derived. We see that is of comparable magnitude to , indicating that the dynamics is in the electromagnetic regime. Significant magnetic fluctuations of over can be seen in in this snapshot. We also show in Figure 6 the time and spatially averaged profile of magnetic fluctuations vs , which shows that, on average, we observe magnetic fluctuations of the order of . This radial profile, and similar ones that will follow, is computed by averaging in and using data near the midplane over the period of ms.
In Figures 8 and 8 we show projections of the three-dimensional magnetic field line trajectories. These plots are created by integrating the field line equations for the total (background plus fluctuation) magnetic field. In Figure 8, each field line starts at m and either m or m for a range of values and is traced to m. The starting points (at m ) are marked with circles, while the ending points (at m) are marked with crosses. The trajectories have been projected onto the plane, and we have also plotted the ion density at m in the background. From left to right, we show a short time series of snapshots, with and s. At s, a blob is starting to emerge from the source region at m. The field lines that start at m are beginning to be stretched radially outward as the blob emerges. In the s snapshot, we see that the blob is now propagating radially outward into the SOL region and the m field lines have been stretched further. The field lines that start at m are now also starting to be stretched near m, and they are stretched even more in the s snapshot as the blob continues to propagate. We can also see the remnants of another blob that was ejected near m in previous frames. In the s snapshot, the field lines have been stretched by this blob, but by s the field lines in this region have returned closer to their equilibrium position. This behaviour of blobs bending and stretching the field lines is an inherently full- phenomenon. The blobs have a higher density and temperature than the background, so they raise the local plasma as they propagate. This causes the field lines to move with the plasma, allowing the fields lines to be deformed and stretched by the radially propagating blobs and ultimately leading to larger magnetic fluctuations.
In Figure 8 we show a slightly different view of the field-line trajectories at s. Field lines are still traced from the bottom ( m) to the top ( m), but now each field line starts at m for a range of . The starting points are again marked with circles and the ending points are marked with crosses. We have projected the three-dimensional trajectories onto the plane in Figure 8, and onto the plane in Figure 8. In we again plot the ion density at m in the background; in the ion density has been averaged over m. As can be seen in Figure 8, the blob propagating near m has stretched several field lines radially outward near the midplane. These bowed-out field lines originate from a range of values, 1.3\m m, and have all been dragged along with the blob as it was ejected from the source region and propagated radially outward. We also see some degree of line tying in these plots, with many of the field lines ending at a similar point in space to where they began, despite being stretched near the midplane. The field lines are not perfectly line tied, however; if they were, the crosses would perfectly align with their corresponding circles in the projections. Because our sheath boundary condition allows current fluctuations at the sheath interface, we can model the finite resistance of the sheath, which makes line tying only partial (Kunkel & Guillory, 1966). This allows the footpoints of the field lines to slip at the sheath interface (Ryutov, 2006). Examining figures 8 and 8, we see evidence of this in the simulation, with most of the end points moving slowly and smoothly in the vicinity of their origin, especially at larger . In the source region, however, there are other field lines whose end points suddenly jump further away from their origin. This suggests that we are seeing line breaking (reconnection) due to electron inertia effects and numerical diffusion, as field lines are pushed close together by large perturbations in the source region.
We have also run a corresponding electrostatic simulation in this configuration for direct comparison. This simulation is identical in configuration to the m case from Shi et al. (2019) except for the increased particle source rate and lack of cross-species collisions. In Figure 10 we show a comparison of radial profiles of density, temperature and for the electromagnetic and electrostatic cases. The profiles for the electromagnetic case are shallower in the SOL region and steeper in the source region. This suggests that there is less radial particle and heat transport into the SOL region in the electromagnetic case. This is in part confirmed by the profile of the radial particle flux in Figure 10, showing approximately 40% less particle transport in the electromagnetic case. The total particle flux includes the particle flux, , with . In the electromagnetic case, also includes the particle flux due to magnetic flutter, , with . The tilde indicates the fluctuation of a time-varying quantity, defined as with the time average of . The brackets denote an average in , (near the midplane), and time. We also note that the electromagnetic profiles might be even shallower in the SOL region if not for the floor on the source used to prevent positivity issues in the distribution function at large .
We also compare fluctuation statistics between the electromagnetic and electrostatic cases in Figure 11. Statistics of the electron density are shown on the top row, the middle row shows statistics of electrostatic potential fluctuations and the bottom row shows statistics of the radial electron particle flux. Despite the fact that the electromagnetic case shows less particle transport, the root-mean-square relative density fluctuations are larger in the electromagnetic case by up to a factor of two. The electromagnetic case also has higher skewness and excess kurtosis, indicating that the density fluctuations in the electromagnetic case are more intermittent. The skewness and kurtosis of the particle flux also indicate more intermittency in the electromagnetic case. This contributes to the reduced transport shown in the electromagnetic case, since the transport events are rarer even if the fluctuation levels are larger. Meanwhile, the fluctuation statistics for the potential are relatively similar between the electromagnetic and electrostatic cases. The statistics of the density and potential appear more coupled in the electrostatic case, consistent with the electrons behaving adiabatically when electromagnetic effects are neglected. In the electromagnetic case the density fluctuations are more intermittent and higher amplitude than the potential fluctuations.
Finally, we note that in terms of computational cost, the electromagnetic simulation is less than twice as expensive as the corresponding electrostatic simulation on a per-time-step basis. On 128 cores, the time per time step was 0.41 s/step for the electrostatic simulation and 0.68 s/step for the electromagnetic simulation. The increased cost is due to the additional field solves required for Ohm’s law, along with additional terms in the gyrokinetic equation. However, due to time-step restrictions on an electrostatic simulation due to the electrostatic shear Alfvén mode (also known as the mode) (Lee, 1987), the electromagnetic simulation makes up some of the additional cost by taking slightly larger time steps. The total wall-clock time (on 128 cores) for the electrostatic simulation was approximately 65 h, and the electromagnetic simulation took about 82 h. Altogether, the cost of these simulations is relatively modest, and the addition of electromagnetic effects only makes the simulations marginally (%) more expensive. We also note that the new version of Gkeyll, which uses a quadrature-free modal DG scheme, is approximately 10 times faster than the previous version of Gkeyll used in Shi et al. (2019), which used nodal DG with Gaussian quadrature. More details about the improvements from the quadrature-free modal scheme will be reported elsewhere.
7 Summary & Conclusion
In this paper we have presented an energy-conserving scheme for the full- electromagnetic gyrokinetic system. We choose the symplectic formulation of EMGK, which uses the parallel velocity as an independent variable. This leads to the time derivative of the parallel vector potential, , appearing explicitly in the gyrokinetic equation. We handle this term directly by solving an Ohm’s law. We presented the conservation properties of the EMGK system.
We described the discontinuous Galerkin scheme used to discretize the EMGK system in phase space. We proved that the scheme preserves particle conservation, and that the scheme also preserves energy conservation provided that the discrete Hamiltonian is continuous. This is achieved by using the (continuous) finite-element method for the field solves. We also detailed a basic forward-Euler time-stepping scheme to be used in the stages of a multi-stage high-order SSP-RK scheme. The time-stepping scheme updates the gyrokinetic equation in two parts, with the result of the first part [denoted ] being used in Ohm’s law to solve directly for so that it can then be used in the second part of the gyrokinetic update.
We have implemented the scheme in the gyrokinetic solver of the Gkeyll computational plasma framework. We provided two linear benchmarks to validate the electromagnetic scheme: a kinetic Alfvén wave calculation and a local kinetic ballooning mode instability calculation. In both cases results from Gkeyll agree well with analytic results. The success of these calculations, especially in cases with , indicates that the scheme avoids the Ampère cancellation problem. Further, the discontinuous Galerkin nature of the scheme enables (but does not enforce) exact cancellation of the components of , allowing the system to capture the MHD limit with .
We presented a nonlinear electromagnetic full- gyrokinetic simulation of turbulence on helical, open field lines as a rough model of the tokamak scrape-off layer. This simulation is the first nonlinear electromagnetic gyrokinetic simulation on open field lines. We showed data illustrating the interplay between blobs propagating into the SOL and the resulting bending and stretching of magnetic field lines. We also made quantitative comparisons between the electromagnetic simulation and a corresponding electrostatic simulation. Notably, the electromagnetic simulation exhibits less transport and shallower density and temperature profiles in the SOL, as well as larger root-mean-square density fluctuations with more intermittency. Future work will examine the particle and energy balance in these nonlinear simulations to confirm the conservation properties proved in Section 3.2 after accounting for sources and losses to the walls.
A number of extensions have been left to future work. The capability to simulate more realistic magnetic geometries using a non-orthogonal field-aligned coordinate system is currently in progress, so that effects of magnetic shear and non-constant curvature can be included. The inclusion of both open and closed-field-line regions, including the X-point in diverted geometries, is an additional complication that must be addressed. Gyroaveraging is another important effect that must be implemented to improve fidelity. Further, a robust solution to the issue of maintaining positivity of the distribution function has been implemented for Hamiltonian terms and is in progress for the collision operator. This could, for example, alleviate the need to use a source floor in the nonlinear simulations presented in Section 6, which could further enhance the differences between the electromagnetic and electrostatic profiles.
The modest cost of the nonlinear full- gyrokinetic simulations that we have presented make the prospect of using Gkeyll for whole-device modelling a feasible goal. The inclusion of electromagnetic effects will be crucial to the fidelity of these efforts. Such a tool would be invaluable to future studies of turbulent transport in fusion devices, both from a theoretical perspective and also as a model for predicting and analysing the performance of current and future experiments.
Acknowledgements
We would like to thank Tess Bernard, Petr Cagas, James Juno and other members of the Gkeyll team for helpful discussions and support, including the development of the postgkyl post-processing tool which facilitated the creation of many figures in this paper. Research support came from the U.S. Department of Energy: NRM was supported by the DOE CSGF program, provided under grant DE-FG02-97ER25308; AH and GWH were supported by the Partnership for Multiscale Gyrokinetic Turbulence (MGK) Project, part of the DOE Scientific Discovery Through Advanced Computing (SciDAC) program, via DOE contract DE-AC02-09CH11466 for the Princeton Plasma Physics Laboratory; MF was supported by DOE contract DE-FC02-08ER54966. Computations were performed on the Eddy cluster at the Princeton Plasma Physics Laboratory.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/S0022377820000070.
Appendix A Ampère cancellation problem
To understand the root of the Ampère cancellation problem, we examine the simple Alfvén wave case from Section 5.1. Recall that in this simple case, the gyrokinetic system is given by
[TABLE]
Note that now we have inserted two constants, and , in the integrals in Eq. (155). We will use these constants to represent small errors that could arise in the numerical calculation of these integrals. As in Section 5.1, we can calculate the dispersion relation for this system, but now we will take the limit , so that the dispersion relation reduces to
[TABLE]
where recall that . This reduces to the correct result if . However, if , there will be large errors for modes with . This means the two integrals in Eq. (155) must be calculated carefully and consistently so that any errors exactly cancel.
A.1 Hamiltonian () formulation
We now briefly discuss the cancellation problem in the Hamiltonian () formulation of gyrokinetics. In this case, the simple system above becomes
[TABLE]
where we have again included constants and to represent numerical errors in the integrals. The resulting dispersion relation is identical to Eq. (156), so the two integrals must again be calculated carefully and consistently so that errors exactly cancel. This could be slightly more challenging than in the symplectic case. Suppose the grid does not extend to infinity but has some finite limits which are expected to be large enough in practice. These limits are used when numerically computing the integrals. Since depends on a time-dependent quantity in , it is possible that the limits may need to be time-dependent if fluctuations are large in order to consistently compute the integrals.
Appendix B The discrete weak form of Ohm’s law
To obtain the discrete weak form of Ohm’s law, we start by taking the time derivative of Eq. (37):
[TABLE]
Now, note that, analogously to Eq. (15), we can write the discrete weak form of the gyrokinetic equation as
[TABLE]
where
[TABLE]
Substituting in Eq. (161) and summing over velocity cells, we obtain
[TABLE]
Note that, for , the surface term on the right-hand side vanishes because is continuous across cell interfaces when is included in the basis, resulting in cancellations. However, for this term is not continuous, and we must keep this surface term; further, the last term on the right-hand side vanishes for since . We can now substitute this result into the right-hand side of Eq. (160), giving
[TABLE]
In Eq. (164), is the piecewise-constant projection of .
Appendix C Semi-discrete dispersion relation for Alfvén wave
In this appendix we derive a semi-discrete Alfvén wave dispersion relation using a piecewise-linear DG discretization for only the coordinate. The main purpose of this appendix is to show how our discrete scheme avoids the Ampère cancellation problem. We will also show how the integrals in the DG weak form are computed analytically in our modal scheme.
The semi-discrete gyrokinetic weak form for this system is
[TABLE]
We begin by mapping each cell to via the transformation , where is the cell centre of cell
[TABLE]
Taking an orthonormal piecewise-linear basis in , , we expand on the basis in cell as
[TABLE]
(Note that in the fully discretized case all coordinate dependence would be contained in multi-variate basis functions.) We can then analytically integrate the weak form for each to obtain the modal evolution equation for each DG ‘mode’ :
[TABLE]
Finally, we will make the ansatz and linearize:
[TABLE]
We now turn to the field equations. The Poisson equation is
[TABLE]
Expanding and using the ansatz, this becomes
[TABLE]
where we will define so that by definition. For Ohm’s law, we must use the form from Eq. (39), which gives
[TABLE]
where
[TABLE]
Again expanding and using the ansatz, Ohm’s law becomes
[TABLE]
Analogously to Appendix A, we can rewrite this equation as
[TABLE]
where we have defined
[TABLE]
Clearly , which allows us to move the term to the right-hand side, giving a term proportional to the total parallel electric field, :
[TABLE]
This is essential for avoiding the cancellation problem because if we instead had , we would have had a leftover term proportional to on the left-hand side. This leftover term would then lead to the spurious term proportional to in Eq. (156).
In order to compute the integral quantities in the field equations, we use Eq. (171) to compute
[TABLE]
where we have expanded in the limit . Now we can calculate
[TABLE]
Substituting these integral quantities into the field equations, the Poisson equation becomes
[TABLE]
and Ohm’s law becomes
[TABLE]
where we have substituted the definition of . We can now combine Eqs. (185) and (186) by multiplying Eq. (185) by , multiplying Eq. (186) by and summing the two equations to get
[TABLE]
with . This then yields the dispersion relation
[TABLE]
To evaluate and the other sum, we need to project the background onto the basis in each cell. Taking , we project onto the basis in cell as
[TABLE]
where we have taken the cell centre to be . Now we can evaluate integrated quantities such as
[TABLE]
where now note that we have finite limits on the sum to indicate finite extents of the grid. As we alluded to before, we will define so that by definition, which means
[TABLE]
Note that quickly approaches 1 with increasing , so that for example when , . We can also calculate
[TABLE]
where is the sign of the upwind velocity, and the boundary terms that result from the finite limits on the sum are small for . Thus we have as expected, although it does not need to be exactly equal to unity to eliminate the cancellation problem. Instead, it was sufficient that on either side of Eq. (178).
One can also show that
[TABLE]
Now substituting these results into the dispersion relation from Eq. (188), we obtain
[TABLE]
after again taking the limit and assuming . This finally gives
[TABLE]
which is the expected dispersion relation.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Arnold & Awanou (2011) Arnold, D. N. & Awanou, G. 2011 The serendipity family of finite elements. Found. Comput. Math. 11 (3), 337–344.
- 2Bao et al. (2018) Bao, J., Lin, Z. & Lu, Z. X. 2018 A conservative scheme for electromagnetic simulation of magnetized plasmas with kinetic electrons. Phys. Plasmas 25 (2), 022515.
- 3Beer et al. (1995) Beer, M. A., Cowley, S. C. & Hammett, G. W. 1995 Field–aligned coordinates for nonlinear simulations of tokamak turbulence. Phys. Plasmas 2 (7), 2687–2700.
- 4Bernard et al. (2019) Bernard, T. N., Shi, E. L., Gentle, K. W., Hakim, A., Hammett, G. W., Stoltzfus-Dueck, T. & Taylor, E. I. 2019 Gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak. Phys. Plasmas 26 (4), 042301.
- 5Brizard & Hahm (2007) Brizard, A. J. & Hahm, T. S. 2007 Foundations of nonlinear gyrokinetic theory. Rev. Mod. Phys. 79 , 421–468.
- 6Cagas et al. (2017) Cagas, P., Hakim, A., Juno, J. & Srinivasan, B. 2017 Continuum kinetic and multi–fluid simulations of classical sheaths. Phys. Plasmas 24 (2), 022118.
- 7Candy & Waltz (2003) Candy, J. & Waltz, R. E. 2003 An Eulerian gyrokinetic–Maxwell solver. J. Comput. Phys. 186 (2), 545–581.
- 8Chen & Chacon (2015) Chen, G. & Chacon, L. 2015 A multi-dimensional, energy- and charge-conserving, nonlinearly implicit, electromagnetic Vlasov–Darwin particle-in-cell algorithm. Comp. Phys. Comm. 197 , 73–87.