Energy conservation tests of a coupled kinetic plasma–kinetic neutral transport code

The design of the ITER¹ magnetic fusion experiment is based on extrapolations of empirical scaling relations for the plasma confinement time developed from present-day devices [1]. Since the fusion gain (ratio of fusion to heating power) achieved by ITER is very sensitive to deviations from those scalings, the magnetic confinement fusion community has made a concerted effort to understand the physics underlying the scalings so as to increase our confidence in their extrapolation to the dimensionless parameters anticipated for ITER.

The principal plasma confinement mode targeted for ITER, the so-called 'high confinement mode' or 'H-mode', is characterized by a substantial barrier to outward radial plasma transport just inside the last closed magnetic flux surface (the magnetic separatrix in a diverted configuration like ITER). The resulting plasma density and temperature profiles possess extremely steep radial gradients just inside the last closed surface, but are much flatter between there and the magnetic axis. Consequently, this steep gradient region is referred to as the 'H-mode pedestal'.

The associated steep gradients in plasma energy and current are potential sources of free energy for magnetohydrodynamic modes of various scales. The most common instabilities in this region, called 'edge localized modes' or ELMs, act to periodically flatten the pedestal gradients. As long as the plasma remains in the H-mode, the pedestal builds up again rapidly and the cycle repeats.

The effort to understand the H-mode has proceeded on multiple fronts, attempting to answer key questions:

• (i)  
  What physics controls the transition from 'low confinement' or 'L-mode' into the H-mode? In particular, what determines the observed threshold heating power for this transition?
• (ii)  
  What instabilities drive the wide variety of observed ELMs and how can we control them?
• (iii)  
  What transport processes govern the pedestal buildup and the subsequent near-equilibrium phase of the H-mode?

The research effort and codes described in this paper are targeted at the third of these questions. Our understanding of plasma transport in magnetic confinement devices is poorest with regard to the anomalous radial transport associated with plasma microturbulence. However, in the H-mode pedestal this microturbulence is much less virulent than in the L-mode, and radial transport fluxes are reduced to levels consistent with 'neoclassical' theory. Neoclassical transport sets an upper limit on our ability to confine an assembly of plasma particles orbiting in a toroidal magnetic field and undergoing self-collisions. These effects result in unequal losses of electrons and ions, giving rise to plasma currents and electric fields. Because neoclassical transport theory is a relatively mature topic, an examination of its implications for the H-mode pedestal is an obvious target for making progress toward a greater understanding.

Typical H-mode pedestal parameters, however, violate the assumptions of most analytic neoclassical theories: the pedestal gradient scale lengths can be narrower than the neoclassical ('banana') ion orbit widths, the plasma ion distribution is expected to deviate significantly from a Maxwellian, and the role of the magnetic separatrix, i.e. magnetic topology, is significant. For these and other reasons that have become apparent over the last several years, accurately simulating neoclassical behavior in the H-mode pedestal requires a kinetic, particle-based simulation technique. This is the approach taken here, which builds on the initial study by Chang et al [2] examining the effects of a self-consistent radial electric field, collisions and neutral particles on the pedestal characteristics. The need for a kinetic treatment was confirmed by the observation of significant non-Maxwellian ion distributions. Recycled neutral atoms were found to result in a buildup of the pedestal density, provided anomalous radial transport was indeed less than that due to neoclassical effects.

We describe here improvements to the neutral transport routine in the code developed by Chang (XGC0) that yield a more realistic model of plasma–neutral interactions and the ability to incorporate plasma–wall interactions. First, we resolve plasma and neutral quantities poloidally, as well as radially, throughout the entire vacuum vessel volume, allowing the neutral sources to be placed at material surfaces. The spatial variation of those sources is determined directly from the distribution of ions lost to the walls in a physically realistic way. To allow the resulting neutral profile to evolve consistently, the neutral transport is simulated in a time-dependent manner and synchronized with the plasma time evolution. Both plasma and neutral species have three-dimensional flow velocity vectors so that the effects of neutral particles on the transport of plasma momentum are treated. The atomic physics interactions between plasma and neutral particles are characterized in detail, including the density dependence of the electron impact ionization of hydrogen and the velocity dependence of the charge exchange cross section. Molecules and impurities will be incorporated in a subsequent version of the code.

The algorithm used for these calculations is described in section 2. The results of a standard test case, comparable to that used in [2], are presented in section 3. The conservation properties of the code in these example runs are examined in section 4, documenting an essential aspect of the verification of this algorithm. Finally, some conclusions and plans for future improvements are described in section 5.

The treatment of ion–ion collisions in the XGC0 code [2] is complicated by the nonlinearity of the collision operator. In fact, the development of techniques for accurately and efficiently simulating these collisions remains an active area of research. The collision operator governing plasma–neutral interactions is likewise effectively nonlinear due to the tight coupling between plasma ions and atoms in most regions of the tokamak boundary.

Our approach for dealing with this nonlinearity is closely related to the 'test particle Monte Carlo' method attributed to Haviland and Lavin [3]. In Haviland's original technique, and in that of subsequent practitioners (e.g. [4–6]), individual test particles collide with a background characterized by a specified distribution function. In the purely nonlinear case, that distribution is typically updated in an iterative fashion from the test particle trajectories until a converged (for steady-state problems) result is obtained. The accuracy of this method hinges on the adequacy of the chosen representation of the background distribution function, e.g. a drifting Maxwellian.

The next level of sophistication, direct simulation Monte Carlo [7], bypasses this limitation by colliding the current test particle with other test particles present in the same computational cell. Because this effectively represents a Monte Carlo discretization of six dimensions of velocity space, versus three for test particle Monte Carlo, the computational resources required to achieve a specified precision in the results are significantly greater. While we do not consider this approach in this work, the continued and rapid increase in available computational speed will eventually render it practical.

A significant difference between our problem and those addressed in the aforementioned references is that the transport of the plasma (i.e. charged) and neutral species is fundamentally different in character, and the two phenomena have frequently been simulated with separate codes. The most familiar examples to magnetic fusion researchers are coupled fluid plasma and Monte Carlo neutral transport, e.g. B2-EIRENE [8, 9] and UEDGE–DEGAS2 [10]. The coupling between such codes has largely been explicit and iterative or time-dependent. The two-dimensional plasma density, flow velocity and temperature fields produced by the plasma code are input as a background to the neutral transport code; the plasma fluxes to the material boundary provide the source of neutral particles via plasma–material interactions, referred to as 'recycling' [11]. Execution of the neutral transport code yields volumetric sources (or sinks) of plasma mass, momentum and energy that are transferred back to the plasma code for use as the 'right-hand sides' of its transport equations in its next iteration or time step advance. Adapting this method to the present application is problematic since all of the kinetic details of the plasma–neutral interactions simulated by the neutral routine are lost in the process of compiling the volumetric sources.

We instead follow Chang et al's [2] approach of extending the test particle Monte Carlo method and employ two complementary plasma–neutral operators, one embedded in the plasma code and one in the neutral transport routine. The neutral collision operator works in a manner very similar to that used in coupling to fluid plasma codes: the kinetic neutral species collide with a plasma background characterized by a particular distribution function, the parameters of which are specified by the plasma code. Instead of returning volumetric sources, the neutral routine provides the parameters needed to specify the background neutral distribution in the plasma code's collision operator.

The accuracy and realism of the overall simulation hinges on the two operators being consistent with each other and the representation of the distributions being appropriate to the problem at hand. In this work, the plasma transport calculation is performed by the XGC0 code originally described by Chang et al [2] and subsequently extended by others [12]; the neutral transport routine is based on the DEGAS2 Monte Carlo neutral transport code [13]. Henceforth, we will refer to the implementations of these two complementary plasma–neutral collision operators as the 'XGC0' (kinetic plasma particles colliding on a neutral background) and 'DEGAS2' (kinetic neutral particles colliding on a plasma background) collision routines.

The XGC0 code reduces the ion phase space from six dimensions to five by simulating only the motion of the guiding centers of the ion orbits, averaging over the much faster gyromotion. That is, the ion velocity is written as

$$\begin{equation} \bi{v} = \bi{v}_{\rm gc} + v_{\perp} ( -\sin \theta \hat{\bi{e}}_{1} - \cos \theta \hat{\bi{e}}_{2}),\end{equation} \tag{ 1 }$$

where v_gc is the guiding center velocity. The second term represents the ion gyromotion with speed v_⊥, gyrophase θ and unit vectors perpendicular to the local magnetic field, $\hat {\bi {e}}_{1} \times \hat {\bi {e}}_{2} = \hat {\bi {b}}$, where $\hat {\bi {b}} = \bi {B}/B$ is the magnetic field direction and B = |B| is its magnitude.

The Hamiltonian equations of motion [2, 14] for the ions are integrated by XGC0 for a discrete set of marker particles; the ion guiding center distribution function is written as

$$\begin{equation} f(\bi{x}, \bi{v}) = \sum_{k = 1,N} w(k) \delta[ \bi{x} = \bi{x}(k)] \delta[ \mu = \mu(k) ] \delta[ \rho = \rho(k) ], \end{equation} \tag{ 2 }$$

where the sum is over the N marker particles, each of which is assigned a statistical weight w(k). The marker particles are localized in phase space first by the normalized velocity of the ions parallel to the magnetic field, ρ = mv_∥/qB, with m and q being the ion mass and charge, respectively. The second coordinate is the magnetic moment, defined as μ = mv²_⊥/(2B) with v_⊥ defined by (1). Each marker particle effectively represents a portion of phase space occupied by a large number of actual particles, as indicated by the statistical weight.

Conservation of the total weight of the XGC0 marker particles, $w_{\mathrm { tot}} \equiv \sum _{k} w(k)$, can be checked by computing, at time t_i, the left-hand side of

$$\begin{equation} w_{\rm tot}(t_{i}) + \sum_{j=1}^{i} L(t_{j}) - \sum_{j=0}^{i} I(t_{j}) = w_{\rm tot}(t_{0}) \end{equation} \tag{ 3 }$$

and comparing with the total weight at the beginning of the run, w_tot(t₀). The quantity L represents the ion losses to the material boundary (the 'limiter') at each time step, and I represents the sources (or sinks) due to neutral–plasma interactions, e.g. due to ionization of atoms. Expressions analogous to (3) can be written for the energy and for each of the components of the momentum vector.

A strict test of conservation would evaluate I using the actual kinetic collisions in the XGC0 collision routine, I_p. But, the DEGAS2 collision routine also yields its own value for I, which we denote by I_n. Since the neutral distribution produced by DEGAS2 reflects the effects of the source (or sink) I and since this neutral distribution is going to be used as the background for subsequent execution(s) of the XGC0 collision routine, any inequality between I_p and I_n represents a violation of mass conservation. The magnitude of |I_p − I_n|/[(|I_p| + |I_n|)/2] quantifies any deficiencies in our plasma–neutral interaction model. The amount of error introduced into (3) by a difference between I_p and I_n provides a quantitative measure of its impact on the overall plasma evolution.

Because the two collision routines run independently, conservation errors may not manifest themselves as obvious numerical problems. The physics results produced by such a simulation would be in error, but we would in general be unable to detect that error since we do not have in hand a comparably detailed, but independent, calculation. We can instead only insist that the algorithm be consistent with the expected conservation laws by explicit evaluation of (3) and monitoring the relative values of I_p and I_n.

2.1. Ion distribution function

The present implementation of the coupling algorithm utilizes a drifting Maxwellian distribution to characterize the background species in a collision. Given the kinetic characteristics of both plasma and neutral species, highlighted above, this represents an approximate treatment. Its principal benefits are simplicity, ease of implementation and relatively modest computational requirements. In section 4, we quantify the consequences of this approximation by examining energy conservation in the coupled codes, both for neutral–plasma exchanges and for the plasma as a whole.

The expressions used to evaluate the parameters of the ion distribution have been constructed so as to ensure that the treatment of charge exchange collisions is equivalent in the two routines. The DEGAS2 collision routine samples candidate charge exchange ions from the ion distribution,

$$\begin{equation} f_{\rm i}(\bi{v}) = \frac{n_{\rm i}}{\left( 2 \pi T_{\rm i} / m \right)^{3/2}} \exp \left[ \frac{-m}{2 T_{\rm i}} ( \bi{v} - \bi{v}_{\rm f,i})^{2} \right], \end{equation} \tag{ 4 }$$

where n_i, T_i and v_f,i are the ion density, temperature and flow velocity.

The ion density is computed in an obvious manner,

$$\begin{equation} n_{\mathrm{i}}(i_{\rm tri}) = \frac{1}{V(i_{\rm tri})} \sum_{k \in i_{\rm tri}} w(k), \end{equation} \tag{ 5 }$$

where i_tri is the index of a (triangular) mesh cell and V is its volume. The sum is over all marker particles k in cell i_tri.

We then require that the ion flow velocity used in (4) equal the average of the ion velocity (1) over the marker distribution (2):

$$\begin{equation} \bi{v}_{\mathrm{f,i}}(i_{\rm tri}) = \frac{1}{n_{\mathrm{i}}(i_{\rm tri}) V(i_{\rm tri})} \sum_{k \in i_{\rm tri}} w(k) \bi{v}_{\mathrm{gc}}(k). \end{equation} \tag{ 6 }$$

The ion temperature can be expressed in terms of the average ion energy, 〈E〉,

$$\begin{equation} T_{\rm i} = \frac{2}{3} \left ( \langle E \rangle - \frac{1}{2} m \bi{v}_{\rm f,i}^{2} \right ). \end{equation} \tag{ 7 }$$

Averaging the ion energy over the marker distribution results in an expression analogous to (6) so that

$$\begin{equation} T_{\mathrm{i}}(i_{\rm tri}) = \frac{m}{3} \left [ \frac{1}{n_{\mathrm{i}}(i_{\rm tri}) V(i_{\rm tri})} \sum_{k \in i_{\rm tri}} w(k) v^{2}(k) - v_{\mathrm{f,i}}^{2}(i_{\rm tri}) \right ]. \end{equation} \tag{ 8 }$$

To minimize statistical fluctuations associated with the gyrophase in v (1), we also use the definition of μ and write

$$\begin{equation} v^{2}(k) = \frac{2 \mu(k)B}{m} + v_{\mathrm{gc}}^2(k) - 2 v_{\perp}(k) \left[ \bi{v}_{\mathrm{gc}}(k) \cdot \hat{\bi{e}}_{1} \sin \theta (k) + \bi{v}_{\mathrm{gc}}(k) \cdot \hat{\bi{e}}_{2} \cos \theta (k) \right] \end{equation} \tag{ 9 }$$

in the expression for T_i. The gyrophase terms in this expression will average out upon integration over the distribution function.

2.2. Vacuum vessel filling mesh

The computational mesh used by both collision routines in compiling and exchanging information about the neutral and plasma distributions begins with a specification of the tokamak vacuum vessel as a closed polygon. Also required is a characterization of the magnetic equilibrium for the plasma discharge and time when that discharge is being simulated.

The most common representation of such equilibria is as a set of poloidal magnetic flux values ψ on a regular, rectangular grid spanning the poloidal plane of the tokamak. Such equilibria are usually obtained as output files from the EFIT code [15]. As XGC0 runs, the plasma pressure and electric field evolve to satisfy force balance for that magnetic equilibrium. In other applications of XGC0, the Grad–Shafranov equation is solved during the simulation to allow the magnetic equilibrium to evolve in response to transient events [16]. The procedure described in the remainder of this section would need to be automated before DEGAS2-based neutral transport could be included in those calculations.

To ensure adequate spatial resolution in the mesh cells adjoining the boundary and to provide a convenient basis for compiling particle and heat fluxes to all segments of the boundary, we first divide it so that no segment is longer than a specified maximum length.

Most of the plasma heat and particle fluxes will be deposited, however, in a narrow (of the order of centimeters wide) band adjacent to the points at which the magnetic separatrix intersects the boundary; these locations are the 'strike points' [11]. To adequately resolve these fluxes and to allow for the fact that the strike point locations vary with the equilibrium, we allow the user to specify a grid of poloidal magnetic flux values, normalized to the magnetic flux of the separatrix surface ψ_n ≡ ψ/ψ_sep, that can be mapped onto the boundary, further refining it.

This normalized poloidal magnetic flux grid is characterized by the width of the two grid segments closest to and farthest from the separatrix, as well as the total width of the grid; this approach is analogous to that used by Marchand and Dumberry [17]. These three parameters allow one to derive the coefficients of a quadratic expression for the individual poloidal magnetic flux values, which can then be mapped onto the boundary by interpolation from the equilibrium. Independent sets of three parameters are utilized for the SOL (ψ_n > 1) and PFR (ψ_n < 1) sides of the separatrix. Any boundary points added in the initial subdivision step that fall within these regions are removed to simplify interpretation of the resulting heat flux profiles. The discretization of the boundary in the divertor region for example runs described in section 3 is shown in figure 1.

Zoom In Zoom Out Reset image size

This boundary is also used in loading and pushing the XGC0 marker particles to provide consistency throughout the code. Determination of whether a particle is inside or outside the boundary is made using the Jordan Curve Theorem [18]. If a just-advanced particle is found to be outside the boundary, the particular segment crossed is identified by interpolating between the particle's previous and current locations. Its weight and energy are added to the boundary and heat flux arrays, and the particle is removed from the problem.

The extreme disparity between plasma transport parallel and perpendicular to the magnetic field results in parameters (e.g. (8)) having gradients much larger in the direction perpendicular to magnetic flux surfaces than along them. Adequately resolving these gradients with an arbitrary spatial discretization would require a prohibitively large number of cells. For this reason, fluid edge plasma transport codes [9, 19] are based on quasi-orthogonal meshes in which one coordinate is aligned with magnetic flux surfaces. Such meshes suit our purposes as well, once the modifications described below have been incorporated. We ensure that the spacing of the aligned surfaces and the second, quasi-orthogonal, are fine enough to resolve spatial variations in the plasma and neutral parameters. Here, we employ the Carre code [17], originally developed to produce meshes for B2 and used previously in DEGAS2 [20] applications.

Carre's algorithm for establishing and optimizing the quasi-orthogonal surfaces in the vicinity of the divertor targets, i.e. the 'ends' of the magnetic flux surface contours, results in a boundary that differs slightly from the one that was provided on input. To ensure that the Carre mesh resides entirely within the vacuum vessel boundary used elsewhere within XGC0 (e.g. to identify lost particles), the boundary input to Carre is shifted inward by 5 mm. This shift is implemented as a 'polygon offset' using the Clipper library [21, 22].

The other three boundaries of the Carre mesh coincide with magnetic flux surfaces. The innermost boundary is a magnetic flux contour located somewhere between the separatrix and the magnetic axis; the particular contour used depends on the physics being studied. The outermost boundary is the 'last' (i.e. encompassing the greatest poloidal magnetic flux) magnetic flux surface that provides a continuous connection between the inner and outer divertor targets. The volume between this flux surface and the separatrix provides the most direct connection to the core plasma and is referred to as the 'scrape-off layer' (SOL) [11]. The volume between the separatrix legs near the strike points is magnetically isolated from the rest of the plasma; this is the 'private flux region' (PFR). The outermost flux surface in the PFR is determined via criteria analogous to those for the SOL boundary.

The resulting mesh fills most of the vacuum vessel volume. We tile the space between it and the vessel boundary, as well as the region interior to the mesh's inner boundary, with triangles via the Triangle utility [23]. Note that these triangles also fill the narrow gap remaining between the 'offset' (i.e. the bottom of the Carre mesh) and actual boundaries. Each of the Carre mesh quadrilaterals is split into two triangles so that the mesh is comprised entirely of triangles. Finally, the triangles are renumbered by a Hilbert space filling curve to make the particle localization procedure more efficient [14]. The mesh used for the calculations described in the remainder of this paper is depicted in figure 2.

Zoom In Zoom Out Reset image size

This method for discretizing the vacuum vessel volume allows all 'boundary conditions' to be established at material surfaces in a manner consistent with our understanding of the physics of plasma– and neutral–material interactions. Relative to the original XGC0 neutral transport model, an increased number of particles is required to obtain adequate statistics because plasma quantities are resolved in two dimensions instead of being averaged over magnetic flux surfaces. The capability for averaging the plasma moments over closed flux surfaces could be added to the XGC0–DEGAS2 code should it be deemed necessary.

The calculation of integrals over the ion marker distribution in other XGC0 subroutines (e.g. for the ion–ion collision operator) is performed via linear interpolation onto flux surfaces or spatial locations in the manner typical of particle-in-cell codes. By instead evaluating (8) as the sum over all simulation particles j in a mesh triangle i_tri, the result becomes equivalent to the volume integrals used in compiling the output of the DEGAS2 routine. Consequently, the ionization rate computed by the XGC0 collision routine would be exactly equal to that produced by the DEGAS2 routine if one were to hold the ions fixed in between the two subroutine calls. The number of actual ionizations processed during a given time interval, i.e. the effective ionization rate, will converge to this rate in the limit of a large number of marker particles or a long time interval. This implies that the two routines should conserve mass, at least in a statistical sense. In practice, the marker particles are advanced one time step between the calls to the DEGAS2 and XGC0 collision routines. For the simulations described in section 4, the time step and number of particles are such that the difference between the ionization rates cannot be discerned from the statistical error. While runs with many more particles could, in principle, render this difference statistically significant, the neutral density profile used in the XGC0 collision routine is automatically rescaled (section 2.5) to account for such deviations so as to conserve mass globally.

2.3. Plasma collision algorithm

The null collision method [24] is used to process plasma–neutral collisions in the XGC0 routine. The principal idea behind this technique is to sample collisions along a particle's trajectory using a constant reaction rate ν_max, but do nothing to the particle for 1 − ν(v)/ν_max of these (null) collisions. The algorithm is exact in the statistical sense and is made particularly simple if ν_max can be taken as the maximum reaction rate over the problem volume. The total reaction rate for an ion marker particle with velocity v is the sum over all pertinent neutral collision processes,

$$\begin{equation} \nu(\bi{v}) = \sum_{j} \nu_{j}(\bi{v}) \equiv n_{\mathrm{n}} \sum_{j} \langle \sigma v_{\rm rel} \rangle_{j}(\bi{v}), \end{equation} \tag{ 10 }$$

where n_n is the local neutral density and v_rel = |v − v_n| is the relative velocity of the ion and a colliding neutral. The angle brackets indicate an average of the cross section for process j, σ_j, over the neutral distribution function.

Neutral–plasma collisions in XGC0 are sufficiently infrequent that the neutral collision routine is usually invoked only every N_ion time steps with N_ion > 1. This approach represents an approximation to the original null collision algorithm in that the collisions are processed with the particles at the final x(t = t₀ + N_ionΔt) and v(t = t₀ + N_ionΔt) of an uncollided trajectory rather than at a sampled t₀ ⩽ t_c ⩽ t₀ + N_ionΔt.

The simulations described in sections 3 and 4 use N_ion = 20 to obtain reasonable statistics for the conservation analysis. We do not consider larger values so as to ensure that (N_ionΔt)ν_max < 1. When this inequality is violated, marker particles could potentially undergo multiple real collisions during N_ionΔt, implying significant modification of the local ion distribution function in those regions. Note that our routine does allow for the marker particles to experience multiple collisions during the interval N_ionΔt to handle transients with large ν_max, particles with very short sampled collision times, etc.

To quantify the error incurred by having N_ion > 1, we consider the change in the volume integrated core particle and energy source rates, which are of interest in tokamak fueling studies [25]. In reducing N_ion from 20 to 1, the core particle source rate falls by 1.3%, and the core energy source rate by 0.2% (both computed by the DEGAS2 routine). The former figure is small in comparison with other uncertainties and approximations made in these calculations. The latter is smaller than the statistical error in the total energy source rate and is thus insignificant.

2.4. Treatment of ionization and charge exchange

Only deuterium ions and atoms are simulated in this initial version of the coupled plasma–neutral code. As noted previously, the relevant interactions are ionization and charge exchange. The ions are used as a quasi-neutral proxy for the electrons in simulating ionization. That is, ionization is processed using the ion locations and a fixed electron temperature profile specified on input. The ionization rate, a function of the electron temperature and density, is computed using DEGAS2's atomic physics routines, which in turn rely on data obtained from a collisional radiative model [26, 27]. The ionization rate associated with a given ion thus depends on its location in the problem, but not its velocity.

The charge exchange cross section, in contrast, depends only on the relative velocity of an ion and colliding neutral. The probability for a given ion undergoing a charge exchange collision is based on the corresponding Maxwellian averaged reaction rate ν_cx, as in (10). This rate is a function of the local neutral temperature and the ion marker particle's velocity relative to the neutral flow velocity. We randomly sample a gyrophase for the ion marker particle and use (1) to obtain its velocity when evaluating this reaction rate.

The velocity dependence of the cross section is enforced by sampling the colliding neutral via a rejection technique, as depicted in figure 3. A neutral velocity v_n is sampled from the local (Maxwellian) distribution, and the charge exchange cross section σ_cx(v_rel) is computed. A uniform deviate ξ_cx is obtained and the sampled neutral accepted if $\xi _{\mathrm { cx}} < \sigma _{\mathrm { cx}} v_{\mathrm { rel}} / \max (\sigma _{\mathrm { cx}}v_{\mathrm { rel}})$, where $\max (\sigma _{\mathrm { cx}}v_{\mathrm { rel}})$ is the maximum value of σ_cxv_rel for all values of v_rel.

Zoom In Zoom Out Reset image size

The charge exchange cross section and rate data are accessed via DEGAS2's atomic physics routines, which in turn interpolate from the original data computed by Krstic and Schultz [28] for quantum mechanically 'indistinguishable' particles. A common misconception is that charge exchange and elastic scattering must be simulated separately for a system such as D⁺ + D. However, as is explained in [28], elastic scattering and charge exchange represent the small and large angle scattering limits, respectively, of the same process. For interaction energies above roughly 1 eV, the differential scattering cross section peaks strongly near 0 and 180° in the center-of-mass frame. The latter peak is equivalent to classical charge exchange in which the atom and ion effectively swap velocity vectors. The peak near $0^\circ$ represents elastic scattering. Because those collisions yield small collision angles, the momentum transferred between the ion and atom is negligible and the process can be ignored in a transport simulation. Since the emphasis of the present XGC0–DEGAS2 calculations is on physics close to the separatrix, i.e. well above 1 eV, we take this approach here and employ the 'spin exchange' cross section determined by Krstic and Schultz. For energies <1 eV, the differential scattering cross section is significant for a range of scattering angles, blurring the distinction between 'elastic scattering' and 'charge exchange'. A more complicated sampling procedure would be required in this case [29].

We set the value of ν_max in the XGC0 collision routine, figure 3, to be the maximum value of the sum of the ionization and charge exchange rates over the entire volume for simplicity. In the problems run to date, the time spent in the plasma–neutral collision routines is too small compared with the rest of the code (see section 2.5) to justify the implementation of a more efficient algorithm, e.g. as in [30].

2.5. Consistent source and time dependence

These example simulations are based on DIII-D discharge 96333 [33]. The EFIT equilibrium for this shot at 3300 ms has served as a standard reference discharge for XGC simulations [34, 35]. As in those papers, we assume initial H-mode-like profiles with a pedestal density of 5 × 10¹⁹ m⁻³ and a temperature of 1 keV. However, for the present simulations we employ a somewhat lower electron temperature, figure 4(b), consistent with DIII-D H-mode profiles at this density (e.g. [36]). The only other adjustable parameters in these simulations are a collisionless gyroviscosity coefficient [37] of 5 × 10⁻² m² s⁻¹ used by XGC0 in pushing its ions and a recycling coefficient, noted in section 2.5, of 90%. The latter is known to be close to or equal to unity in most fusion devices; the value of 90% was otherwise chosen somewhat arbitrarily.

Zoom In Zoom Out Reset image size

The ion marker particles are tracked for 20 ion transit times (10⁴ time steps, 1.56 ms) until all transients have died off. The density pedestal builds up and the gradients steepen, as in [2]. Note that the ion temperature drops from its initial value of 1 keV at the top of the pedestal since we have, for simplicity, not included a heat source from the core to offset the ion heat loss to the boundary and neutral cooling. More detailed XGC0–DEGAS2 simulations of particular experimental discharges would include an appropriate heat source, kinetic electrons, impurity species, as well as turbulent diffusion with anomalous diffusion and recycling coefficients calibrated to match experimental profiles.

The ion fluxes to the divertor floor are plotted in figure 5 as a function of distance along the boundary. These ion losses are solely due to ion motion along and across the open field lines, as well as to ion orbits intersecting the material boundary, with the latter occurring primarily in the low poloidal field region around the X-point [2].

Zoom In Zoom Out Reset image size

XGC0 scales very efficiently on massively parallel computers, up to peta-flop levels. However, these simulations are performed on a smaller, local Linux cluster utilizing 32 cores on 8 nodes, a practical amount of computational power given that we are not using XGC0's kinetic electron capability. About 19 h are required to process the 10 000 time steps in these runs. All of the neutral-related computations, including both the XGC0 and DEGAS2 collision routines, occupied about 20% of the time. In comparison, the routines that advance the ion marker positions and velocities typically consume more than half of the time. A total of 7.6 million ions and 0.96 million neutral particles are tracked through the 17 325 mesh cells in the geometry used for the plasma–neutral calculation. Less than 1 GB of memory is needed on each core.

To assess the conservation properties of the plasma–neutral coupling algorithm, we compare the global integrals of the neutral–plasma exchange rates computed by the two collision routines (I_p and I_n in section 2). The time variation of the ion particle and energy source rates are shown in figure 6(a). The period prior to 0.4 ms is dominated by transient relaxation from the input profiles. Note, in particular, the spike around 0.16 ms (1000 time steps) corresponding to the initial update of the particle flux poloidal profile (section 2.5). Charge exchange results in no net production of ions so that the ion source rate represents only the effects of ionization. As expected from the design of the algorithm (section 2), the source rates from the two collision routines match to within the statistical variations.

Zoom In Zoom Out Reset image size

Note that the XGC0 rates are computed via the equivalent of a 'collision estimator' [27], i.e. based entirely on the actual collisions of the XGC0 marker particles. Since only a small fraction of the XGC0 ions experience such collisions in the time interval N_ionΔt, the resulting rates are relatively noisy. Note, however, that these rates are computed only for diagnostic purposes, and the statistical variations have no impact on the outcome of the simulation. The DEGAS2 rates, in contrast, are compiled with a 'track length estimator' [27] in which quantities are integrated along the neutral particle trajectories. This approach yields more precise results than the collision estimator in the low collision frequency limit.

The negative values of the ion energy source rates are expected since, on average, the warmer plasma ions are transferring energy to the cooler neutral atoms via the charge exchange process. Ionization serves as an energy source since it adds new ions, albeit cold ones, to the population. The cooling associated with electron impact excitation and ionization of atoms is not accounted for in these rates since the electron temperature profile is held fixed.

The principal result of interest in figure 6(a) is that the energy source rates from the two collision routines differ significantly, by ∼30% at the end of the run. To assess the global impact of this error, we compute the energy analogue of (3), shown in figure 6(b), with one curve representing the XGC0 collision routine (I = I_p) and one representing the DEGAS2 collision routine (I = I_n). The steady divergence of the two curves in time is expected since one is integrating the roughly constant difference between the energy source rates in 6(a). Note, however, the suppressed zero in this plot so that even at the end of the simulation, the time integrated discrepancy represents <2% of the initial plasma energy.

We first show that this discrepancy in the energy source rates is due entirely to charge exchange and not ionization. In figure 7(a) we plot the charge exchange rate integrated over the entire plasma (left vertical axis) and the energy 'source' only due to charge exchange (right axis). The difference between the XGC0 and DEGAS2 energy source curves in figure 7(a) is on average within 3% (relative to the total ion energy source) of the difference between the XGC0 and DEGAS2 energy source curves in figure 6(b). This is less than the estimated statistical error in the energy source rates of 5%.

Zoom In Zoom Out Reset image size

We next argue that the discrepancy is associated with the velocity dependence of the charge exchange cross section and the non-Maxwellian character of the distribution functions of the native kinetic species in the XGC0 and DEGAS2 codes. To demonstrate the latter point, we examine the expected rates for the charge exchange process and the associated energy exchange in the XGC0 collision routine for the given neutral background. The charge exchange reaction rate evaluated at the local neutral temperature and flow velocity is

$$\begin{equation} \langle \sigma_{\rm cx}v \rangle_{\rm M,n}(\bi{v}; T_{\rm n},\bi{v}_{\mathrm{f,n}}) = \int \, \mathrm{d}\bi{V} f_{\rm M,n}(\bi{V}; T_{\mathrm{n}}, \bi{v}_{\mathrm{f,n}}) \sigma_{\rm cx}(|\bi{v} - \bi{V}|) |\bi{v} - \bi{V}| \end{equation} \tag{ 11 }$$

$$\begin{equation} = I_{0,0}(|\bi{v} -\bi{v}_{\mathrm{f,n}}|,T_{\mathrm{n}}), \end{equation} \tag{ 12 }$$

where

$$\begin{eqnarray} I_{\ell,\mathrm{n}}(u,T)& \equiv& \frac{\exp(-u^{2} / u_{t}^{2})}{\sqrt{\pi} u u_{t}^{n+1}} \int_{0}^{\infty} \, \mathrm{d}v_{\rm rel} v_{\rm rel}^{2+n} \sigma_{\ell}(v_{\rm rel}) \left \{ \exp \left [ \frac{ -v_{\rm rel}(v_{\rm rel} - 2 u)}{u_{t}^{2}} \right ] \right. \nonumber \\ & & - (-1)^{n} \left. \exp \left [ \frac{ -v_{\rm rel}(v_{\rm rel} + 2 u)}{u_{t}^{2}} \right ] \right \} \end{eqnarray} \tag{ 13 }$$

is an integral that arises in taking the moments of atomic physics cross sections over a Maxwellian distribution function [38]. Only the total cross section, corresponding to ℓ = 0, is of interest for present purposes.

We can then compute in each XGC0 mesh triangle the velocity space integral of this rate over the local ion marker distribution, (2),

$$\begin{equation} \langle \langle \sigma_{\rm cx} v \rangle_{\rm M,n} \rangle_{\mathrm{i}} = \int \, \mathrm{d}\bi{v} f_{\mathrm{i}}(\bi{v}) \langle \sigma_{\rm cx} v \rangle_{\rm M,n}(\bi{v},T_{\rm n},\bi{v}_{\mathrm{f,n}}). \end{equation} \tag{ 14 }$$

The sum of this over all mesh triangles yields the charge exchange rate curve labeled 'Expected' in figure 7(a).

We can similarly compute the expected energy exchange rate

$$\begin{eqnarray} \langle \langle \Delta E_{\rm cx} \rangle_{\rm M,n} \rangle_{\mathrm{i}} &=& \frac{m}{2} \int \, \mathrm{d}\bi{v} f_{\mathrm{i}}(\bi{v}) \int \, \mathrm{d}\bi{V} f_{\rm M,n}(\bi{V}; T_{\mathrm{n}}, \bi{v}_{\mathrm{f,n}}) \nonumber \\ &&\times \sigma_{\rm cx}(|\bi{v} - \bi{V}|)|\bi{v} - \bi{V}|(V^{2} - v^{2}). \end{eqnarray} \tag{ 15 }$$

The average deviation of the 'Expected' charge exchange rate (energy exchange rate) from the 'XGC0' rate over the plotted interval is 1% (4%), comparable to the estimated statistical error of 1% (5%). This confirms that the null collision algorithm is working properly and that the ions are undergoing collisions at the expected rate.

If we instead replace f_i with a Maxwellian distribution at the local T_i and v_f,i, we obtain

$$\begin{equation} \langle \langle \sigma_{\rm cx} v \rangle_{\rm M,n} \rangle_{\rm M,i} = I_{0,0}[w_{\mathrm{f}},(T_{\rm i}+T_{\rm n})], \end{equation} \tag{ 16 }$$

where w_f ≡ |v_f,i − v_f,n|, and

$$\begin{eqnarray} \langle \langle \Delta E_{\rm cx} \rangle_{\rm M,n} \rangle_{\rm M,i} &= & \frac{2(\bi{v}_{\mathrm{f,n}} -\bi{v}_{\mathrm{f,i}})}{w_{\mathrm{f}}^{2}} \cdot (\bi{v}_{\mathrm{f,i}} T_{\mathrm{n}} + \bi{v}_{\mathrm{f,n}} T_{\mathrm{i}}) \nonumber \\ & & \times \left(\frac{w_{\mathrm{f}}}{\sqrt{2(T_{\mathrm{n}} + T_{\mathrm{i}})/m}} I_{0,1} - \frac{I_{0,0}}{2}\right) + (T_{\mathrm{n}} - T_{\mathrm{p}}) I_{0,2}, \end{eqnarray}\noindent \tag{ 17 }$$

where all of the I_ℓ,n are evaluated as I_ℓ,n[w_f,(T_i + T_n)]. The corresponding curves in figure 7(a) are labeled 'Maxwellian'.

The 'Maxwellian' charge exchange and energy source rates differ from both of the corresponding XGC0 and DEGAS2 curves, implying that the non-Maxwellian character of both atoms and ions is significant.

The other half of the argument begun above is that such kinetic effects enter into these rates via the velocity dependence of the charge exchange cross section. To demonstrate this explicitly, we repeat the simulation with a cross section σ_cx∝1/v_rel so that the kinetic reaction rate σ_cxv_rel is independent of the relative particle velocity. In particular, we set σ_cxv_rel = 3.8 × 10⁻¹⁴ m³ s⁻¹. The resulting energy source rates for the XGC0 and DEGAS2 routines, shown in figure 7(b), now match, with an average deviation of 5% that is comparable with the estimated statistical error of 5%.

We can use this simulation to further demonstrate, as we did with figure 6(b), that the global characteristics of the XGC0 simulation are insensitive to the details of the neutral–plasma interactions and to the energy source rate discrepancy in particular. First, the density profile from this run, included in figure 4(a), is virtually indistinguishable from that of the baseline simulation. The density at the top of the pedestal is only 2% higher than in the baseline simulation. Moreover, the volume integrated particle and energy source rates in the two runs (figures 6(a) and 7(b)) differ by 8 and 3%, respectively.

Because the discrepancy in energy source rates is associated with the non-Maxwellian character of the ion and neutral distribution functions, it cannot be mitigated by increasing the number of simulation particles, decreasing the time step, or any other straightforward means. Rather, one would need to use a more general representation for the background distribution functions than the assumed drifting Maxwellian (see section 2).

We anticipate that the momentum exchange rates in these simulations exhibit discrepancies similar to that described above for the energy exchange. However, they cannot be discerned from the analogous plots due to greater noise levels.

We have described a new coupled kinetic plasma–kinetic neutral transport code that can be used to address critical problems in edge tokamak plasmas. The principal approximation in the algorithm is the use of drifting Maxwellian distribution functions to exchange information between the two codes. The non-Maxwellian character of both the neutral and plasma distribution functions combines with the velocity dependence of the charge exchange cross section to cause energy non-conservation in the neutral–plasma exchanges by as much as 30%. However, these differences represent <2% of the total energy in the system. Moreover, a simulation with a velocity-independent charge exchange rate, which does conserve energy, yields almost indistinguishable core plasma profiles. We thus consider this energy exchange discrepancy a kinetic subtlety that will not have a significant impact on physics results, e.g. on studies of H-mode pedestal dynamics. Consequently, we are beginning to apply the code and continuing to improve its physics capabilities.

Planned upgrades include compiling the energy distribution of ions lost to the material boundary and using that distribution in determining the kinetic character of the neutral source, as described in section 2.5. This will allow molecules to be introduced into the problem in a realistic manner. We also plan to incorporate XGC0's kinetic electrons [12] into the XGC0 collision routine by using their temperature to determine the ionization rate; the electron cooling associated with the excitation and ionization of atoms can then be accounted for in the simulations, a significant effect. The addition of impurities to the coupled XGC0–DEGAS2 code is also being pursued.

Over a longer time frame, or when physics studies demand greater fidelity, we will investigate ways of improving energy conservation in the XGC0–DEGAS2 coupling. The most obvious approach would be to employ more general representations for the distribution functions passed between the codes, e.g. spline fits in velocity space.

This work was supported by US DOE contract no. DE-AC02-09CH11466.