Here, we report on a 3D MHD simulation of multiple sawteeth occurring in a tokamak plasma with moderate. We employ a modern, massively parallel, implicit 3D MHD code, which uses high-order finite elements in all three dimensions.These features enable high-resolution and long timescale calculations of MHD activity in a tokamak. Unlike the papers cited in Sec. II A , we do not start the configuration in an unstable state and watch its transition into a stable state. Rather, we define particle and energy sources and transport coefficients and run forofof Alfvén times and look for a repeating cycle.

A. Sawtoothing simulation

Figs. 4–6 R / a = 3.2 , ellipticity κ = 1.3 , triangularity δ = 0.2 , β ≳ 2 % , and edge safety factor q a = 4.3 . Figure 4 q 0 never falls below unity, and the oscillations and fast crashes are not primarily due to the (1, 1) mode. At the start of each temperature crash, a large number of localized (m, n) modes with m = n > 1 grow up and cause only the central region to become stochastic. We show the results inof a M3D-C1 simulation of a canonical tokamak discharge with aspect ratio, ellipticity, triangularity, and edge safety factorshows that the simulation develops quasi-periodic oscillations in which the central temperature slowly rises and abruptly crashes, as is the case in sawtooth oscillations. Consistent with the new model described above, in this simulation,never falls below unity, and the oscillations and fast crashes are not primarily due to the (1, 1) mode. At the start of each temperature crash, a large number of localized () modes withgrow up and cause only the central region to become stochastic.

Fig. 5 Fig. 4 q-profiles at the same two times, which are seen to stay essentially unchanged. Also shown in the two graphs, in dashed lines, are the results we obtained in a 2D axisymmetric calculation with the same transport coefficients and heating sources. Of course, there is no (1, 1) mode activity and no crashes in 2D, and the central temperatures are higher and the q-profile is lower because of this. Shown in(left) is the midplane temperature profile just before and just after a crash at the times indicated by the vertical lines in. The central electron temperature is seen to have decreased by about 25% during the crash. Not shown is that the central density also decreased but only slightly, presumably because it is originally less peaked than the temperature due to the lack of central fueling and because the density equation does not contain a large parallel diffusion term. Shown on the right are the-profiles at the same two times, which are seen to stay essentially unchanged. Also shown in the two graphs, in dashed lines, are the results we obtained in a 2D axisymmetric calculation with the same transport coefficients and heating sources. Of course, there is no (1, 1) mode activity and no crashes in 2D, and the central temperatures are higher and the-profile is lower because of this.

Fig. 6 E × B convective velocities from the many unstable modes and parallel transport in the stochastic region contribute to the temperature flattening. We show Poincaré plots at the same two time slices in. It is seen that before the crash, the magnetic surfaces are mostly good everywhere, with some small islands at rational surfaces. Just after the crash, most of the surfaces are still good, but those near the center have been destroyed. Both theconvective velocities from the many unstable modes and parallel transport in the stochastic region contribute to the temperature flattening.

Figure 7 Fig. 4 T e , as a function of time. The bottom curve shows the kinetic energy in each of the first eight (8) toroidal harmonics during this time period. It is seen that as multiple modes become unstable and their kinetic energy peaks, the maximum T e begins to decrease. For the event shown, the n = 4 mode reaches the largest amplitude, and n = 2 persists, seemingly in a state of marginal stability. This is also the case for other crashes in this sequence and may be a consequence of the elongation. is a closeup of one of the “crash” periods, delineated with dashed vertical lines in. The top frame shows the maximum electron temperature,, as a function of time. The bottom curve shows the kinetic energy in each of the first eight (8) toroidal harmonics during this time period. It is seen that as multiple modes become unstable and their kinetic energy peaks, the maximumbegins to decrease. For the event shown, the4 mode reaches the largest amplitude, and n = 2 persists, seemingly in a state of marginal stability. This is also the case for other crashes in this sequence and may be a consequence of the elongation.

5,39 5, 014002 (2012). 5. S. C. Jardin, N. Ferraro, J. Breslau, and J. Chen, Comput. Sci. Discovery, 014002 (2012). https://doi.org/10.1088/1749-4699/5/1/014002 39. N. M. Ferraro, S. Jardin, M. Shephard, A. Bauer, J. Breslau, J. Chen, F. Delalondre, X. Luo, and F. Zhang, “ Fluid modeling of fusion plasmas with M3D-C1 ,” in Proceedings of the SciDAC 2011 , Devner, CO, 10–14 July (2011). η ∼ T e − 3 / 2 , but the resistivity was uniformly enhanced so that the central Lundquist number was S = 10 6 . It used a uniform viscosity with a value 10 times the central resistivity (dimensionless code units). The perpendicular thermal conductivity varied, increasing radially from 18 (center) to 36 (edge) times the central resistivity. The parallel thermal conductivity was 106 times greater than the perpendicular one. Sufficient beam heating was applied to maintain β at about 2.5 % . The neutral beam model also drove a centrally peaked sheared toroidal velocity, which never exceeded 10% of the Alfvén velocity. A loop voltage was applied at the boundary in a feedback loop to keep the total toroidal current constant in time. The code uses a 3D finite element mesh. This calculation had 32 Hermite cubic finite elements in the toroidal direction and an unstructured mesh in the poloidal plane with 4th order Bell elements 40 40. D. Braess, Finite Elements ( Cambridge University Press , Cambridge , 2002). 5 × 10 5 τ A requiring ≈ 10 6 processor-hours of computer time using 2.4 GHz Xeon processors with Infiniband interconnect. The nonlinear calculations shown here were performed with the M3D-C1 code.A time and space varying Spitzer resistivity profile was used,, but the resistivity was uniformly enhanced so that the central Lundquist number was. It used a uniform viscosity with a value 10 times the central resistivity (dimensionless code units). The perpendicular thermal conductivity varied, increasing radially from 18 (center) to 36 (edge) times the central resistivity. The parallel thermal conductivity was 10times greater than the perpendicular one. Sufficient beam heating was applied to maintainat about. The neutral beam model also drove a centrally peaked sheared toroidal velocity, which never exceeded 10% of the Alfvén velocity. A loop voltage was applied at the boundary in a feedback loop to keep the total toroidal current constant in time. The code uses a 3D finite element mesh. This calculation had 32 Hermite cubic finite elements in the toroidal direction and an unstructured mesh in the poloidal plane with 4th order Bell elementsof typical linear size 6 cm. This was a single-fluid simulation in which the temperature, density, and all components of the magnetic field were advanced. The calculation ran forrequiringprocessor-hours of computer time using 2.4 GHz Xeon processors with Infiniband interconnect.

Fig. 1 ( q 0 , β p 1 ) space of a sawtooth oscillation as calculated and depicted in Figs. 4–6 m, n) modes with m = n > 1 is exceeded. As these modes grow, they locally increase the pressure gradient and excite other (m, n) modes with m = n > 1 . These many interchange modes destroy the surfaces in the center as shown on the right in Fig. 6 q-profile as shown in Fig. 5 q 0 will initially decrease due to resistive diffusion. Once the (1, 1) stability boundary is crossed (C), the associated central dynamo voltage will act to raise q 0 and stabilize its drop. When the increasing central pressure causes the stability threshold for several (m, n) modes with m = n > 1 to be crossed again, the process will repeat. Also shown inis a schematic trajectory in thespace of a sawtooth oscillation as calculated and depicted in. At location (A), the instability threshold for several () modes withis exceeded. As these modes grow, they locally increase the pressure gradient and excite other () modes with. These many interchange modes destroy the surfaces in the center as shown on the right in, collapsing the central pressure without changing the-profile as shown in. Once the central pressure is flattened (B), these modes become stable and the magnetic surfaces re-form. As central heating is applied, the central pressure will again increase andwill initially decrease due to resistive diffusion. Once the (1, 1) stability boundary is crossed (C), the associated central dynamo voltage will act to raiseand stabilize its drop. When the increasing central pressure causes the stability threshold for several () modes withto be crossed again, the process will repeat.