In fusion reactor design, steels under consideration for the blanket are ferromagnetic, so the steel's effect on the plasma physics must be examined. For efficient calculation of these fields, we can exploit the fact that the ferromagnetic material gives a small perturbation relative to the fields from the electromagnetic coils and plasma. Moreover the magnetization is saturated due to the strong fields in typical fusion systems. These approximations significantly reduce the nonlinearity of the problem, so the magnetic materials can be described by an array of point dipoles of known magnitude, oriented in the direction of the coil and plasma field. The approach is verified by comparison to finite-element calculations with commercial software and shown to be accurate. As no linear or nonlinear solve is required, only evaluation of Biot-Savart-type integrals, the method here is significantly simpler to implement than other methods, and extremely fast. The method is compatible with arbitrary CAD geometry, and also allows rapid computation of the magnetic forces. We demonstrate adding the ferromagnetic effects to free-boundary magnetohydrodynamic equilibrium calculations, assessing the effect on plasma physics properties such as confinement and stability. Moreover, it is straightforward to differentiate through the model to get the derivative of the field with respect to the electromagnet parameters. We thereby demonstrate gradient-based coil optimization for a quasi-isodynamic stellarator in which the field contribution from a ferromagnetic blanket is included. Even a significant steel volume is found to have little impact on the plasma physics properties, with the main effects being a slight destabilization of ballooning modes and a radial shift of the edge islands due to decrease in rotational transform. Both of these issues are corrected by the minor reoptimization of the coil shapes to account for the field from the steel.

The magnetohydrodynamic (MHD) equilibrium and stability properties of the Infinity Two fusion pilot plant baseline plasma physics design are presented. The configuration is a four-field period, aspect ratio $A = 10$ quasi-isodynamic stellarator optimised for excellent confinement at elevated density and high magnetic field $B = 9\,T$ . Magnetic surfaces exist in the plasma core in vacuum and retain good equilibrium surface integrity from vacuum to an operational $\beta = 1.6 \,\%$ , the ratio of the volume average of the plasma and magnetic pressures, corresponding to $800\ \textrm{MW}$ deuterium-tritium fusion operation. Neoclassical calculations show that a self-consistent bootstrap current of the order of ${\sim} 1\ \textrm{kA}$ slightly increases the rotational transform profile by less than 0.001. The configuration has a magnetic well across its entire radius. From vacuum to the operating point, the configuration exhibits good ballooning stability characteristics, exhibits good Mercier stability across most of its minor radius and it is stable against global low-n MHD instabilities up to $\beta = 3.2\,\%$ .

We provide an assessment of the Infinity Two fusion pilot plant (FPP) baseline plasma physics design. Infinity Two is a four-field period, aspect ratio $A = 10$ , quasi-isodynamic stellarator with improved confinement appealing to a max- $J$ approach, elevated plasma density and high magnetic fields ( $ \langle B\rangle = 9$ T). Here $J$ denotes the second adiabatic invariant. At the envisioned operating point ( $800$ MW deuterium-tritium (DT) fusion), the configuration has robust magnetic surfaces based on magnetohydrodynamic (MHD) equilibrium calculations and is stable to both local and global MHD instabilities. The configuration has excellent confinement properties with small neoclassical transport and low bootstrap current ( $|I_{bootstrap}| \sim 2$ kA). Calculations of collisional alpha-particle confinement in a DT FPP scenario show small energy losses to the first wall ( ${\lt}1.5 \,\%$ ) and stable energetic particle/Alfv & eacute;n eigenmodes at high ion density. Low turbulent transport is produced using a combination of density profile control consistent with pellet fueling and reduced stiffness to turbulent transport via three-dimensional shaping. Transport simulations with the T3D-GX-SFINCS code suite with self-consistent turbulent and neoclassical transport predict that the DT fusion power $P_{{fus}}=800$ MW operating point is attainable with high fusion gain ( $Q=40$ ) at volume-averaged electron densities $n_e\approx 2 \times 10<^>{20}$ m $<^>{-3}$ , below the Sudo density limit. Additional transport calculations show that an ignited ( $Q=\infty$ ) solution is available at slightly higher density ( $2.2 \times 10<^>{20}$ m $<^>{-3}$ ) with $P_{{fus}}=1.5$ GW. The magnetic configuration is defined by a magnetic coil set with sufficient room for an island divertor, shielding and blanket solutions with tritium breeding ratios (TBR) above unity. An optimistic estimate for the gas-cooled solid breeder designed helium-cooled pebble bed is TBR $\sim 1.3$ . Infinity Two satisfies the physics requirements of a stellarator fusion pilot plant.

Transport characteristics and predicted confinement are shown for the Infinity Two fusion pilot plant baseline plasma physics design, a high field stellarator concept developed using modern optimization techniques. Transport predictions are made using high-fidelity nonlinear gyrokinetic turbulence simulations along with drift kinetic neoclassical simulations. A pellet-fuelled scenario is proposed that enables supporting an edge density gradient to substantially reduce ion temperature gradient turbulence. Trapped electron mode turbulence is minimized through the quasi-isodynamic configuration that has been optimized with maximum-J. A baseline operating point with deuterium-tritium fusion power of $P_{{fus,DT}}=800$ MW with high fusion gain $Q_{{fus}}=40$ is demonstrated, respecting the Sudo density limit and magnetohydrodynamic stability limits. Additional higher power operating points are also predicted, including a fully ignited ( $Q_{{fus}}=\infty$ ) case with $P_{{fus,DT}}=1.5$ GW. Pellet ablation calculations indicate it is plausible to fuel and sustain the desired density profile. Impurity transport calculations indicate that turbulent fluxes dominate neoclassical fluxes deep into the core, and it is predicted that impurity peaking will be smaller than assumed in the transport simulations. A path to access the large radiation fraction needed to satisfy exhaust requirements while sustaining core performance is also discussed.

Vico et al. suggest a fast algorithm for computing volume potentials, beneficial to fields with problems requiring the solution of the free-space Poisson's equation, such as beam and plasma physics. Currently, the standard is the algorithm of Hockney and Eastwood, with second order in convergence at best. The algorithm proposed by Vico et al. converges spectrally for sufficiently smooth functions, i.e., faster than any fixed order in the number of grid points. We implement a performance portable version of the traditional Hockney-Eastwood and the novel Vico-Greengard Poisson solver as part of the Independent Parallel Particle Layer (IPPL) library. For sufficiently smooth source functions, the Vico-Greengard algorithm achieves higher accuracy than the Hockney-Eastwood method with the same grid size, reducing the computational demands of high-resolution simulations since one could use coarser grids to achieve them. Additionally, we propose an improvement to the Vico-Greengard method which further reduces its memory footprint. This is important for GPUs, which have limited memory, and should be taken into account when selecting numerical algorithms for performance portable codes. Finally, we showcase performance through GPU and CPU scaling studies on the Perlmutter (NERSC) supercomputer, with efficiencies staying above 50% in the strong scaling case. To showcase portability, we also run the scaling studies on the Alps supercomputer at CSCS, Switzerland and the GPU partition of the Lumi supercomputer at CSC, Finland.

The multi-region relaxed magnetohydrodynamics (MRxMHD) model partitions plasma into constant-pressure regions separated by ideal MHD interfaces. In each region, magnetic fields satisfy the Beltrami equation with constant parameters. This work extends the boundary integral equation solver for Taylor states (BIEST) to compute stepped-pressure equilibria in axisymmetric geometries. The solution involves iterative Taylor state computations and interface position updates using a fixed-point algorithm enhanced by adaptive step-length control and Anderson acceleration. High accuracy in Taylor state computations after interface updates, along with minimization of total force balance conditions at the interfaces, is achieved through reparametrization and filtering techniques. This approach establishes a robust framework for extending the BIEST solver to axisymmetric MRxMHD equilibria.

In this work, we present a detailed assessment of fusion-born alpha-particle confinement, their wall loads and stability of Alfv & eacute;n eigenmodes driven by these energetic particles in the Infinity Two Fusion Pilot Plant baseline plasma design, a four-field-period quasi-isodynamic stellarator to operate in deuterium-tritium fusion conditions. Using the Monte Carlo codes, SIMPLE, ASCOT5 and KORC-T, we study the collisionless and collisional dynamics of guiding-centre and full-orbit alpha-particles in the core plasma. We find that core energy losses to the wall are less than 4 %. Our simulations shows that peak power loads on the wall of this configuration are approximately 2.5 MW m- $<^>2$ and are spatially localised, toroidally and poloidaly, in the vicinity of x-points of the magnetic island chain $n/m = 4/5$ outside the plasma volume. Also, an exploratory analysis using various simplified walls shows that shaping and distance of the wall from the plasma volume can help reduce peak power loads. Our stability assessment of Alfv & eacute;n eigenmodes using the STELLGAP and FAR3d codes shows the absence of unstable modes driven by alpha-particles in Infinity Two due to the relatively low alpha-particle beta at the envisioned 800 MW operating scenario.

In a recent publication (Toler et al. 2023), we demonstrated that for axisymmetric geometries, the Kapur-Rokhlin quadrature rule provided an efficient and high-order accurate method for computing the normal component, on the plasma surface, of the magnetic field due to the toroidal current flowing in the plasma, via the virtual-casing principle. The calculation was indirect, as it required the prior computation of the magnetic vector potential from the virtual-casing principle, followed by the computation of its tangential derivative by Fourier differentiation, in order to obtain the normal component of the magnetic field. Our approach did not provide the other components of the virtual-casing magnetic field. In this letter, we show that a more direct and more general approach is available for the computation of the virtual-casing magnetic field. The Kapur-Rokhlin quadrature rule accurately calculates the principal value integrals in the expression for all the components of the magnetic field on the plasma boundary, and the numerical error converges at a rate nearly as high as the indirect method we presented previously.

We introduce a novel particle-in-Fourier (PIF) scheme based on [1], [2] that extends its applicability to non-periodic boundary conditions. Our method handles free space boundary conditions by replacing the Fourier Laplacian operator in PIF with a mollified Green's function as first introduced by Vico-Greengard-Ferrando [3]. This modification yields highly accurate free space solutions to the Vlasov-Poisson system, while still maintaining energy conservation up to an error bounded by the time step size. We also explain how to extend our scheme to arbitrary Dirichlet boundary conditions via standard potential theory, which we illustrate in detail for Dirichlet boundary conditions on a circular boundary. We support our approach with proof-of-concept numerical results from two-dimensional plasma test cases to demonstrate the accuracy, efficiency, and conservation properties of the scheme. By avoiding grid heating and finite grid instability we are able to show an order of magnitude speedup compared to the standard PIC scheme for a long time integration cyclotron simulation.

We perform a scaling and performance portability study of the particle-in-cell scheme for plasma physics applications through a set of mini-apps we name "Alpine", which can make use of exascale computing capabilities. The mini-apps are based on Independent Parallel Particle Layer, a framework that is designed around performance portable and dimension independent particles and fields. We benchmark the simulations with varying parameters such as grid resolutions (512^3 to 2048^3) and number of simulation particles (10^9 to 10^11) with the following mini-apps: weak and strong Landau damping, bump-on-tail and two-stream instabilities, and the dynamics of an electron bunch in a charge-neutral Penning trap. We show strong and weak scaling and analyze the performance of different components on several pre-exascale architectures such as Piz-Daint, Cori, Summit and Perlmutter. While the scaling and portability study helps identify the performance critical components of the particle-in-cell scheme in the current state-of-the-art computing architectures, the mini-apps by themselves can be used to develop new algorithms and optimize their high performance implementations targeting exascale architectures.

This paper presents the new 2D electrostatic particle-in-cell code FENNECS developed to study the formation of magnetized non-neutral plasmas in geometries with azimuthal symmetry. This code has been developed in the domain of gyrotron electron gun design, but solves general equations and can be applied in other domains of plasma physics. FENNECS is capable of simulating electron-neutral collisions using a Monte Carlo approach and considers both elastic and inelastic (ionization) processes. It is also capable of solving the Poisson equation on domains with arbitrary geometries with either Dirichlet or natural boundary conditions. The Poisson solver is based on a meshless Finite Element Method, called web-splines, based on b-splines of any order, and used for the first time in the domain of plasma physics. In addition, the effect of fast ions colliding with the electrodes and causing ion induced electron emission at the electrode surfaces has been implemented in the code. In this paper, the governing equations solved by FENNECS and the numerical methods used to solve them are presented. A number of verification cases are then reported. Finally, the parallelization scheme used in FENNECS and its parallel scalability are presented.

In a recent publication (Toler et al., J. Plasma Phys., vol. 89, issue 2, 2023, p. 905890210), we demonstrated that for axisymmetric geometries, the Kapur–Rokhlin quadrature rule provided an efficient and high-order accurate method for computing the normal component, on the plasma surface, of the magnetic field due to the toroidal current flowing in the plasma, via the virtual-casing principle. The calculation was indirect, as it required the prior computation of the magnetic vector potential from the virtual-casing principle, followed by the computation of its tangential derivative by Fourier differentiation, to obtain the normal component of the magnetic field. Our approach did not provide the other components of the virtual-casing magnetic field. In this letter, we show that a more direct and more general approach is available for the computation of the virtual-casing magnetic field. The Kapur–Rokhlin quadrature rule accurately calculates the principal value integrals in the expression for all the components of the magnetic field on the plasma boundary, and the numerical error converges at a rate nearly as high as the indirect method we presented previously.

Many applications in magnetic confinement fusion require the efficient calculation of surface integrals with singular integrands. The singularity subtraction approaches typically used to handle such singularities are complicated to implement and low-order accurate. In contrast, we demonstrate that the Kapur–Rokhlin quadrature scheme is well-suited for the logarithmically singular integrals encountered for a toroidally axisymmetric confinement system, is easy to implement and is high-order accurate. As an illustration, we show how to apply this quadrature scheme for the efficient and accurate calculation of the normal component of the magnetic field due to the plasma current on the plasma boundary, via the virtual-casing principle.

We present a robust optimization algorithm for the design of electromagnetic coils that generate vacuum magnetic fields with nested flux surfaces and precise quasi-symmetry. The method is based on a bilevel optimization problem, where the outer coil optimization is constrained by a set of inner least-squares optimization problems whose solutions describe magnetic surfaces. The outer optimization objective targets coils that generate a field with nested magnetic surfaces and good quasi-symmetry. The inner optimization problems identify magnetic surfaces when they exist, and approximate surfaces in the presence of magnetic islands or chaos. We show that this formulation can be used to heal islands and chaos, thus producing coils that result in magnetic fields with precise quasi-symmetry. We show that the method can be initialized with coils from the traditional two stage coil design process, as well as coils from a near axis expansion optimization. We present a numerical example where island chains are healed and quasi-symmetry is optimized up to surfaces with aspect ratio 6. Another numerical example illustrates that the aspect ratio of nested flux surfaces with optimized quasi-symmetry can be decreased from 6 to approximately 4. The last example shows that our approach is robust and a cold-start using coils from a near-axis expansion optimization.

This work introduces meta estimators that combine multiple multifidelity techniques based on control variates, importance sampling, and information reuse to yield a quasi-multiplicative amount of variance reduction. The proposed meta estimators are particularly efficient within outer-loop applications when the input distribution of the uncertainties changes during the outer loop, which is often the case in reliability-based design and shape optimization. We derive asymptotic bounds of the variance reduction of the meta estimators in the limit of convergence of the outer-loop results. We demonstrate the meta estimators, using data-driven surrogate models and biasing densities, on a design problem under uncertainty motivated by magnetic confinement fusion, namely the optimization of stellarator coil designs to maximize the estimated confinement of energetic particles. The meta estimators outperform all of their constituent variance reduction techniques alone, ultimately yielding two orders of magnitude speedup compared to standard Monte Carlo estimation at the same computational budget.

We extend the single-stage stellarator coil design approach for quasi-symmetry on axis from (Giuliani et al 2020) to additionally take into account coil manufacturing errors. By modeling coil errors independently from the coil discretization, we have the flexibility to consider realistic forms of coil errors. The corresponding stochastic optimization problems are formulated using a risk-neutral approach and risk-averse approaches. We present an efficient, gradient-based descent algorithm which relies on analytical derivatives to solve these problems. In a comprehensive numerical study, we compare the coil designs resulting from deterministic and risk-neutral stochastic optimization and find that the risk-neutral formulation results in more robust configurations and reduces the number of local minima of the optimization problem. We also compare deterministic and risk-neutral approaches in terms of quasi-symmetry on and away from the magnetic axis, and in terms of the confinement of particles released close to the axis. Finally, we show that for the optimization problems we consider, a risk-averse objective using the conditional value-at-risk leads to results which are similar to the risk-neutral objective.

We propose a new method to compute magnetic surfaces that are parametrized in Boozer coordinates for vacuum magnetic fields. We also propose a measure for quasisymmetry on the computed surfaces and use it to design coils that generate a magnetic field that is quasisymmetric on those surfaces. The rotational transform of the field and complexity measures for the coils are also controlled in the design problem. Using an adjoint approach, we are able to obtain analytic derivatives for this optimization problem, yielding an efficient gradient-based algorithm. Starting from an initial coil set that presents nested magnetic surfaces for a large fraction of the volume, our method converges rapidly to coil systems generating fields with excellent quasisymmetry and low particle losses. In particular for low complexity coils, we are able to significantly improve the performance compared with coils obtained from the standard two-stage approach, e.g. reduce losses of fusion-produced alpha particles born at half-radius from $17.7\,\%$ to $6.6\,\%$ . We also demonstrate 16-coil configurations with alpha loss ${<}1\,\%$ and neoclassical transport magnitude $\epsilon _{\text {eff}}^{3/2}$ less than approximately $5\times 10^{-9}$ .

With the funding provided by this award, we developed numerical codes for the study of magnetically confined plasmas for fusion applications. Accordingly, our work can be divided into two separate categories: 1) the design and analysis of novel numerical methods providing high accuracy and high efficiency; 2) the study of the equilibrium and stability of magnetically confined plasmas with some of these numerical codes, as well as the study of the nature of the turbulent behavior which may arise in the presence of instabilities. We first developed new numerical schemes based on integral equation methods for the computation of steady-state magnetic configurations in fusion experiments, providing high accuracy for the magnetic field and its derivatives, which are required for stability and turbulence calculations. We employed different integral formulations depending on the application of interest: axisymmetric or non-axisymmetric equilibria, force-free or magnetohydrodynamic equilibria, fixed-boundary equilibria or free-boundary equilibria. While efficient, these methods do not yet apply to plasma boundaries which are not smooth, a situation which is fairly common in magnetic confinement experiments. To address this temporary weakness, we also constructed a new steady-state solver based on the Hybridizable Discontinuous Galerkin (HDG) method, which provides full geometric flexibility. In addition to these numerical tools focused on steady-states, we also contributed to the improvement of the speed and accuracy of codes simulating the plasma dynamics of fusion plasmas, by developing a novel velocity space representation for the efficient solution of kinetic equations, which most accurately describe the time evolution of hot plasmas in fusion experiments. Using the tools discussed above, we studied several questions pertaining to the equilibrium and stability of magnetically confined plasmas. In particular, we derived a new simple model for axisymmetric devices called tokamaks, to predict how elongated a fusion plasma can be before it becomes unstable and collapses. We also looked at the effect of the shape of the outer plasma surface on key properties of the steady-state magnetic configurations, and how these properties impact turbulence in fusion plasmas, and the corresponding transport of momentum. Likewise, we studied the role of large localized flows on the steady-state magnetic configurations, and how they may influence plasma stability and turbulence. Non-axisymmetric steady-state magnetic configurations are inherently more complex than axisymmetric steady-state configurations, and the subject of ongoing controversies regarding the regularity of the equations determining such steady-states, and their solutions. Implementing an existing NYU code in a new geometry, we studied the nature of the singularity of the solutions observed in the code, and methods to eliminate them. Our main conclusion is that by appropriately tailoring the plasma boundary, it is possible to eliminate the singularities otherwise appearing in our simulations, and to obtain steady-states which appear to be smooth. To gain further insights on incompletely understood turbulence phenomena, we proposed a new reduced model capturing most of these phenomena, which is simple enough to not require expensive numerical simulations on massive supercomputers to investigate them. We demonstrated the strong similarity between our simulations and published results obtained from computationally expensive simulations, and plan to rely on our reduced model to identify the key mechanisms determining the evolution and strength turbulent driven transport in fusion plasmas. Finally, we proposed a new framework for tokamak reactor design studies, enabling us to consider the relative merits of steady-state versus pulsed fusion reactors. We found that pulsed fusion reactors may benefit most from recent advances in magnet technology, and the availability of very high field magnets. As such, they may become more desirable than steady-state tokamak reactors for cost efficient electricity generation.

Alpine consists of a set of mini-apps that makes use of exascale computing capabilities to numerically solve some classical problems in plasma physics. It is based on IPPL (Independent Parallel Particle Layer), a framework that is designed around performance portable and dimension independent particles and fields. In this work, IPPL is used to implement a particle-in-cell scheme. The article describes in detail the following mini-apps: weak and strong Landau damping, bump-on-tail and two-stream instabilities, and the dynamics of an electron bunch in a charge-neutral Penning trap. We benchmark the simulations with varying parameters such as grid resolutions ($512^3$ to $2048^3$) and number of simulation particles ($10^9$ to $10^{11}$). We show strong and weak scaling and analyze the performance of different components on several pre-exascale architectures such as Piz-Daint, Cori, Summit and Perlmutter.