{'fraction_non_alphanumeric': 0.022247745448358006, 'fraction_numerical': 0.005530032329419772, 'mean_word_length': 4.95013920526449, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing nonlinear magnetohydrodynamic edge localized instabilities in fusion plasmas

. The onset and nonlinear evolution of Edge Localized Modes (ELMs) in toroidally conﬁned plasmas are known to shed thermal energy from the edge of the conﬁnement region, and may also aﬀect the core plasma through nonlinear mode coupling. The physics of this process is not well understood, although the concomitant large bursts of thermal energy transport are a major concern for future burning plasma experiments. The evolution of ELMs is inherently nonlinear and analytic approaches are limited by the complexity of the problem. Save a handful of recent important theoretical works, the nonlinear consequences of ELMs are mainly unexplored. Recent developments in the NIMROD code [http://nimrodteam.org] have enabled the computational study of ELMs in tokamaks in the extended magnetohydrodynamic (MHD) framework, and a new initiative was formed to understand the basic physics of their nonlinear evolution. The results of these investigations are presented for both model equilibria and accurate reconstructions from the DIII-D experiment at General Atomics [http://fusion.gat.com/diii-d/]. These results show a ﬁlamentary high temperature structure propagating radially outward, which is strongly damped by experimentally relevant toroidal ﬂow shear. Two ﬂuid and gyroviscous terms are included linearly as a preliminary indication of these important physical eﬀects, and stabilization of higher wave number modes is observed.


Introduction
With recent advances in computer hardware and numerical algorithm efficiency, it has become viable for large-scale computational modeling of plasma dynamics to play an important role in the design and analysis of magnetic confinement fusion devices. An example of an important fusion problem that is being addressed by this approach is the onset and nonlinear evolution of Edge Localized Modes (ELMs) and their effect on global confinement. ELMs are high frequency bursting perturbations observed in experimental discharges relevant to fusion. These modes shed thermal energy from the edge of the confinement region and can release enough energy to damage or destroy plasma-facing components of future burning plasma experiments. They may also drive instabilities in the core plasma through nonlinear mode coupling. Recently, the NIMROD Team (http://nimrodteam.org) has begun a numerical study of the global dynamics of ELMs in tokamaks.
Toroidal magnetic confinement systems are characterized by a helical magnetic field where the field winding ratio of toroidal to poloidal turns (the safety factor q) is a function of magnetic flux, typically increasing from q ≈ 1 in the core to q ≈ 4 at the edge, with an inreasing gradient toward the edge. Resonance in perturbed fields occurs at rational surfaces in q, where m/n = q with m the poloidal mode number and n the toroidal mode number, m and n being integers. Near the edge of the plasma the density of rational surfaces increases as the gradient in q increases. Thus ELMs typically have a large number of interacting, overlapping linear modes contributing to their nonlinear evolution.
Previous studies have found the onset of ELMs to be well described as a linear MHD instability termed peeling-ballooning [1,2]. The peeling component of the instability is the free-energy drive due to the current density gradient, and the ballooning component is the free-energy drive due to the pressure gradient. At the edge of the confinement region of tokamaks, poloidal flows reduce energy transport and allow a steep pressure gradient to develop. The strong pressure gradient drives highly localized charge-current density. The pressure and current gradients are also coupled through the MHD equilibrium force-balance (described below), but they effect distinct characteristics in the linearly unstable modes that are excited. An individual ELM will typically display a combination of both sets of linear characteristics.
However, the evolution of ELMs is inherently nonlinear, the consequences of which are relatively unexplored [3]. Analytic methods for the study of the nonlinear evolution of ELMs are limited by the complexity of the problem. This study focuses on computational stability analyses using the NIMROD code [4]. NIMROD solves linear and nonlinear extended-MHD equations as initial-value computations with a mesh of finite elements for the poloidal (R-Z) plane and finite Fourier series for the toroidal direction. The time-advance is described in Ref. [5]. Among numerous options one can choose to solve either the single fluid MHD equations or include Hall and gyroviscous terms to solve the two fluid MHD equations. The single fluid form of the equations is ∂n ∂t where n is the number density, V the fluid velocity, ρ the mass density, J the current density, B the magnetic field, p the pressure, Π the stress tensor (including a numerical −ρν∇ 2 V ), E the electric field, η the electric resistivity, T the temperature, Γ the ratio of specific heats, q the heat flux, and κ the thermal diffusivity. Within the code there are various options for the functional form of the diffusive coefficients. Here, we make use of either constant or temperature-dependent resistivity (η ∼ T − 3 2 from the Spitzer result [6]), while using time-independent thermal diffusivity coefficients and viscosity. A more accurate model can be necessary to address the subtle physics effects of differences in the electron and ion behaviors, such as electric field responses and the gyro-motion around magnetic field lines. Some of these effects can be critically important to stability, and can generate flows in the evolution of the plasma. A well known, more accurate ordering scheme than the MHD ordering used above yields the two fluid terms in the momentum equation and Ohm's law: where e and i refer to the electron and ion species, and and ⊥ refer to parallel and perpendicular to the magnetic field. For the Braginskii model [7] the gyroviscous stress tensor is whereb = B/B. More details can be found in Ref. [5].
The semi-implicit advance used in the NIMROD code allows numerically stable computations at time-steps that vastly exceed explicit numerical stability limits. For example, the nonlinear simulation of magnetic island evolution presented in Ref. [4] was run with a CFL number of 10 4 . While the nonlinear dynamics of the ELM computations occur over relatively short time-scales, the more accurate two-fluid equations have dispersive waves that further extend the range of time-scales. The time-advance requires the solution of sparse ill-conditioned linear systems at each time-step, and the parallel SuperLU library [8] provides an effective means for solving and preconditioning these systems.
In Section 2 we describe the two types of equilibria used in this study. Next in Section 3 we present results from linear computations using the single fluid equations and results including the Hall and Gyroviscous terms. Nonlinear results using the single fluid equations are presented in Section 4. Finally, in Section 5 the results are summarized and discussed.

Equilibria
The purpose of these computations is to gain intuition into experimental observations. The basic model for interpreting experimental data is ideal MHD, which are equations 1 with all diffusivities set to zero. To lowest order the plasma state is well represented by the steady state force balance solution, or equilibrium, of these ideal MHD equations. Accurate equilibrium reconstructions from the DIII-D experiment are used as initial conditions in these computations. To explore the context of experimental configurations, model equilibria are constructed with parameterized pressure and current, the two sources of free energy to the ELMs.
A model equilibrium configuration was constructed to be unstable only to a dominantly ballooning edge instability with low shear and shaping, and yet to reside in the zero dimensional parameter space of typical DIII-D discharges. The profiles for electron temperature and safety factor q are specified parametrically as initial settings for the solution. The equilibrium magnetic flux solution is computed keeping the total current and pressure fixed. This produces a reduction of shear in the near edge q profile. The localized edge current density and pressure gradient are clearly evident in the profiles shown in Fig. 1. This ballooning-dominant configuration is relatively easy to resolve radially because of the low edge shear.
Accurate equilibrium reconstructions using data from a high β ELMing DIII-D discharge, designed to study edge localized modes, were also used in this study. In ELMing discharges, the pressure and current profiles near the edge are constantly changing as the gradients increase and decrease due to competition between heating and instability induced transport. Typically, reconstructions are done at times in between the strong signal perturbations associated with ELMs. Thus, if the reconstructions and stability analysis are done to extreme accuracy, the results tend to be stable to ELMs. The typical approach to the study of edge modes is to take the equilibrium reconstructions as a reference point and parameterize the pressure and current profiles such that the stability boundaries can be mapped within the parameter space traversed by the experiment. An alternative approach that is used for our study is to take reconstructions at times just before the ELM becomes unstable, which tends to produce unstable equilibria directly.
The parallel current density, magnetic winding ratio (or safety factor, q), and pressure profiles for the equilibrium are shown in Fig 2. The poloidal magnetic field configuration has a single null below the confinement region. The non-dimensional geometric and profile properties are also similar to high performance tokamak operation (H-mode) in Alcator Cmod (http://www.psfc.mit.edu) and advanced scenarios for ITER (http://www.iter.org). This equilibrium represents an average of the data between several ELM cycles just before onset of the mode.

Linear Stability
With the model equilibrium shown in Fig. 1, the resistive MHD computations find that growth rates increase monotonically with n as seen in Fig. 3. This is consistent with the MHD theory of ballooning modes. The linear results from NIMROD compare favorably with results linear ideal MHD results for ballooning-dominant instabilities with eigenfunctions and growth rates that are in reasonable agreement. This implies that toroidal resolution in a nonlinear resistive MHD computation cannot be achieved when starting from an initial condition near equipartition in energy, as the highest mode numbers will be the fastest growing. Nonetheless, this equilibrium provides a suitable benchmark case for comparing linear results obtained with different numerical approaches.
In the specification of the resistivity profile in the case of the model equilibrium, for simplicity we use a constant resistivity in the core which quickly transitions to a large constant resistivity to model a cold, low density region surrounding the plasma, as shown in Fig. 1(a). The growth rates computed by NIMROD with non-ideal MHD depend weakly on the resistivity and viscosity in the vicinity of the mode, where there is a large gradient in resistivity, but the variation of growth rate with n agrees qualitatively with the linear ideal MHD results from the ELITE [9] and GATO [10] codes, as shown in Fig. 3. Varying the electrical resistivity shaping parameter by three orders of magnitude affects the growth rate values by less than a factor of 2, but it does not change the spectrum qualitatively. Here, the magnetic Prandtl number Pm, the ratio of kinematic viscosity to electrical diffusivity is held fixed while both the core and surrounding resistivity values are varied.
The effect of dissipation on the linear behavior of ideal modes can be either stabilizing or destabilizing, depending on whether damping or the breaking of the frozen field line constraint is dominant. The resistivity should increase near the edge and into the vacuum region, which can be important in the nonlinear evolution. But linearly the scaling is weak.
Although there is a weak scaling of the growth rates with the resistivity, the resistivity varies strongly across this region, and the vacuum resistivity setting can be important. In Fig. 3 the resistivity everywhere is decreased by a factor of 10, along with viscosity, from S = 1 × 10 6 to S = 1 × 10 7 , this decreases the growth rates. When the vacuum resistivity alone is increased the growth rates increase, and conversely we see an increase of growth rates with decreasing core resistivity while the vacuum resistivity is held fixed. This is due to the reduction of the viscosity while the Prandtl number is held fixed. Reduction of the viscosity to levels where the linear growth rates are weakly affected has not proven viable in the computation of the nonlinear evolution. In the experiment, there is some viscous damping to velocity fluctuations, which we approximate with diffusive viscosity.
For the experimental equilibrium reconstruction, the linear growth rate spectrum with respect to toroidal mode number for the single fluid MHD case is peaked at n = 12 with modes 5 − 20 being unstable, as shown in Fig. 4. The peak growth rate is 2.3 × 10 5 s −1 . This indicates that the respective equilibrium is well above the threshold for ideal linear instability. In the actual experiment the stability boundary is crossed and growth rates increase from zero with time. The ELM becomes nonlinear and modifies the equilibrium and the cycle repeats. With the equilibrium reconstructions taken just before the ELM with time averaging in the data, the reconstruction is linearly unstable.
Recent advances in NIMROD allow drift effects associated with the Hall electric field and gyroviscosity to be considered a described in Section 1. Conventional wisdom holds that two-  Figure 4. The growth rates as a function of toroidal mode number n for the DIII-D equilibrium reconstruction. Nonlinear growth rates are the exponential growth phase driven by coupling between exponentially growing linear modes, before the growth rates slow. fluid effects are stabilizing at large n-values [11], and this is often cited to explain experimental observations that very high mode numbers are not observed in ELM instabilities. This stabilization is due to the the drift velocity which is derived from the two fluid equations which is peaked in the edge region as shown in Fig. 5(a). NIMROD is based on a fluid model, which does not eliminate the high frequecy waves and accurately tracks the evolution. High frequency waves can be eliminated in the drift ordering, which can be convenient when used to investigate long time-scale dynamics.
We have calculated the linear growth rates of the equilibria shown in Fig. 2 with Hall and gyroviscous effects, and preliminary results are summarized in Fig. 5(b). Here we plot the linear growth rate as a function of the toroidal mode number for two different models (MHD and Hall-MHD with gyroviscosity) with anisotropic thermal diffusivity and resistivity relevant to experiment. The drift effect is clearly stabilizing to higher mode numbers while the lower mode numbers are weakly affected. The eigenfunctions are similar, but include a resonance condition associated with the poloidal flow, which is also seen with toroidal flow shear as described below in Fig. 8. A poloidal shift in the peak amplitude of higher mode number eigenfunctions has also been observed [5].
ELMs are characterized by a broad linear instability spectrum, which depends on the employed physics model. As these growing modes couple nonlinearly the spectral energy transfer will display rich dynamics that are important to understand for experimental interpretation. Differences in the model will cause differences in the linear spectrum, but how these differences affect the mode coupling and evolution will be even stronger because of long time-scales involved and additional nonlinear terms. Much of this physics remains unexplored and yet is critical for interpreting experimental results. We begin this exploration with nonlinear computations using these equilibria.

Nonlinear Resistive MHD
The nonlinear ELM computations solve the resistive MHD model with small amplitude perturbations in the initial conditions which excite the unstable ELMs. The computations follows the evolution well into the nonlinear phase and show the initial linear mode structure causing a large radial perturbation to the pedestal region with wavelength of the dominant linear mode. These form vortices as they move outward, curling in on themselves. Finally the structures break off and separate from the main plasma, forming filaments. These filamentary thermal structures extend from the separatrix and propagate radially to the wall as the computation progresses. Radial propagation of filamentary structures has been observed over the past few years on DIII-D [3] and MAST [12].
For the model equilibrium with ballooning dominant instability, the early nonlinear evolution of the edge mode shows an energy spectrum peaked at high and low n. The high n ∼ 20 modes have the largest growth rates and the low n ∼ 1 modes being nonlinearly driven by nearest neighbor coupling of high n modes. The mode structure shows fingers of high temperature flowing outward, as shown in Fig. 6(a). These computations lose toroidal resolution as the high n modes become large, and the computation terminates long before the filaments propagate to the wall.
For the computation using the realistic equilibrium, the electrical resistivity is temperaturedependent (Spitzer resistivity, based on the evolving toroidally averaged temperature). The magnetic Reynolds number S = 2 × 10 7 while the Prandtl number is P m = 10. The perpendicular thermal diffusivity is 25m 2 /s, and the parallel thermal diffusivity is 10 5 times larger. In general the dependence of the growth rates on thermal diffusivity is weak in the linear phase. In the nonlinear phase, thermal transport can be extremely important. We use larger values than the experiment for the convenience of computation.
The temporal evolution of the kinetic energy associated with each of the toroidal Fourier components for the realistic equilibrium is shown in Fig 7. There is a linear growth phase for 0.01ms < t < 0.03ms, followed by nonlinear saturation. Only modes in the range 5n20 exhibit linear growth, as shown in Fig. 4. Higher and lower mode numbers are nonlinearly driven by mode coupling and exhibit larger growth rates γ ∼ e (γ 1 +γ 2 )t for each mode driven by two other modes with growth rates γ 1 and γ 2 .
As described in the previous section, the growth rate spectrum is peaked around n = 12. This is attributed to the anisotropic thermal conduction, some viscous dissipation, and the large current-gradient (peeling component) at the edge of the confinement region in this equilibrium. Numerical convergence tests indicate that while greater resolution is needed to achieve quantitative accuracy, the peaked linear spectrum is qualitatively correct. Other modes may be linearly unstable but grow too slowly to show independent activity before nonlinear coupling becomes significant. Of particular importance is that the nonlinear coupling drives low n-fluctuations (including n = 1), in addition to high-n fluctuations, and low-n activity in the The evolution is more dynamic than what is evident from still figures. The nonlinear evolution of this mode drives a rapid loss of internal energy with approximately 70kJ( 10%) of the internal energy being lost within 50µs. Laboratory measurements indicate 15 − 20% energy loss during large ELM events in DIII-D, and we note that the numerical computation has not completed the ELM cycle. The internal energy is still decreasing at the end of the computation while pressure is lost over the entire pedestal region.
Analysis of the results finds the heat flux on the wall to be nonuniform with greatest intensity occurring in spots on the top and bottom of the chamber. The convective heat flux contribution to the heat loss is greater than that of the conductive heat loss, in contrast to NIMROD computational results on plasma disruption due to an internal mode [13]. The fact that the primary loss channel is convective rather than conductive is not inconsistent with laboratory measurements [14].
Toroidal flow shear can stabilize, destabilize, increase or decrease the linear growth rate spectrum of an ELM. A toroidal flow profile was included in the computations to explore the effect on the nonlinear evolution. The profile has constant flow across the core and a linearly decreasing flow with magnetic flux from ψ = 0.85 to 1.0 as shown in Fig. 8(a). The value of the flow at the very core was given as V c = 2 × 10 5 m/s, which is experimentally relevant. A series of computations with varying flow conditions is planned, but is prohibitively expensive to compute readily, especially nonlinearly, and we report here on a single case for comparison with no flow. When toroidal flow is included in the evolution, the growth rates increase by approximately 20% linearly with no qualitative changes in the spectrum. The linear eigenfunctions show a phase reversal of the perturbation radially as shown in Fig. 8. This is in qualitative agreement with ideal MHD results.
Nonlinearly, the effects of toroidal flow shear are extreme. There is a marked decrease in the amplitude of the perturbation poloidally as shown in Fig. 9. The filaments which propagated well out into the vacuum region for the no flow case are sheared off by the flow and do not propagate far into the vacuum region.

Discussion and Summary
Because of the deleterious effects of ELMs on plasma performance and operation, it is important to understand both their onset and the manner in which they deposit heat on the wall. A better understanding of onset will lead to operations and control techniques that avoid the most harmful form of ELMs. Knowledge of heat transport can be used to engineer configurations that are robust to ELM activity. We have only begun the efforts to understand these processes.
At present, the results must be considered preliminary, since more numerical investigation is needed. The results are sufficiently mature, however, to provide a useful guide for the requirements of nonlinear computations. We expect that a two-fluid Ohms law, FLR effects (gyro-viscosity), and anisotropic thermal diffusivity will become effective stabilization mechanisms at pedestal parameters relevant to the true experimental crossing of the stability boundary, where the linear MHD growth rate is lower. The nonlinear computations presented here are some of the first to show significant plasma-wall interactions as a result of an ELM instability over a global computational domain. The results appear to reproduce some important experimental observations such as the energy loss time-scale and the heat loss mechanism. This suggests that nonlinear fluid based computational models have the potential to provide insight into how ELMs evolve and deposit heat onto the wall. Further refinements will include more accurate geometry, improved boundary conditions, and more refined physics models, such as nonlinear two-fluid and gyro-viscous effects.