Nonlinear gyrokinetic predictions of SPARC burning plasma profiles enabled by surrogate modeling

As we are interested in the macroscopic behavior of the core plasma (temperature and density profiles to determine fusion performance), the transport equations to solve are written in the form of flux-surface average velocity moments of the electrons and ions distribution functions, and closed by providing a proper model for particle and energy fluxes. In this work, fluxes will be computed with high-fidelity collisional and turbulence models, as it will be described in detail in section 3.3. The separation between macroscopic (transport) and microscopic (turbulence) evolution is possible by virtue of δf gyrokinetic scale separation in devices with sufficiently small ρ_*, which is the case in many of current experiments and certainly in upcoming burning plasmas.

By integrating over the interior to the flux surface at radius r and assuming steady-state, the system of coupled equations to solve has the following form:

$$\begin{align}\hfill \langle {\mathbf{q}}_{\text{e}}\cdot \nabla r\rangle & =\frac{1}{{V}^{\prime }}{\int }_{0}^{r}\left(\langle {P}_{\text{aux},\text{e}}\rangle +\langle {P}_{\alpha ,\text{e}}\rangle -\langle {P}_{\text{ei}}\rangle -\langle {P}_{\text{rad}}\rangle \right){V}^{\prime }\,\mathrm{d}r\hfill \\ \hfill \langle {\mathbf{q}}_{\text{i}}\cdot \nabla r\rangle & =\frac{1}{{V}^{\prime }}{\int }_{0}^{r}\left(\langle {P}_{\text{aux},\text{i}}\rangle +\langle {P}_{\alpha ,\text{i}}\rangle +\langle {P}_{\text{ei}}\rangle \right){V}^{\prime }\,\mathrm{d}r\hfill \\ \hfill \langle {\mathbf{\Gamma }}_{\text{e}}\cdot \nabla r\rangle & =\frac{1}{{V}^{\prime }}{\int }_{0}^{r}\langle {S}_{\text{e}}\rangle {V}^{\prime }\,\mathrm{d}r,\hfill \end{align} \tag{ 1 }$$

where the r coordinate used here is the half-width of the mid-plane intercept, ⟨⋅⟩ indicates flux-surface averaging, and V' = ∂V/∂r where V is the volume interior to the flux surface. Sources and sinks are described in section 3.2, and the flux-surface averaged transport fluxes in section 3.3. The volume integrals in equation (1) result in the coupling of radial locations (i.e. a change in deposited power and particles in inner locations has a direct effect in outer locations), which is a consequence of the steady-state condition used in this work.

The elements in the left-hand side of equations (1) are descriptive of the transport fluxes, and will be referred to as ${Q}_{\text{e}}^{\text{tr}}$ and ${Q}_{\text{i}}^{\text{tr}}$ (in units of power per unit area), and ${{\Gamma}}_{\text{e}}^{\text{tr}}$ (in units of particles per second and unit area). The right-hand side includes target fluxes, and will be referred to as ${Q}_{\text{e}}^{\text{target}}$, ${Q}_{\text{i}}^{\text{target}}$ and ${{\Gamma}}_{\text{e}}^{\text{target}}$. A steady-state solution of the problem (also called flux-matched solution) is the set of electron temperature, ion temperature and electron density radial profiles that satisfy ${Q}_{\text{e}}^{\text{tr}}={Q}_{\text{e}}^{\text{target}}$, ${Q}_{\text{i}}^{\text{tr}}={Q}_{\text{i}}^{\text{target}}$ and ${{\Gamma}}_{\text{e}}^{\text{tr}}={{\Gamma}}_{\text{e}}^{\text{target}}$ at all radial locations. Therefore, a converged solution will satisfy equations (1) (i.e. will ensure energy and particle balance). The ansatz used in this work and the calculation of both transport and target quantities are based on the implementation used by the TGYRO code [20].

Auxiliary electron heating ⟨P_aux,e⟩ and ion heating ⟨P_aux,i⟩ are assumed unchanged during the transport evolution, together with the particle sources ⟨S_e⟩. This assumption is based upon the fact that we do not expect strong variations of the deposited power and particle profiles from electromagnetic waves and wall sources and recycling when the core temperature and density profiles vary. Nevertheless, the flexibility of the Bayesian framework presented here allows for the self-consistent variation of auxiliary power profiles, but it is left for future work.

The classical energy exchange between electron and ion species is defined as:

$$\begin{equation}\begin{aligned}\hfill \langle {P}_{\text{ei}}\rangle ={c}_{\text{ex}}\cdot \sum\limits _{\mathrm{i}}\frac{{Z}_{\text{i}}^{2}{n}_{\text{i}}{m}_{\text{i}}^{1/2}}{{\left({m}_{\text{e}}{T}_{\text{i}}+{m}_{\text{i}}{T}_{\text{e}}\right)}^{3/2}}\cdot \mathrm{ln}\,{\Lambda}\cdot {n}_{\text{e}}\cdot \left({T}_{\text{e}}-{T}_{\text{i}}\right)\end{aligned},\end{equation} \tag{ 2 }$$

where we use $\mathrm{ln}\,{\Lambda}=24.0-\mathrm{ln}\left(\sqrt{{n}_{\text{e}}}/{T}_{\text{e}}\right)$ and c_ex is a constant of proportionality. Detailed definitions and proper units can be found in reference [36]. The energy exchange power, as it enters in the electron and ion energy balance equations (equations (1)) will be varied self-consistently during the convergence process.

The total alpha heating, ⟨P_α ⟩ is calculated using the parameterization of the D–T fusion cross-sections proposed by Bosch and Hale [37]. The partition between electrons ⟨P_α,e⟩ and ions ⟨P_α,i⟩ is calculated as described by Stix [38]. In this model, the local power density that goes into the ions is given by:

$$\begin{equation}\begin{aligned}\hfill \langle {P}_{\alpha ,\text{i}}\rangle & =\langle {P}_{\alpha }\rangle .{\hat{W}}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}{\int }_{0}^{1/{\hat{W}}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}\frac{1}{1+{x}^{3/2}}\,\mathrm{d}x\hfill \\ \hfill {\hat{W}}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}& ={c}_{\alpha }\cdot {T}_{\text{e}}{\left(\frac{1}{{n}_{\text{e}}}\sum\limits _{\mathrm{i}}\frac{{Z}_{\text{i}}^{2}{n}_{\text{i}}}{{m}_{\text{i}}}\right)}^{2/3},\hfill \end{aligned}\end{equation} \tag{ 3 }$$

where c_α is a constant of proportionality [38].

Radiation ⟨P_rad⟩ as a sink for the electron energy channel will be calculated accounting for Bremsstrahlung, line and synchrotron radiations. The impurity assumptions used in this work are identical to those reported in reference [34], and each ion radiation contribution is computing using a Chebyshev polynomial representation of ADAS data [39, 40] for the average charge state of the impurity ion, as implemented in the TGYRO code [20]. Following this methodology, the cooling rate per impurity is given by:

$$\begin{equation}\mathrm{l}\mathrm{n}\;{L}_{Z}=\sum\limits _{k=0}^{n}{c}_{k}\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(k\cdot \mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{2\;\mathrm{l}\mathrm{n}({T}_{\text{e}}/{T}_{\mathrm{m}\mathrm{i}\mathrm{n}})}{\mathrm{l}\mathrm{n}({T}_{\mathrm{m}\mathrm{a}\mathrm{x}}/{T}_{\mathrm{m}\mathrm{i}\mathrm{n}})}-1\right)\right),\end{equation} \tag{ 4 }$$

where c_k are the Chebyshev coefficients (n = 11) and T_max and T_min define the fit range. The Bremsstrahlung and line radiated power density contribution per impurity is then ⟨P_rad,i⟩ = n_e n_i L_Z,i. Synchrotron radiation is calculated following reference [41].

Figure 1 depicts a comparison between the high-fidelity TRANSP [42] simulation of the SPARC PRD [34] and the targets model described in this section. Differences of less than 10% are found, despite the several orders of magnitude differences in run time. In particular, the predictions from the Stix model (equations (3)) are in remarkable agreement with the expensive, Monte-Carlo simulations with the NUBEAM code [43] in TRANSP.

Zoom In Zoom Out Reset image size

In summary, we have that each of the power densities can be written as nonlinear functions of the local kinetic profiles: ⟨P_ei⟩ = f(T_e, T_i, n_e), ⟨P_α,e/i⟩ = f(T_e, T_i, n_e) and ⟨P_rad⟩ = f(T_e, n_e). These local power densities enter in equation (1) in the form of volume integrals, which leads to radial and channel coupling (i.e. non block-diagonal Jacobian) when solving for the steady-state solution. Here, the normalized inverse gradient scale lengths are the free parameters to be used to find the flux-matching solution, and the local kinetic profiles are reconstructed using the trapezoidal rule, following the TGYRO methodology [20]. Likewise, volume integrals of the local power densities are performed using a realistic geometry parameterization [44].

The uniqueness of the solver presented here is the use of very high fidelity transport simulations. Instead of reduced quasi-linear transport models typical of flux-matching solvers, here we employ direct non-linear gyro-kinetic plus kinetic neoclassical particle and energy fluxes.

Turbulent fluxes were calculated using the continuum gyrokinetic code CGYRO [2]. CGYRO is a widely used code that has been optimized for modern computing architectures, and has been specifically designed for enabling multi-scale and near-edge and pedestal gyrokinetic simulations. Here, nonlinear simulations spanning much of the plasma minor radius (ρ_N = 0.25–0.85, where ρ_N is the normalized square root of toroidal flux coordinate) are performed. Unlike previous work that utilized nonlinear gyrokinetics for the prediction of plasma profiles (described in section 1), the simulations performed here included high physics fidelity, allowing for the simultaneous determination of particle and heat transport. Each simulation was performed utilizing six gyrokinetic species: deuterium, tritium, fluorine (Z = 9, A = 18, representing a lumped low-Z ion species), tungsten (partially ionized with Z = 50, A = 184), and a fast ³He population at 5% of the local electron density, representing the minority species for the ion-cyclotron heating. A population of alpha particles was notably not included in this simulations as it was determined that they represented a small fraction of the overall ion density and thus could be ignored without significant loss of fidelity.

The simulations performed in this work built off of previously reported nonlinear gyrokinetic simulations [45] which outlined the model fidelity and numerical setup needed for accurate simulation of the SPARC PRD. All simulations performed for this work were nonlinear and simulated electromagnetic turbulent fluctuations (evolving both δϕ and δA_|| contributions). The inclusion of electromagnetic effects was reported as potentially playing a significant role in the stabilization of ion temperature gradient (ITG) modes in the inner regions of plasma simulated here, and due to the minor increase in computational cost, they were included. Additionally, recent results [46] demonstrated a non-negligible impact of electromagnetic contributions to density peaking in predictive TGLF simulations.

The target simulation box sizes were L_x , L_y ∼ 80 × 80ρ_s. However, the exact box sizes that are attainable are determined by local equilibrium values. As a result, the minimum box size utilized was $\sim 80\times 80{\rho }_{\text{s}}$ with box sizes at inner radii generally being larger $(\sim 110{\rho }_{\text{s}})$ than this target. We note that the chosen box sizes where found to be sufficient for capturing the turbulence in the conditions studied, as discussed in reference [45]. Each simulation utilized approximately 300 radial modes and 16 toroidal modes, capturing ion-scale turbulence and zonal flows in the range k_θ ρ_s = 0.0–1.2 while also using N_ξ = 24 (pitch angle grid points), N_θ = 24 (poloidal grid points), and N_energy = 8 (energy grid points). Here we note that k_θ is the poloidal wavenumber of the turbulence and ρ_s is the ion Larmor radius evaluated with the sound speed, c_s. Nonlinear convergence tests were reported in a previous paper (reference [45]) on the initial profiles used in this work. Increases in box size and the number of radial modes yielded no statistically meaningful changes in the fluxes. Linear checks in other resolutions, (energy, pitch-angle, radial modes, and poloidal grid points) were performed on both the initial and converged profiles to provide confidence in the convergence of the results. These tests yielded small changes in the linear growth rates (most changes < 2%) with a 50% increase in each numerical resolution. The chosen resolutions for these simulations were an attempt to achieve accurate results while also minimizing the computing time utilized in this work. Plasma shaping properties were modeled using a MXH equilibrium [44] and collisions were modeled using the Sugama collision operator [47]. The effects of rotation were not included in these simulations but future work will include an estimated rotation profile to assess the impact of rotation and E × B shearing effects on the predicted profiles. However, it should be noted that such effects are generally expected to improve confinement and increase performance.

The simulation times required to obtain converged heat and particle fluxes at each of the radii simulated varied substantially, consistent with the large variation in the linear growth rates observed from ρ_N = 0.25 compared with 0.85. At the innermost radial position, nonlinear simulations were typically run for $\sim 1200a/{c}_{\text{s}}$ with averaging periods of approximately 400 − 500a/c_s, while the larger growth rates and therefore quicker saturation of the turbulence, allowed for approximately 600a/c_s simulations at ρ_N = 0.85 with average windows still targeting approximately 400a/c_s. Similar averaging windows were typically used throughout this work, however, there was variation due to the large variations in profiles and gradients that were studied and conditions that resulted in stable turbulence were sometimes terminated earlier.

The large variation in the profiles studied in this work resulted in dramatically different required time steps. This work took advantage of the adaptive time-stepping that was recently implemented in CGYRO to ensure that the appropriate numerical time step was utilized regardless of the profiles studied. Each simulation reported here utilized 72 nodes on the NERSC Cori KNL partition. Due to the variations in the length of the simulations, the exact computing time usage of these nonlinear simulations varied by location and profile studied. However, due to the large number of species simulated and resolutions utilized, on average we estimate that approximately 100 k CPU-hours were required per radial location, per profile studied. A number of methods (such as early stopping and reduction of number of ion species) can be utilized to reduce the overall computational time which are speculated to further reduce the computational resources required by approximately a factor of 2. These methods will be investigated as part of future work.

The use of ion-scale simulations in this work is motivated both by practical and physical considerations. Practically, the computational resources required to simultaneously resolve ion and electron-scales in a tokamak plasma would make profile prediction intractable, with the only existing simulations typically requiring 10–20 M CPU-hours per radial location [48–50]. Despite significant contributions to electron heating from fusion alphas, the conditions in the SPARC scenario studied here have been found to exhibit a large ratio of ion to electron heat flux Q_i/Q_e in the range 3–5 [34]. This balance of heat fluxes has been determined to occur in large part due to the strong collisional coupling of electrons and ions in the high density regime used in SPARC and due to the significant radiative losses that occur at the high operational temperatures. Large Q_i/Q_e ratios are a hallmark of dominant ion-scale turbulence, as electron-scale contributions (high-wavenumber trapped electron and electron temperature gradient modes) are responsible for near exclusively driving electron heat transport [51]. Furthermore, linear gyrokinetic and nonlinear electron-scale simulations were reported for the SPARC PRD profiles (that represent the initial profiles used in this work). These investigations suggested that the ratio of the electron to ion-scale growth rates (γ_high-k /γ_low-k ) [52] and the electron heat fluxes obtained by electron-scale simulations, generally indicate a negligible role of high-k turbulence in the SPARC PRD [45]. As a result, the use of ion-scale simulations for prediction of SPARC appears robust and should be expected to lead to an accurate description of profiles and performance.

Even though turbulence dominates the energy and particle losses in the cases presented here, the predictions also include second-order neoclassical transport fluxes as calculated with the multi-species kinetic NEO [53] solver.

Under the modeling framework of Gaussian process (GP) regression, energy and particle transport fluxes as a function of input parameters (i.e. kinetic profile gradients and relevant plasma quantities) are considered to follow independent joint Gaussian distribution functions. The GP that models such distribution is considered to have a linear mean, μ(x), and a covariance described by a Matern-5/2 kernel, $\mathcal{K}(\mathbf{x},{\mathbf{x}}^{\prime })$:

$$\begin{equation}{F}^{\text{tr}}(\mathbf{x})\sim \mathcal{G}\mathcal{P}(\mu (\mathbf{x}),\mathcal{K}(\mathbf{x},{\mathbf{x}}^{\prime })),\end{equation} \tag{ 5 }$$

where F^tr represents each of the local transport fluxes (i.e. electron energy flux, ${Q}_{\text{e}}^{\text{tr}}$, ion energy flux, ${Q}_{\text{i}}^{\text{tr}}$, and electron particle flux, ${{\Gamma}}_{\text{e}}^{\text{tr}}$), and x is the set of local gradients.

The solution of the transport flux evaluation with CGYRO is assumed to be polluted by heteroscedastic Gaussian noise, to account for the uncertainties in each of the evaluations of the true gyrokinetic transport fluxes (feature of initial value solvers such as the one used here):

$$\begin{equation}{y}^{\text{tr}}(\mathbf{x})={F}^{\text{tr}}(\mathbf{x})+\mathcal{N}(0,{\sigma }_{y}^{2}(\mathbf{x})).\end{equation} \tag{ 6 }$$

In this context, we assume that the standard deviation of the flux, σ_y (x), is known for each evaluation, and is calculated individually per profile and radial location simulated in CGYRO. Here we employ an approach that uses the auto-correlation time of the time series to determine the standard deviation of the fluxes, similar to previous work [54]. Each simulation time series was decomposed into windows that were three times wider than the calculated auto-correlation time and assumed to represent an independent sample. The standard deviation of the mean fluxes to be used in equation (6) was then reported as the standard deviation of independent samples divided by the square root of the number of samples.

This implementation of the likelihood function will be convenient when, in future work, a multi-fidelity approach to flux-matching is used. Similar to that proposed in reference [31], early stopping of the initial value simulations (before turbulence reaches full saturation) in those cases that are estimated to be far from targets could be used to save computational resources and can be assigned a larger noise to account for the evaluation uncertainty.

4.1.1. Physics-informed modeling

As described in section 3, three transport channels are solved for using n_r uniformly distributed radial locations in the ρ_N coordinate. This results in 3n_r design variables of the problem: electron temperature (∇T_e), ion temperature (∇T_i) and electron density (∇n_e) radial gradients at each of the simulated locations. This 3n_r–D optimization problem can, however, be simplified by taking into account that the flux evaluations in the context of this work are local, i.e. the flux at a given location can be described solely from the solution of transport equations at such location, and that long-range interactions only appear as a result of the imposition of the steady-state condition (i.e. in time-scales longer than turbulence decorrelation times). This means that the flux-matching technique takes care of radial coupling of the transport flux evaluations by ensuring a converged solution at the end, but each transport flux evaluation can be taken to be independent from each other at each iteration, and fully defined locally. By taking this into account, each local transport flux is assumed to be a function of six input quantities:

$$\begin{equation}{F}^{\text{tr}}(\mathbf{x})=f(\nabla {T}_{\text{e}},\nabla {T}_{\text{i}},\nabla {n}_{\text{e}},{T}_{\text{e}},{T}_{\text{i}},{n}_{\text{e}}),\end{equation} \tag{ 7 }$$

where the set of input parameters of the models are fully defined locally. At the end of this paper, a discussion on the locality assumption is provided.

To further simplify the problem and ease model training, we make the physics-informed assumption that turbulence within the nonlinear simulations is described in terms of the following normalized quantities:

$$\begin{equation}{F}^{\text{tr}}(\mathbf{x})={C}_{\text{GB}}\cdot \hat{{F}^{\text{tr}}}(a/{L}_{{T}_{\text{e}}},a/{L}_{{T}_{\text{i}}},a/{L}_{{n}_{\text{e}}},{T}_{\text{i}}/{T}_{\text{e}},{\hat{\nu }}_{\text{ei}}),\end{equation} \tag{ 8 }$$

where $a/{L}_{{T}_{\text{e}}}$, $a/{L}_{{T}_{\text{i}}}$ and $a/{L}_{{n}_{\text{e}}}$ are the normalized inverse gradient scale lengths. C_GB represents the gyro-Bohm normalization, either energy flux ${Q}_{\text{GB}}\propto {n}_{\text{e}}{T}_{\text{e}}^{5/2}$ or particle flux ${{\Gamma}}_{\text{GB}}\propto {n}_{\text{e}}{T}_{\text{e}}^{3/2}$. And ${\hat{\nu }}_{\text{ei}}$ is the electron–ion collision frequency normalized to c_s/a, so that ${\hat{\nu }}_{\text{ei}}\propto \mathrm{ln}\,{\Lambda}\cdot {n}_{\text{e}}{T}_{\text{e}}^{-2}$. The functions that map normalized inputs to normalized outputs, $\hat{{F}^{\text{tr}}}$, will be the ones to be modeled with the GPs.

Figure 2 shows an example of fitting independent GP models with Matern kernel and linear mean to transport fluxes at several core radial locations in 1D scans of ion temperature gradient from nonlinear CGYRO simulations. GPs provide a good way to capture different behaviors (including critical gradient) without the need to provide functional forms and associated hyper-parameters. The non-parametric nature and the ability to include uncertainties in the observed data points make GP models suitable for the problem at hand. Even though the burning plasmas of interest in this work are often dominated by ITG-type turbulence [34, 45, 56], the non-parametric surrogates employed here do not preclude other types of regimes, or even transitions in between dominant turbulence modes (likely at the expense of more iterations in such complex situations), as, in principle, the GP models could include any arbitrary non-linear behavior with respect to the input parameters. While this work utilizes standard kernels and mean functions used widely in the BO literature, it is possible that functions specific to fusion applications would lead to even better performance, following the experience of past work with GP regression in the field [57], but this and the investigation of whether GPs are suitable to reproduce transport in any arbitrary turbulence regime is left for future work.

Zoom In Zoom Out Reset image size

Equation (8) indicates that the surrogate modeling problem translates into 3n_r independent 5D GP models. As it will be shown next, this choice works well for the context of this problem, but it could be generalized by the addition of further normalized quantities, at the expense of having to perform more observations to properly fit the GP model ⁴ . Figure 3 depicts the schematic of the surrogate model that results from implementing the physics-informed GP model as described earlier. For GP training and evaluation, the GPyTorch package [58] is used. The entire model described by equation (8) and every operation depicted in figure 3 are implemented with PyTorch [59] tensors, which allows to leverage automatic differentiation [60] for both training of the GP models and the evaluation of the full Jacobian of the system for optimization and root finding.

Zoom In Zoom Out Reset image size

Note that the physics-informed surrogates as presented in equation (8) are only a function of the number of channels to simulate, n_p, and not a function of the number of radial locations, n_r. With a finer radial grid, the model would require more input gradients to find the solution, but the accuracy of the local surrogates would not be compromised with the same number of complete profiles simulated (i.e. the surrogate models remain 5D). This suggests that the framework presented here has excellent scalability and parallelization properties that will be exploited in future work.

Despite leveraging the use of normalized gyro-Bohm units to fit the turbulent transport surrogates, the minimization of the transport residual is done in real units. This approach is advantageous when extending the predictive core region to the far edge. The gyro-Bohm energy flux is very sensitive to the temperature, ${Q}_{\text{GB}}\propto {T}_{\text{e}}^{5/2}$, which results in extreme differences between edge and core (in particular when solving for L-mode-type plasma conditions). This is typically very challenging for flux-matching solvers, as the edge transport residual dominates the flux-matching process, often requiring extremely small step sizes in, e.g., Newton solvers. Using real units to minimize the residuals thus provides a clear advantage, as the fluxes to match are typically of the same order of magnitude between radial locations.

To find the predicted plasma profiles in steady-state, the transport fluxes must match the target fluxes at each radial location. This condition can be written as:

$$\begin{equation}{\mathbf{x}}^{\ast }=\underset{\mathbf{x}{\epsilon}[{\mathbf{x}}_{\text{L}},{\mathbf{x}}_{\text{U}}]}{\mathrm{arg min}}\vert {\mathbf{F}}^{\text{tr}}(\mathbf{x})-{\mathbf{F}}^{\text{target}}(\mathbf{x})\vert ,\end{equation} \tag{ 9 }$$

where F^tr and F^target represent the set of electron energy, ion energy and electron particle transport and target fluxes at all radial locations, x is the set of all gradients and [x_L, x_U] are the allowed bounds. This results in a multi-objective minimization problem of the 3n_r transport residuals. We remind the reader that the target fluxes F^target are not fixed. The expected exchange power, radiation, and alpha heating are affected by the changes in kinetic profiles and are updated at each iteration. Sometimes this is referred to as a 'moving targets' calibration problem, and contributes to the high cost of finding flux-matching conditions with standard optimization methods and finite-differences-based Jacobian calculations.

The framework built here to solve for the flux-matching solution is based on standard BO techniques. A small set of kinetic profiles (typically 5 to 10) is used to train the initial set of GP models to reproduce each of the transport fluxes. This training set is produced by means of Latin hypercube sampling to ensure efficient filling of the parameter space even when a low number of samples are used. Even though 5 to 10 initial training samples may seem very low to fit 5D surrogate models, they are typically sufficient to drive the system near marginal stability in stiff plasma conditions. Further refinement to the surrogate models occurs later on to understand where the critical gradient sits at exactly. At present, the method requires that the final converged gradients are within the original training bounds, so that the optimization method only searches in regions of the multi-dimensional parameter space where the model has been trained. Future work will focus on devising methods to find solutions outside of the original bounds.

With the GP models in hand, the flux-matching solution of the GP representation can be obtained by leveraging standard optimization techniques. This work performed optimization based on the mean of the posterior distribution. A genetic algorithm (GA) based on the DEAP library [61] is employed to generate a set of initial solutions by minimizing the mean of the posterior distribution of the average residual of all channels and all radial locations. A multi-start local minimizer based on the L-BFGS-B algorithm [62] as implemented in the BoTorch toolset [63] is used to refine the GA solutions. Finally, a multi-variate root-finder based on the modified Powell method, as implemented in scipy [64] is used to further refine the solutions. The use of the exact Jacobian without the need of finite difference approximations greatly reduces the computational cost of the optimization process, and eliminates the need to specify step sizes for finite differences and prevents oscillatory behavior.

Once a new flux-matching solution is found for the surrogate GP models, the best set of kinetic profiles is evaluated with the high-fidelity transport simulations. If convergence (i.e. actual flux-matching) has not been attained, the profile is added to the database of training samples for the next iteration. Because of this, the surrogate models become more accurate around the region of the parameter space that is expected to contain the final converged solution (around critical gradients). Further details on similar surrogate-based optimization methods for physics and engineering applications can be found in references [65–68].

The critical gradient and stiff transport nature of plasmas often complicates convergence of flux-matching solvers. A very small variation of gradients can lead to large variations in predicted fluxes, and thus reducing the residual to very small levels often requires unnecessary long simulations if the purpose is predicting performance. As shown in figure 4, fusion gain (strong function of density and temperature) oscillated only around $\sim 5\%$ of the final value for the last five iterations, even though the residual changed by over 100%.

For this reason, as a stopping criterion for our simulations, we not only accounted for a low enough average residual, but we also developed a metric based on the stiffness of transport:

$$\begin{equation}{\Delta}a/{L}_{X}^{\ast }\equiv \vert {\mathbf{F}}^{\text{tr}}({\mathbf{x}}_{\text{i}})-{\mathbf{F}}^{\text{target}}({\mathbf{x}}_{\text{i}})\vert \cdot {\left[\frac{\partial {\mathbf{F}}^{\text{tr}}}{\partial a/{L}_{X}}\right]}_{({\mathbf{x}}_{\text{i}})}^{-1},\end{equation} \tag{ 10 }$$

evaluated at each iteration optimum point, x_i, after model training. Here, a/L_X refers to each of the driving gradients ($a/{L}_{{T}_{\text{e}}}$, $a/{L}_{{T}_{\text{i}}}$ and $a/{L}_{{n}_{\text{e}}}$) for each channel.

In short, equation (10) refers to the average relative change in normalized gradients that is predicted to provide flux matching by a linearized diagonal version of the transport model. This metric is valid only if the change in gradient is small enough, which is the case toward the end of the flux-matching process. By the end, as shown in figure 4(l), the average of this metric for all channels and radii is 3.4% and it was reached at profile #17. An extra profile (#18) was evaluated to demonstrate that the fusion gain remains unchanged. Unfortunately, this extra evaluation did not reduce the residual nor the convergence metric. For pragmatic reasons, we decided to stop the flux-matching process, as we did not consider that further reduction of the residual would lead to any meaningful change to the profiles, as demonstrated in figure 5, and even less to the fusion performance prediction.

Zoom In Zoom Out Reset image size

Although in our current research, we strive to have near-exact flux matches, if one wants to use this technique to project fusion power, confinement time and plasma performance in conditions of stiff transport such as the one considered here, the convergence criteria can be relaxed. Figure 5 clearly highlights the stiff conditions in SPARC PRD plasma. Minimal changes in the gradients for the last few iterations of the present framework led to significant changes in the fluxes that oscillated around the targets, even reaching stable conditions at some radii (e.g. profile #15). The profiles, however, barely changed and therefore the fusion gain remained the same.

Figure 6 compares the prediction from CGYRO to that of different TGLF saturation rules. Although the prediction of ion temperature with SAT2 is in remarkably good agreement with CGYRO, the electron energy transport is overpredicted for the inner half of the plasma minor radius. All rest of saturation rules predicted lower levels of ion energy transport in the near edge. However, the strongest and most important difference between CGYRO and TGLF predictions is in the electron density profile. All saturation rules resulted in predictions of density profiles that are less peaked in the plasma core than with CGYRO, consequence of the overall lower gradients at around mid-radius. As a result of the very different density peakings, the fusion gain widely varies from one prediction to another, despite having similar levels of the global energy confinement factor, H_98,y2. However, it is encouraging that the first-principles prediction with CGYRO yields a fusion gain of Q ≈ 8, greatly surpassing SPARC's core mission (Q > 2), even ignoring the favorable effect of intrinsic toroidal rotation in SPARC [45, 72], possible stabilization effects of fast alpha populations [73], and stabilization of ITG potentially provided by simulation of δB_|| perturbations [74]. The evaluation of the effect of rotation, fast ions, and the stabilization of ITG due to compressional magnetic perturbations (δB_||) in the prediction of SPARC scenarios is left for future work.

Zoom In Zoom Out Reset image size

Even though POPCON, TGLF and CGYRO profiles all feature similar quality of energy confinement (H_98,y2 = 0.95–1.05), the predictions of fusion gain can differ by a factor of 2. We expect that such differences will be even larger when simulating reactor plasmas that operate somewhat near ignition conditions. This is due to two important elements: (1) density peaking, and (2) temperature gradient scale length profile shape. This work suggests that predicting or extrapolating energy confinement (or pedestal height) alone is not enough for reliable predictions of fusion power, and calls for the need to use first-principles core transport simulations or to develop reduced models around the region of the parameter space of interest for a reactor study, as well as extensive experimental validation studies of profile predictions in present devices.

The new prediction of density peaking factor (evaluated at ρ_ψ = 0.2, where ρ_ψ is the square root of the normalized poloidal flux, following the nomenclature of references [75, 77]) with CGYRO is closer to the experimental observations in present devices than any of the TGLF predictions, as depicted and described in figure 7. A database of JET and ASDEX Upgrade experimental profiles was used to derive an empirical scaling for the density peaking based on the effective collisionality, normalized pressure and a core particle source parameter [75]. Inclusion of Alcator C-Mod data, which had no neutral beam injections, was demonstrated to also be captured reasonably well by the same empirical scaling [76]. Evaluation of SPARC parameters yields ν_ne ≡ n_e(ρ_ψ = 0.2)/⟨n_e⟩ ≈ 1.48. CGYRO predictions produced a density profile with ν_ne ≈ 1.45, which is in remarkable agreement with the empirical scaling. Although not conclusive, the improved agreement with empirical density peaking results provides additional confidence that the SPARC PRD profiles predicted in this work are representative of those expected during SPARC operation.

Zoom In Zoom Out Reset image size