{'fraction_non_alphanumeric': 0.02357857036263244, 'fraction_numerical': 0.009998507685420087, 'mean_word_length': 4.847731239092496, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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}}
Massively parallel modeling of the sawtooth instability in tokamaks

The sawtooth instability is one of the most fundamental dynamics of an inductive tokamak discharge such as will occur in ITER. The sawtooth occurs when the current peaks in a tokamak, creating a region in the core where the safety factor is less than unity, q<1. While this instability is confined to the center of the plasma in low-pressure, low-current, large-aspect-ratio discharges, under certain conditions it can create magnetic islands at the outer resonant surfaces or set off a sequence of events that leads to a major disruption. Sawtooth behavior is complex and remains incompletely explained. The SciDAC Center for Extended MHD Modeling (CEMM) has undertaken an ambitious campaign to model this periodic motion as accurately as possible using the most complete fluid-like description of the plasma - the Extended MHD model. The multiple time and space scales associated with the reconnection layer and growth time make this an extremely challenging computational problem. The most recent simulation by the M3D code used over 500,000 elements for 400,000 partially implicit time steps for a total of 2×1011 space-time points, and there still remain some resolution issues. However, these calculations are providing insight into the nonlinear mechanisms of surface breakup and healing. We have been able to match many features of a small tokamak and can now project to the computational requirements for simulations of larger, hotter devices such as ITER.


Introduction
The resistive internal kink instability in tokamaks leads to the phenomenon of the sawtooth crash [1], which mixes plasma from the core and outer regions, cooling the center. It is therefore desirable to develop a model of the triggering and consequences of such events that provides a quantitative predictive capability so that they can be avoided or controlled. Analytic theories and quasi-empirical theory-based models such as that of Porcelli [2] provide a good starting point, but more fundamental numerical models stand a better chance of giving accurate numbers for particular discharges in particular physical devices and for testing and optimizing various control techniques.
The M3D (Multilevel 3D) code [3] was designed for the study of such macroscopic instabilities. It employs any of a hierarchy of fluid and fluid-particle hybrid models in tracking the nonlinear evolution of magnetically confined plasmas in toroidal topologies. The code uses a partially implicit time integration algorithm in which the fast wave and diffusive terms are treated implicitly but the shear Alfvén, slow wave, and convective terms are all explicit. For the purposes of the initial study described here, the resistive MHD model was used. Future studies will employ a more complete Extended MHD model.
It is well-known that resistive MHD phenomena involving magnetic reconnection, such as the sawtooth crash, are associated with thin sheets of high current density whose thickness scales as S where the Lundquist number S is the ratio of the resistive diffusion time to the Alfvén transit time in the plasma. Because S is roughly proportional to temperature T 3/2 , its value in large tokamak experiments is typically of the order of 10 8 or higher, which would result in current sheets too thin to be practically resolved in a calculation that makes use of a partially explicit time discretization and is hence limited in step size by the CFL condition. Such high-temperature discharges also require the inclusion of two-fluid or kinetic effects in addition to MHD. In contrast, smaller, colder tokamaks with S of 10 5 or lower may have current sheets thick enough to be resolved in a long calculation, and so are better candidates for realistic MHD physics simulation on present-day computers. The resistive MHD model gives a reasonable description of magnetic reconnection in these devices, although there will be some quantitative differences when the more complete extended MHD model is used. The Current Drive Experiment Upgrade (CDX-U) [4] is a small, low-aspect ratio university-scale tokamak at the Princeton Plasma Physics Laboratory. While its discharges are much colder and briefer than those in larger experiments, it is reasonably well-diagnosed, and exhibits the same MHD instabilities that are of interest in reactor-relevant devices. In particular, it has a sawtooth regime with an easily characterizable period of approximately 125 µs. With S≈10 4 , it is therefore an ideal candidate for modeling with M3D.
In Section 2, we set out the choices of physical parameters and details of the numerical model for the CDX-U sawtooth simulation. In Section 3, results of a linear study are presented. Section 4 describes the results of the nonlinear study, including attempts at convergence in mesh resolution. Section 5 contains our conclusions.

Statement of the problem
The study begins with a sequence of equilibria that occur as time slices from a reconstruction of a typical CDX-U discharge using the 2D (axisymmetric) transport-timescale code TSC [5]. Stability to sawteeth depends sensitively on the instantaneous profile of the "safety factor" q(r)~B T r/B p R, where B T and B p are respectively the toroidal and poloidal components of the local magnetic field, and r and R are the minor and major radii. The central value of q is inversely proportional to the current density on axis, q(0) ∝ J(0) -1 . During the course of the discharge, J(0) increases sufficiently to cause q(0) to fall below unity, the condition for instability of the sawtooth mode. The reference case for the simulations described in the rest of this paper is one for which q(0) has dropped to 0.92, increasing monotonically outward, with the q=1 surface approximately a third of the way out from the axis. A description of this equilibrium state is then imported into M3D as an initial condition.
The magnetic field and the pressure and density profiles are read into M3D from the equilibrium file and used to initialize both the fixed-boundary geometry (with a flux-surface-aligned unstructured triangular mesh) and the fields themselves. The initial resistivity is taken to be the Spitzer value based on the initial temperature profile, which gives a peak S of 1.94×10 4 , with the minimum value constrained to be no smaller than 1.94×10 2 . The resistivity is not evolved in time. The viscosity is set to ten times the resistivity at the peak location and is constant in space and time. The unstructured mesh used 79 zones in the radial direction, giving it a resolution of approximately 9200 vertices in each poloidal cross-section. Because the initial equilibrium is axisymmetric, linear modes with different toroidal mode number n decouple, although each contains a wide range of poloidal modes 0 ≤ m ≤ m max represented on the finite element grid.
The first objective was to run the M3D code in its linear mode in order to determine the eigenmodes and linear growth rates of the starting equilibrium for the first few toroidal mode numbers, n=1-7. The dependencies of these modes and growth rates on the choice of transport coefficients were also investigated. The primary studies are full nonlinear runs of the code, initialized by superimposing the n=1 eigenmode found in the previous study on the equilibrium state at an amplitude such that the maximum B p in the perturbation is 10 -4 of the maximum B T in the equilibrium. At various mesh resolutions, the plasma is allowed to evolve in the presence of volumetric source terms that tend to restore the current density and temperature to their equilibrium values. The two highest-n resolvable toroidal modes are filtered out to avoid aliasing problems; the others are kept and couple to each other in a fully nonlinear way. The behavior of the current and temperature profiles and of the field topology are closely monitored for indications of m=1 activity as the equations are solved out to several predicted sawtooth periods.

Linear results
The linear search for an unstable n=1 eigenmode of the provided equilibrium revealed that such a mode does exist ( figure 1a). It has a predominantly m=1 character centered at the q=1 rational surface, consistent with the sawtooth model, and a normalized linear growth rate of γτ A ≈ 8.61×10 -3 . The search for unstable higher-n modes unexpectedly yielded positive results as well. As shown in figures 1b,c, they tend to occur near the boundary. These modes have growth rates higher than that of the n=1 mode, said rates increasing with n and approximately in proportion to η 3/5 . This scaling and the location of the modes both suggest that these should be identified as resistive ballooning modes [6]. Their growth rates are reduced by the addition of a realistic level of parallel heat conduction, but can be linearly stabilized only by the additional imposition of a large isotropic heat conduction: approximately 200 m 2 /s, consistent with the boundary value in the TSC reconstruction. It is conjectured that the self-consistent thermal transport associated with these modes nonlinearly stabilizes them at a small amplitude in the experiment [7]. a.

Sawtooth cycle
To initialize the nonlinear study, the n=1 eigenmode described in section 3 was added as a small initial perturbation to the unstable equilibrium. Both parallel and perpendicular heat conduction were employed as described above, stabilizing the higher-n modes. As shown in figure 2a, this nonlinear run begins with a long phase of essentially linear behavior during which the n=1 mode grows at the linearly predicted rate, while nonlinearly coupling to and destabilizing higher n modes. These higher-n modes have predominantly m=n character and are at the q=1 surface, and so should not be confused with the higher-n eigenmodes seen near the boundary in the linear study. a.
b. The sequence of Poincaré plots in figure 3 shows the evolution of the magnetic flux surfaces during a typical sawtooth cycle. The nonlinear m=1, n=1 mode results in a complete sawtooth crash, with formation of a new magnetic axis, and reconnection of the original one, as in the Kadomtsev model [8]. Because the n=1 mode and the nonlinearly driven higher-n modes also have higher m components at rational surfaces with q>1, the (1,1) island merging coincides with the growth and overlap of other island chains at larger minor radii, resulting in general stochasticity after the "crash". As a current source term gradually drives the central safety factor back below unity (and a temperature source reheats the cooled core region), the stochastic regions heal and again become good flux surfaces. The process then repeats. The cycle time of 395 τ A is equivalent to about 100 µs, in reasonably good agreement with observations. a.

Convergence study
In order to assess numerical convergence, the nonlinear study was repeated with the same poloidal resolution but with twice the toroidal resolution, allowing modes up to n=22 to be resolved. The behavior of the higher-resolution run, however, was qualitatively different from that of the initial one. The initial crash happened sooner and with greater violence, and subsequent crashes followed with a shorter period but with diminishing energy (figure 2b). This failure to converge with more planes appears to be caused by the presence of large amounts of energy in the additional modes. We are performing additional studies to determine whether these short-wavelength modes are physical or are artifacts of converging in only one dimension or of the simplified resistive MHD model.

Conclusion
In a series of high-resolution, massively parallel nonlinear MHD runs totaling approximately 100,000 CPU-hours, the sawtooth cycle in the CDX tokamak has been successfully tracked with the M3D code. The lower-resolution case showed good agreement with the observed sawtooth period (though not with the observed crash time), but the higher-resolution check failed to agree with this period. It is likely that still higher resolutions, as well as greater fidelity to the physics of the actual experiment, will be required to produce good quantitative agreement. Areas for improvement include the use of extended (two-fluid) MHD terms which are likely to give faster reconnection and to preferentially damp the higher-n modes; and the use of ohmic heating and current drive rather than the volume source terms now in use for these quantities, which do not have direct analogues in the device. Better transport models, including an evolving resistivity profile and a perpendicular heat transport profile with much smaller interior values, will also be required.
The results so far are encouraging, but the task of scaling such studies up to ITER-scale devices remains a formidable challenge. The vastly higher temperatures and, to a lesser extent, the much larger physical size, will necessitate an increase of several orders of magnitude in spatial resolution, which, due to the partially explicit nature of the code, will also require orders of magnitude more time steps to reach the same physical time. In order to achieve this, we will not be able to rely solely on increases in hardware capability, but must also improve the software, making use of techniques such as adaptive mesh refinement and making the codes more efficient, more scalable, and more implicit.