Computational Challenges and Optimization in Kinetic Plasma Physics

Plasmas are composed of a large number of charged and neutral particles that interact through electromagnetic forces and collisions over a wide range of temporal and spatial scales to produce collective dynamics. The discrete nature of the plasma particles, their large number, and the wide range of interaction scales make a many-body treatment computationally prohibitive. A statistical treatment evolves the probability density function for each particle species in phase space (x,v) and produces the more manageable but six-dimensional kinetic model, which is typically considered to provide the highest physical fidelity. Solutions to the kinetic model are possible by phase space discretizations using high-order finite element methods that exploit GPU-based HPC systems. Solutions using the discontinuous Galerkin method will be described, where the governing equations are expressed in a balance law form, an approximate Riemann solver is developed for evaluating inter-element numerical fluxes, and Runge-Kutta methods perform the temporal advance. Tailored initializations using kinetic eigenfunctions are shown to access saturated physics regimes with computational efficiency. Assumptions related to thermodynamic relaxation reduce the kinetic model through velocity moments to retain more limited information about the distribution function at each position (x). The moment models, e.g. 5N, 10N, and 13N, are three dimensional and thereby offer substantial computational acceleration compared to the kinetic model. The validity of the assumptions and the resulting moment models depends on local plasma parameters such as collisionality, charge neutrality, and magnetization. Such considerations are necessary when selecting the appropriate model for simulations. Expressing the hierarchy of continuum plasma models in a consistent formulation facilitates hybridization. The computational methods are applied to the GEM collisionless magnetic reconnection problem, 2D electrostatic and electromagnetic kinetic instabilities, and the magnetized Kelvin-Helmholtz instability. Spatially localized deviations from thermodynamic equilibrium are investigated, as well as the resulting agreement and disagreement of global solution metrics between plasma models of differing fidelity. The results suggest opportunities to hybridize by applying the simplest, locally valid plasma model and developing strategies to couple the models across subdomain boundaries.

One of the primary requirements of a fusion power plant is the sufficient confinement of the energetic fusion products. Recently, “precisely quasisymmetric” configurations have been obtained through numerical optimization, demonstrating excellent confinement of the guiding center trajectories of fusion-born alpha particles. There is, however, the potential for enhanced alpha losses due to resonant wave-particle interactions. Modeling the dynamics of energetic particles (EPs) presents challenges due to the multi-scale nature of the problem, coupling the shear Alfven waves of the background plasma with drive due to the energetic particle distribution function. Additional physics is introduced in moving from axisymmetric to 3D fields, as complex orbit types and continuum structures arise. We present pathways to model and optimize the alpha transport driven by Alfvenic instabilities in stellarators. Analysis of EP resonances and their impact on saturation mechanisms indicate key departures from the AE-driven transport in tokamaks, such as the avoidance of phase-space island overlap in quasihelical configurations. Taking inspiration from the condensed matter community, we apply spectral density methods to analyze the shear Alfven continuum gap structure and determine the geometric dependence of gap widths.

Plasma confinement in tokamaks is limited by slow particle and energy losses due to turbulence, driven by unstable drift waves triggered by plasma inhomogeneities. Understanding the mechanisms that drive turbulence in burning plasmas is essential in designing next-generation fusion reactors like ITER with optimal confinement and fusion performance. While ion-scale turbulence in the tokamak core has been modeled extensively, gyrokinetic simulations of turbulent transport in the edge pedestal of high-confinement mode regimes are very challenging. In the pedestal, large gradients in density and temperature can drive multiple drift-wave instabilities across a broad range of scales, namely ion and electron scales that differ in size by nearly two orders of magnitude. Simulating this so-called “multiscale turbulence”, which captures the complex interaction between the slow, large scale motion of the hydrogenic ions and the fast, small-scale motion of the electrons, requires both leadership-scale computing resources and advanced numerical algorithms. Recently, a first analysis of the spectral transition of multiscale turbulence in pedestal-like regimes was demonstrated. The simulations were enabled using advances in numerical algorithms for Eulerian gyrokinetic simulations that use a spatial discretization and array distribution scheme that is spectral/pseudospectral in four of the five phase space dimensions and targets scalability on next-generation, exascale HPC systems that use multicore and GPU-accelerated hardware. Performance optimizations for both compute and communication that target capability-computing simulations on Frontier, the world’s first exascale system, will be introduced. Essential components include highly leveraging GPU offloading and GPU-aware MPI, hipFFT libraries that enhance efficiency for very large numbers of radial/binormal wavenumbers, and adaptive time-stepping that gives fast solution for systems with impulse and oscillatory behavior. * Work supported by US DOE Grants DE-FG02-95ER54309, DE-FC02-06ER54873, DE-SC0017992, and DE-SC0024425.

I will give an overview of recent theoretical and numerical calculations of confined magnetic mirror plasmas with the goal of creating a break even class axisymmetric mirror (Forest, JPP 2023). Fokker Planck solutions from the CQL3D code combined with magnetic equilibria tracking anisotropic pressures and self-consistent ambipolar electric fields provide the starting point as an initial condition for stability analysis. Hybrid VPIC is now being used to model both MHD and kinetic stability and GKYLL fluid and gyro kinetic approaches are not far behind. Major challenges for the mirror and building a digital twin will be discussed.

The study of uncertainty propagation is of fundamental importance in plasma physics simulations. However, the construction of numerical methods is challenging due to the high-dimensionality of the problem, the presence of multiple space-time scales, and the constraints imposed by the need to preserve certain relevant physical properties. In the present talk we present a novel class of methods based on a stochastic-Galerkin (sG) approximation of the uncertain particles’ position and velocity. We first describe the method in the case of the Vlasov-Poisson system with a BGK term describing plasma collisions and than consider the challenging case where collisions are characterized by the Landau-Fokker-Planck operator. We show that the sG particle method preserves the main physical properties of the problem, such as conservations and positivity of the solution, while achieving spectral accuracy for smooth solutions in the random space. Furthermore, in the fluid limit we discuss how to design the sG particle solver in order to possess the asymptotic-preserving property, thus avoiding the loss of hyperbolicity typical of conventional sG methods based on finite differences or finite volumes. References 1. A. Medaglia, L. Pareschi, M. Zanella, Stochastic Galerkin particle methods for kinetic equations of plasmas with uncertainties, J. Comp. Phys. Volume 479, 112011, 2023 2. A. Medaglia, L. Pareschi, M. Zanella, Particle simulation methods for the Landau-Fokker-Planck equation with uncertain data, preprint arXiv:2306.07701, 2023 I will arrive on Monday, February 22, in the afternoon. Therefore, I will be able to parti

In this talk, we introduce two new aspects of inverse problems formulated as PDE-constrained optimization. Firstly, while current approaches assume deterministic parameters, many real-world problems exhibit stochastic behavior. We present a novel approach that treats the PDE solver as a push-forward map to recover the full distribution of unknown random parameters. We introduce a gradient-flow equation to estimate the ground-truth parameter probability distribution. Secondly, as problem dimensions increase, Monte Carlo methods regain relevance. However, directly applying them to gradient-based PDE-constrained optimization poses challenges due to the product of forward and adjoint solutions involving Dirac deltas. We propose strategies to rescue Monte Carlo methods and make them compatible with gradient-based optimization. This talk is based on joint work with Qin Li and Yunan Yang.

Gyrokinetic simulations on parallel supercomputers provide the gold standard for theoretically determining turbulent transport in magnetized fusion plasmas. Applications to large and costly future machines, in particular burning plasma devices like ITER, call for a proper Uncertainty Quantification (UQ) in order to assess the reliability of certain predictions. However, the high computational cost of gyrokinetic simulations prevents straightforward applications of conventional UQ approaches. Here, we present a sensitivity-driven dimension-adaptive sparse grid interpolation approach that can enable UQ in such expensive simulations. This approach explores and exploits the fact that in most problems (i) only a subset of the uncertain parameters are important and (ii) these parameters interact anisotropically. The usefulness of our sensitivity-driven approach will be demonstrated in a realistic description of turbulent transport in the edge of fusion experiments. In a non-linear scenario with more than $264$ million degrees of freedom and eight uncertain inputs, it requires a mere total of $57$ high-fidelity simulations. In addition, we will show that a byproduct of our sensitivity-driven approach is that it also enables the construction of surrogate transport models, which is crucial for tasks such as profile predictions or the design of optimized devices. Lastly, we will discuss the limitations of our approach and some ideas for overcoming them.

Plasmas in the reversed field pinch (RFP), a type of toroidal magnetic confinement device, are characterized by a robust turbulent cascade that physically links unstable tearing modes at global scales to unstable drift waves at microscales. Efforts to model this complex and phenomenologically rich multiscale interaction using gyrokinetic simulation are described. The gyrokinetic code GENE is used, a delta-f solver with global nonlinear simulation capability in toroidal geometry. To properly describe tearing-mode instability, which historically has been analyzed in the resistive MHD framework, a nonuniform equilibrium current profile consistent with discharges in the Madison Symmetric Torus RFP is incorporated in GENE via a shifted Maxwellian background distribution with a radially varying shift. This new capability for GENE has been extensively and successfully tested via benchmarks and comparisons with other gyrokinetic codes, fluid codes, and analytic theory. The benchmarks assess the fidelity of numerical modeling of the tearing-mode instability in terms of growth rates, frequencies, and scalings with key parameters. Comparisons between the new global capability and prior flux-tube calculations have also been made of the instability properties of the dominant trapped-electron-mode (TEM) microinstability in these
discharges. Multiscale simulations have examined the effect of tearing modes on the zonal flows excited by TEM turbulence. The magnetic fluctuations of saturated tearing modes are observed to erode the zonal flows. Then, because zonal flows regulate the saturation level of TEM turbulence, the heat flux driven by TEM turbulence is considerably higher in the presence of the tearing-mode turbulence. This confirms earlier theoretical hypotheses and reduced gyrokinetic models that represented the tearing modes as a external magnetic perturbation. It also establishes gyrokinetics as a model for addressing the wealth of physics in the RFP multiscale turbulent fluctuation
spectrum.

Collaborators: T. Jisuk (University of Wisconsin, Madison) and M.J. Puescel (Dutch Institute for Fundamental Energy Research)

As we progress down the stellarator path to a fusion pilot plant, a key need will be the ability to accurately and robustly predict the equilibrium profiles (and hence fusion power) attainable by a given FPP design configuration. In modern optimized stellarators, neoclassical transport can be reduced to the point that turbulent transport becomes the dominant confinement loss mechanism. Thus a transport solution that includes first-principles turbulence modeling is required. T3D is a framework that leverages multi-scale gyrokinetic theory to efficiently model macro-scale profile evolution in fusion plasmas (tokamaks and stellarators) due to micro-scale turbulent and neoclassical processes. To model micro-turbulence we use the GX gyrokinetic code, which has been developed as a GPU-native code that uses an efficient pseudo-spectral discretization scheme to target fast turbulence calculations for reactor design and optimization. Neoclassical transport is modeled by the KNOSOS neoclassical solver, which uses orbit-averaging to solve the drift kinetic equation very efficiently at low collisionality. Additionally, both GX and KNOSOS use a radially-local approach, taking advantage of the reactor-relevant limit of small orbit widths and small turbulent eddies. This enables a series of GX and KNOSOS calculations to be embedded in parallel in the Trinity3D transport solver for tractable fusion profile prediction (and evolution) calculations. An implicit Newton method is used to time-evolve the transport equations, allowing large timesteps on the timescale of the energy confinement time. Ion temperature profile predictions subject to ion temperature gradient turbulence can be completed in less than an hour, while multi-channel predictions including electron temperature and density profiles can be completed in roughly a day. We will demonstrate a validation of the framework via modeling of W7-X plasmas and discuss future plans for using the framework in experimental studies as well as stellarator FPP design and optimization.

The guiding center and gyrokinetic theory of magnetized particle motion can be extended to the regime of large electric field gradients and strong shear flows. A gradient in the electric field directly modifies the oscillation frequency and causes the Larmor orbits to deform from circular to elliptical trajectories. In order to retain a good adiabatic invariant, there can only be strong dependence on a single coordinate at lowest order, so that resonances do not generate chaotic motion that destroys the invariant. When a single gradient direction is dominant, the guiding center drift velocity becomes anisotropic in response to external forces and additional curvature drifts must be included. The electric polarization and magnetization are also modified by the change in gyrofrequency and by the large drift flows. The theory can be applied to electric fields and shear flows that are even stronger than those observed in the edge transport barrier of a high-performance tokamak (H-mode) pedestal, even if the toroidal field is as small as or even smaller than the poloidal field. Yet, the theory retains a mathematical form that is similar to the standard case and can readily be implemented within existing simulation tools.
The theory above motivates the use of exponential integrators for accelerating the computation of particle trajectories. We have improved exponential integrators for small systems of equations by explicitly computing the eigenfrequencies and using the Lagrange-Sylvester
method for the computation of matrix exponential functions. We show that explicit exponential integrators can provide performance enhancements over the best performing explicit methods, the Boris and Buneman method, when gyrodradius accuracy is required. This opens a promising direction for particle-in-cell codes.

*LLNL-ABS-858862 prepared by LLNL for U.S. DOE under Contract DE-AC52-07NA27344.

Collaborators: T. Alex Nguyen (Lawrence Livermore National Laboratory) and Mayya Tokaman (University of California, Merced)

In plasmas created by laser-matter interactions, the electron thermal and magnetic field transport almost always require a kinetic treatment. We will describe progress in coupling a Vlasov-Fokker-Planck description of the electrons to a radiation-hydrodynamics model. The model has been applied to one-dimensional ablation of solid targets by high power lasers, and can be compared to the usual flux-limited diffusion treatments of electron transport which are a standard feature of hydro models. Some of the computational challenges, including designing algorithms that work across a wide variety of plasma conditions and in the presence of sharp gradients, will be discussed.

The quantized tensor network (QTN) format is a low-rank framework which enables the efficient calculation and manipulation of high-dimensional datasets. For QTNs of rank $D$ representing problems of size $N$, computational costs typically scale like $sim mathcal{O}(text{poly}(D) log(N))$. QTNs have been successfully applied to a wide range of numerical simulation problems such as finding extremal eigenvalues, solving linear systems, and even performing time integration of nonlinear partial differential equations; solving them with reasonable accuracy while significantly reducing computational costs. Here, we discuss the application of QTNs for solving the Vlasov-Maxwell equations, and show results for common test problems in 2D3V. In particular, we find that QTNs of rank $D=64$ appears to be sufficient for capturing the expected physics, despite the simulations using a total of $2^{36}$ grid points and thus requiring $D=2^{18}$ for exact calculations. Additionally, we utilize a QTN time evolution scheme based on the Dirac-Frenkel variational principle, which allows us to use larger time steps than that prescribed by the Courant-Friedrichs-Lewy (CFL) constraint. As such, the QTN format appears to be a promising means of approximately solving the Vlasov-Maxwell equations with significantly reduced cost.

Many plasma systems, from here on Earth in fusion reactors and lab plasma experiments to far away compact objects and galaxy clusters, are magnetized. The underlying dynamic magnetic field imposes some structure to the plasma’s evolution, affecting the transport of momentum and energy parallel vs. perpendicular to the magnetic field. A variety of asymptotic models have been derived utilizing the importance of the plasma’s magnetic field for its dynamics. Many of these models are highly non-trivial and their robust, accurate numerical discretization an even grander challenge.

I will present an overview of the Gkeyll simulation framework, which has implemented a number of models for modeling magnetized plasmas for diverse applications, from fusion to astrophysical applications. I will discuss subtleties and developments to create energy-conserving schemes using the discontinuous Galerkin finite element method. I will focus in particular on a novel parallel-kinetic-perpendicular-moment model, which performs a spectral expansion of only the perpendicular degrees of freedom, analogous to spectral methods which have grown in popularity in recent years for gyrokinetics, while retaining the complete dynamics parallel to the magnetic field. Finally, I will demonstrate the utility of Gkeyll’s implementation of these models on a variety of benchmarks and production problems. 

Magnetic confinement fusion reactors feature numerous disparate time-scales, with the gyrofrequency induced by the strong magnetic field among the fastest.  While gyrokinetic models permit stepping over this scale, there is increasing concern that the gyrokinetic ordering may break down in certain key regions of such reactors.  This generates interest in asymptotic-preserving schemes that can recover the gyrokinetic limit when it is valid while still converging to full-orbit dynamics in the limit of small time-step.  We present such a scheme which, in addition to the desired asymptotic properties, preserves the exact energy and charge conservation enjoyed by implicit particle-in-cell (PIC) methods developed over the last decade.  We discuss the scheme’s theoretical development, an efficient and robust nonlinear solver, and numerical results for both single-particle test cases and self-consistent PIC simulations.   * Prepared by LLNL under Contract DE-AC52-07NA27344.