Enhancing predictive capabilities in fusion burning plasmas through surrogate-based optimization in core transport solvers

Abstract

This work presents the PORTALS framework, which leverages surrogate modeling and optimization techniques to enable the prediction of core plasma profiles and performance with nonlinear gyrokinetic simulations at significantly reduced cost, with no loss of accuracy. The efficiency of PORTALS is benchmarked against standard methods, and its full potential is demonstrated on a unique, simultaneous 5-channel (electron temperature, ion temperature, electron density, impurity density and angular rotation) prediction of steady-state profiles in a DIII-D ITER Similar Shape plasma with GPU-accelerated, nonlinear CGYRO. This paper also provides general guidelines for accurate performance predictions in burning plasmas and the impact of transport modeling in fusion pilot plants studies.

Figures (16)

• Figure 1: Visualization of problem geometry and free parameters of a three-channel prediction. a) Illustration of nested flux surfaces for a DIII-D ITER Similar Shape plasma howard_simultaneous_2024, with magnetic axis in dashed black line, 5 selected flux surfaces in blue (cross section in purple) and last closed flux surface in green. Note that the lower x-point region is smoothed out due to Fourier-moments decomposition typical of transport solvers such as Breslau2018. Subplots b), d) and f) depict the piece-wise linear function representations of the normalized logarithmic gradient of electron temperature, ion temperature and electron density with the value at 5 selected flux surfaces (values to be predicted by the transport solver). Subplots c), e) and g) depict the corresponding kinetic profiles of each channel, produced by radial integration of logarithmic gradient from an edge boundary condition (dark blue).
• Figure 2: Performance of at achieving $10\times$ and $100\times$ residual reduction with in two example H-mode plasmas: DIII-D ITER Similar Shape (ISS) and ITER Baseline howard_simultaneous_2024. Random and SR initialization methods are compared, starting from the same initial condition. (top) Residual vs evaluations, with vertical line delimiting initialization phase and purple horizontal lines indicating the $10\times$ and $100\times$ residual reductions. (bottom) Violin plots representing distribution of number of evaluations that were required to achieve residual goals. Both methods were initialized with a random seed, and distributions were obtained for 16 seeds. As expected, SR is not strongly affected by the seed, particularly at the beginning of the flux matching process.
• Figure 3: Summary of effect of radial grid on integrated exchange power calculation on an example L-mode plasma. Electron temperature, ion temperature and electron density are parameterized with piecewise linear logarithmic gradients. The effect of three grids on the calculation of targets are shown ([blue] full profile, [green] 20 grid points and [red] 5 grid points). While kinetic profiles and local power density overlap at grid points, as expected, the integrated exchange power can have significant deviations due to the volume integration.
• Figure 4: Comparison of residual evolution in with random initial training, a standard Newton method (NM) and simple relaxation (SR) schemes. For this exercise, we employed the quasilinear Trapped Gyro-Landau Fluid () SAT0 model Staebler2007. Three example plasmas were used: (a) DIII-D ISS howard_simultaneous_2024, (b) ITER Baseline howard_simultaneous_2024 and (c) JET Ohmic L-mode prf_2023_eps. was initialized with different random seeds, which affect the initial training (5 profiles, vertical dashed line) and the subsequent training of Gaussian processes and acquisition optimization techniques. NM corresponds to method-$1$ and each of the three numerical parameters (relaxation parameter $\eta$, maximum step size $b_{max}$ and finite differences step $\Delta_x$) were varied within reasonable ranges. SR corresponds to method-$6$ with variations of the relaxation parameter $\eta_{SR}$. We note that, even though works with gyro-Bohm normalized residuals, here we plot them in real units, $MW/m^2$, for a direct comparison of performance.
• Figure 5: Characteristic phases during the prediction of three burning plasmas ([left] SPARC Primary Reference Discharge (PRD) Rodriguez-Fernandez2020a, [center] ARC holland_2023_aps and [right] ITER howard_2024_iter) with nonlinear . Phase I is the initial surrogate training phase, Phase II represents the first iterations (often with oscillatory behavior as it learns key parametric dependencies) and Phase III is the final convergence phase. The transition from II to III is smooth and not characterized by a change in the iteration scheme but it represents a transition to a situation where the surrogates are capturing the turbulence dynamics with enough accuracy to drive the system towards steady state. Bottom subplots show the evolution of the residual per radial location separately, shaded from red (core) to blue (edge).