arXiv CSopen access

Cycle-Consistent and Uncertainty-Aware Neural Surrogates for Tokamak Edge Plasmas

Cycle-Consistent and Uncertainty-Aware Neural Surrogates for Tokamak Edge Plasmas A. Diawa,∗, S. De Pascualea , J.-S. Parka , I. Paradela Pereza , J.D. Lorea , S. Dasbachb a Fusion Energy Division, Oak Ridge National Laboratory, Oak Ridge, TN 37932, USA b DIFFER - Dutch Institute for Fundamental Energy Research, De Zaale 20, 5612 AJ Eindhoven, the Netherlands

Abstract

arXiv:2607.21407v1 [physics.comp-ph] 23 Jul 2026

The boundary and divertor plasma play a key role in how a tokamak removes power and particles. They set the heat fluxes, temperatures, densities, and the start of detachment. Predicting these values accurately is crucial for safely running current and future devices. However, detailed edge simulations that resolve these parameters are too slow for tasks such as parameter scans, optimization, or real-time control. To address this, machine learning surrogates are now often used instead of traditional edge-plasma simulations. Still, most standard models can only predict forward and cannot recover input parameters from observed data or check how reliable their predictions are. In this study, we introduce a cycle-consistent neural surrogate for edge plasmas. It combines a conditional U-Net forward model with an optimization-based inverse method that uses the frozen forward network. The forward model takes five control parameters and predicts two-dimensional plasma-state fields on the SOLPS-ITER mesh. The inverse method ensures consistency between forward and inverse predictions, offering a self-supervised quality check that does not need groundtruth labels during inference. We also train a group of multilayer perceptrons to predict electron temperature and density profiles at the outboard midplane and divertor targets, with uncertainty estimates. The variation among the committee members provides a reliability measure for real-time control and helps spot areas where more simulations are needed. The forward model achieves normalized root-mean-square errors below 2.6% and Pearson correlations above 0.95 for all plasma-state fields. Adding cycle-consistency regularization raises the average cyclical R2 from 0.59 to 0.99 without reducing forward accuracy, and it allows recovery of the core fueling rate Γcore , which is hard to determine with forward-only training. The inverse method recovers all five control parameters with Pearson r ≥ 0.97. Using a k-d tree warm start to build the database yields a completion rate above 95%, whereas in a comparable cold-started ensemble roughly 30% of runs failed outright and a further fraction never reached steady state. With about 4 × 106 parameters, the model can generate full two-dimensional predictions in milliseconds, which is five to six orders of magnitude faster than the original SOLPS-ITER runs. This speed is enough for real-time control, thorough parameter scans, uncertainty analysis, and digital-twin applications. Keywords: SOLPS-ITER, edge plasma, neural surrogate, cycle consistency, uncertainty quantification, inverse modeling

1. Introduction Surrogate models have long been a workhorse of computational physics, standing in for first-principles solvers when direct numerical simulation is too slow to explore or optimize a problem at scale [1, 2, 3, 4]. They are most useful precisely where simulation is hardest: multiscale physics, stiff governing equations, and the need for large parameter scans or real-time inference. But this is also where a surrogate is most likely to fail. A surrogate is only trustworthy if it reproduces the system where its behavior is sharpest, near thresholds, bifurcations, and steep ∗ Corresponding author

Email address: [email protected] (A. Diaw) 1 This manuscript has been authored in part by UT-Battelle, LLC,

under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The publisher acknowledges the US government license to provide public access under the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

spatial gradients; these are exactly the regions that are hardest to learn. Scrape-off-layer (SOL) plasma is a textbook case. In the SOL, plasma–neutral and plasma–material interactions are set by kinetic processes, detailed atomic and molecular rates, and sheath physics acting on short spatial and temporal scales. Yet the quantities that matter for divertor design, scenario optimization, and feedback control, namely the target heat fluxes, the density profiles, and the radiation patterns, live at the macroscopic scale of the whole boundary and evolve over milliseconds to seconds. Predictively bridging these scales is essential for current and future devices, yet it remains computationally challenging [5, 6]. For SOL physics, the state-of-the-art tool is SOLPSITER, which couples a two-dimensional Braginskii fluid solver (B2.5) to a kinetic Monte Carlo neutral transport model (EIRENE) on a triangulated vessel mesh [7]. This

multiphysics code has been used extensively to design and interpret experiments on existing devices and to evaluate divertor concepts for future machines. It is also computationally expensive. Each run solves a stiff nonlinear fixedpoint problem whose wall-clock time often exceeds several hours, whose convergence is highly sensitive to initialization, and whose failure rate climbs sharply when one scans broad parameter sets or enters detached regimes. For these reasons, exhaustive sweeps over input power, gas puff, and transport coefficients are computationally prohibitive, and the solver cannot be placed directly inside an optimization loop or a real-time workflow. Data-driven surrogates for edge plasmas have received growing attention, using reduced models or neural networks trained on high-fidelity simulations to interpolate across parameter space [8]. These surrogates already recover the essential features, including detachment, target fluxes, and profile shapes. Two limitations persist, however. Most are confined to one-dimensional profiles at a handful of fixed locations, and they rarely enforce any physical consistency between the forward and inverse mappings. As a result, a prediction arrives with no built-in way to tell whether it can be trusted, which is exactly what one needs in the sharp, undersampled regimes where the surrogate is most likely to be wrong. We present here a cycle-consistent neural surrogate for two-dimensional edge-plasma fields. At its core is a conditional U-Net forward model Fθ that maps five scalar control parameters to the two-dimensional plasma state (Te , Ti , ne , ua ). A gradient-based inverse procedure then recovers the control parameters from a target field by optimizing through the frozen forward model. A cycle-consistency constraint ties the two together: parameters recovered from a field, when pushed back through Fθ , must reproduce that field. This round-trip supplies a self-supervised quality metric that requires no ground-truth labels at inference time. We choose a deterministic, parameter-conditioned UNet. Edge-plasma databases are data-scarce, with a few hundred converged runs at best. In this regime, it is prudent to learn a single-valued operator from control parameters to fields rather than a full distribution. Unlike generative approaches such as variational autoencoders or autoregressive transformers, the U-Net carries no sampling variance and does not have the large data appetite of those models [9, 10]. Our choice is pragmatic, and we believe it is well motivated. The resulting model has roughly 4.3 million parameters and returns a full two-dimensional prediction within milliseconds, fast enough for rapid scenario evaluation and for deployment in real-time plasma control [11, 12] and digital-twin workflows [13, 14]. Alongside this two-dimensional model, we train a queryby-committee ensemble of multilayer perceptrons for the one-dimensional Te and ne profiles at the outboard midplane and the divertor targets. We restrict it to Te and ne as an initial demonstration of the method, and because these are the quantities measured directly by Thom-

son scattering, facilitating future validation. The ensemble supplies predictive uncertainty estimates for closedloop control, and its committee variance flags the undersampled regimes where new SOLPS-ITER runs are most needed, a query-by-committee strategy we have used in earlier work [15, 16]. Together, the two surrogates cover the two deployment modes we care about: design and inverse analysis on one side, control and active-learning acquisition on the other. The structure of this paper is as follows. Section 2 details the simulation database and the k-d tree warmstart strategy. Section 3 describes the surrogate architecture, including the conditional U-Net forward model, the inverse inference procedure, and the cycle-consistency framework. Section 4 evaluates forward model accuracy, inverse parameter inference, and cycle-consistency metrics on held-out test cases. Section 5 presents the uncertaintyquantified one-dimensional profile ensemble for control and active learning. Section 6 interprets the results, and Section 7 summarizes the conclusions and outlines future directions. 2. Dataset The training dataset was generated using SOLPS-ITER [17, 18] simulations of a deuterium plasma in DIII-D, employing the same lower-single-null configuration as used in Ref. [19]. This setup was developed using DIII-D discharge 174310 at t = 3500 ms, has been demonstrated to have excellent numerical stability, and has been used to assess divertor conditions and neutral dynamics [19, 20, 21]. Consequently, it serves as a convenient baseline for constructing surrogate models of the coupled plasma–neutral system and for comparison with other SOL machine-learning studies [22, 23, 24, 25]. To span a representative operational space, we sample the input parameters using Latin hypercube sampling (LHS) over the bounds in Table 1. Each sample corresponds to a five-dimensional input vector: c = (Ptot , Γcore , ΓD2 , D⊥ , χi ), which sets the total power across the separatrix, the target core flux, the deuterium puff rate, and the uniform crossfield particle and heat transport coefficients. While LHS provides good space-filling coverage at coarse resolution, each simulation run effectively solves a stiff nonlinear fixed-point problem whose convergence is highly sensitive to the initial condition. In practice, the chosen restart file strongly influences the solver’s convergence trajectory: depending on the location in parameter space, runs may require many iterations, stall, or fail to converge numerically. Nontrivial failure rates have been reported in previous work. In the large cold-started (fluid-neutral) ensemble of Dasbach and Wiesen [22], about 29% of runs (1198 of 4096) diverged outright, and a further fraction of 2

the surviving runs never reached steady state within the allotted runtime. To mitigate these issues, we adopt a branching k-d tree strategy [26] that both partitions the input space and warm-starts each simulation from the nearest alreadyconverged design point, propagating reliable initial conditions outward from a single seed run. This raises the completion rate above 95%, compared with the roughly 70% completion of the cold-started ensemble above. The construction and distance metric are detailed in Appendix Appendix A. The detailed convergence detection procedure, including autocorrelation analysis and block averaging, is described in Appendix Appendix B. The completed database comprises 762 converged runs, of which two were discarded during quality screening, leaving 760 for surrogate modeling. We split these 80/20 into a training/validation set of 608 runs and a held-out test set of 152 runs. The split is performed at the level of whole simulations, so all mesh cells of a given run fall entirely within either the training or the held-out set, preventing information leakage. All the results shown in this paper are on the test set.

quality metric that does not require ground-truth parameter labels at test time. Each stage is described in detail in the following subsections. 3.1. Conditional U-Net Forward Model The forward surrogate Fθ maps scalar control parameters c and the spatial binary mask m to the multi-channel plasma state

The surrogate framework follows the manifold and cycleconsistency paradigm. At inference, it runs as a threestage cycle (Fig. 1): a forward model, an optimizationbased inverse, and a cycle-consistency self-check. All three stages share a single trained forward model Fθ and a learned pseudo-inverse Gψ , which are produced beforehand by a training-time cycle-consistency regularizer (Sec. 3.4, Fig. 2). The workflow has three components: (1) a conditional UNet Fθ is trained on the SOLPS-ITER database to map scalar control parameters c = (Ptot , Γcore , ΓD2 , D⊥ , χi ) and a spatial binary mask m to the four-channel twodimensional plasma state. Once trained, Fθ defines the learned physics manifold: every output ŷ = Fθ (c, m) is, by construction, a physically plausible edge-plasma state, (2) an inverse model: given a target field y∗ , an optimizationbased inverse procedure G⋆ recovers the control parameters ĉ by minimizing the discrepancy ∥Fθ (ĉ, m) − y∗ ∥ through the frozen forward model via gradient descent, and (3) a cycle-consistency stage in which the recovered parameters are passed back through
Fθ to obtain a reconstructed field ŷcycle = Fθ G⋆ (y∗ ), m . The cycle-consistency loss ∥ŷcycle − y∗ ∥ measures whether the round-trip y∗ → ĉ → ŷcycle preserves the target, providing a self-supervised

Table 1: Control parameters and their sampling ranges for the SOLPS-ITER database.

min max

Γcore (at/s)

Γ D2 (at/s)

D⊥ (m2 /s)

χi (m2 /s)

2 16

1 × 1020 7.5 × 1020

1.3 × 1020 5 × 1021

0.1 2.0

0.1 2.0

=

Fθ (c, m),

(1)

ŷ

=

{Te , Ti , ne , ua },

(2)

on the SOLPS (R, Z) grid, where Te and Ti are the electron and ion temperatures, ne is the electron density, and ua is the parallel velocity of the deuterium ions. Although the SOLPS state arrays are dense in index space (i, j), the physical plasma domain is irregular in (R, Z). It occupies only a subset of the rectangular tensor used by convolutional neural networks (CNNs) [27, 28]. Standard UNet convolutions operate on all grid locations, so without masking, cells outside the valid plasma region would still contribute to normalization and loss, potentially biasing training toward non-physical/background regions. The binary mask, therefore, acts as a domain-of-validity selector: optimization and evaluation are restricted to physically meaningful mesh cells while preserving compatibility with efficient regular-grid convolutions. The scalar parameters c are standardized using the mean and standard deviation of the training set and injected into the network via Feature-wise linear modulation (FiLM) [29]: at each encoder and decoder level, the parameter vector is mapped through a learned affine layer that produces the per-channel scale and shift coefficients applied after group normalization. This yields a singlechannel spatial input (the binary mask alone) while enabling the network to learn a parametric nonlinear operator across multiple plasma regimes. Each output channel is standardized independently, with statistics computed on the training split over valid mesh cells only. Because the temperatures and density are positive and span several orders of magnitude, they are logtransformed before being centered and scaled to unit variance (with a small offset, 10−2 for Te , Ti and 1016 m−3 for ne , to regularize near-zero values). The parallel velocity ua changes sign along the flow and is instead standardized with a symmetric-log transform before centering and scaling. The five input scalars are standardized with a plain mean/standard-deviation z-score. Encoder. The network follows a U-Net architecture [10] with three resolution levels. At each level, an input feature map is processed by two 3 × 3 convolutions, each followed by group normalization and a SiLU activation:
Y′ = σ GN(W ∗ Y) , (3)

3. Surrogate Model Design

Ptot (MW)

ŷ

σ(s) 3

=

s · sigmoid(s).

(4)

Figure 1: Surrogate model architecture: the three-stage inference cycle. (a) Stage 1 (forward model): a conditional U-Net maps a singlechannel geometry mask, with the five scalar control parameters injected via FiLM, to the four-channel plasma state (Te , Ti , ne , ua ); the fields 4 are outputs of the model, not inputs. (b) Stage 2 (inverse model): the control parameters are recovered from a target field by optimizing through the frozen forward model, warm-started by the learned pseudo-inverse Gψ (Fig. 2) evaluated on the observed fields. (c) Stage 3 (cycle consistency): the recovered parameters are passed back through the frozen model; the cycle loss validates that the round-trip preserves the target field, a self-supervised metric requiring no ground-truth parameters. The forward model Fθ and pseudo-inverse Gψ used here are produced beforehand by the training-time regularizer of Fig. 2.

Downsampling is performed using 2 × 2 max pooling with 3.3. Cycle Consistency a stride of 2; at each step, the spatial size halves while the The cycle-consistency constraint links the forward and number of channels doubles. The bottleneck applies the inverse models by requiring that parameters recovered from same double-convolution block at 8× base channels. a target field, when fed back through the forward model, Decoder. Each up-step uses a 2 × 2 transposed convoreproduce that target field. Formally, for a ground-truth lution to increase spatial resolution and reduce channels, field y∗ with mask m: concatenates the corresponding encoder feature map via
a skip connection, and then applies two 3 × 3 convolu(6) Lcycle = Fθ G⋆ (y∗ ), m − y∗ , tions with group normalization and SiLU. After the last where G⋆ (y∗ ) denotes the inverse-recovered parameters decoder block, a 1 × 1 convolution maps the base channels (Eq. 5) and ∥ · ∥ is a masked norm over valid mesh cells. to Cout = 4 output channels corresponding to the four This constraint serves several purposes. First, it reguplasma state fields. Because the targets are continuous, larizes the inverse mapping: even when the forward model’s we use a linear activation at the final layer. Additionally, parameter-to-field mapping is locally flat (so that many because ne ≡ ni in the plasma, we omit the ion density parameter vectors produce similar fields), the cycle loss channel. For thermodynamic consistency, we enforce non-negativity penalizes inverse solutions whose forward reconstructions deviate from the target. Second, it enforces manifold conof the temperature and density outputs (Te , Ti , ne ≥ 0) by sistency [32, 33]: the reconstructed fields ŷcycle = Fθ (ĉ, m) clamping negative values during post-processing. are guaranteed to lie on the learned physics manifold of the forward model, ensuring physically plausible outputs. 3.2. Inverse Model Third, the cycle loss provides a self-supervised quality metThe inverse problem, in which we infer control paramric for the inverse procedure that does not require ground∗ eters c from target plasma fields y , is ill-posed: many truth parameter labels at test time. parameter combinations can yield similar field configurations. Rather than training a separate inverse network, we adopt an optimization-based approach that leverages the differentiability of the frozen forward model Fθ . We denote this optimization-based inverse G⋆ ; unlike Fθ and the pseudo-inverse Gψ (Sec. 3.4), it has no trainable parameters of its own: it recovers ĉ by optimizing the control parameters directly through the frozen forward model. Given a target field y∗ and spatial mask m, we solve for the parameters ĉ that minimize the discrepancy between the forward prediction and the target: ĉ

=

Linv (c′ )

=

3.4. Cycle Consistency as a Training Regularizer The optimization-based inverse of Section 3.2 acts only at inference time and therefore leaves the forward model unchanged. To test whether cyclical consistency can additionally improve the forward surrogate itself, as reported for inertial-confinement-fusion surrogates [33], we introduce a learned pseudo-inverse Gψ : y 7→ c and train it jointly with the forward model under

min ρ Fθ (c, m), y + ρ Gψ (y, m), c θ,ψ | {z } | {z } inverse forward
(7) + λcyc c − Gψ Fθ (c, m), m , | {z }

arg min Linv (c′ ), c′ X
2 ∗ wc mij Fθ (c′ , m)c,i,j − yc,i,j c,i,j

cycle

X

wc mij

c,i,j

+λreg ∥c′ ∥2 ,

(5)

where wc are the per-channel weights and λreg is an L2 regularization coefficient that prevents the recovered parameters from drifting to unphysical extremes. Optimization is performed using Adam [30] with learning rate 10−2 for 1200 steps through the frozen forward model, with gradient-based parameter updates computed via automatic differentiation [31]. The optimization is warm-started from the pseudo-inverse estimate ĉ0 = Gψ (y∗ , m) (Sec. 3.4), which is computed from the observed fields alone. To mitigate local minima, we employ Nr = 5 random restarts about this estimate, perturbed with Gaussian noise (σnoise = 0.2 in standardized space), and retain the solution with the lowest Linv ; no ground-truth parameters are used at any point. 5

where ρ is a masked smooth-L1 discrepancy. Unlike Eq. 6, the cycle term here back-propagates into the forward parameters θ, so consistency can reshape the surrogate. The weight λcyc multiplies only the cycle term; the forward and inverse reconstruction terms are kept at unit weight, so λcyc = 0 recovers the forward-only baseline. The pseudo-inverse Gψ is a compact convolutional encoder that mirrors the forward encoder in reverse. Its input is the four-channel plasma state stacked with the geometry mask as a fifth channel. Three convolutional stages increase the channel width from 24 to 48 to 96; each stage applies two 3 × 3 convolutions, each followed by group normalization and a SiLU activation, and then a 2 × 2 max-pooling that halves the spatial resolution. A global average pooling collapses the final feature map to a 96-dimensional vector, which a two-layer perceptron (hidden width 128) maps to the five control parameters. The

Figure 2: Training-time cycle-consistency regularizer. The forward model Fθ and the learned pseudo-inverse Gψ are trained jointly; the cycle penalty λcyc ∥c − Gψ (Fθ (c, m), m)∥ back-propagates into both networks, so consistency reshapes the forward surrogate. Setting λcyc = 0 recovers the forward-only baseline; the weight is swept in the ablation of Sec. 4.3.

network exists only to seed the inference-time optimization (Sec. 3.2) and to supply the cycle-consistency signal during training; it is not used as a standalone inverse. We quantify self-consistency with an average cyclical R2 score analogous to the metric of Anirudh et al. [33]: for each control parameter, we sweep it linearly across its range (holding the others fixed), push the resulting parameter vectors through Fθ and then Gψ , and measure the coefficient of determination between the swept and recovered values; the score is averaged over the five parameters. A value near unity indicates that the forward and inverse mappings are mutually consistent on held-out single-parameter scans.

and apply a linear warmup from epoch 20 to 60 to stabilize early training. Training uses the Adam optimizer [30] with initial learning rate 3 × 10−4 , ReduceLROnPlateau scheduling (factor 0.5, patience 5 epochs), mixed-precision (AMP) on GPU, gradient clipping at norm 1.0, and early stopping with patience of 80 epochs. The model is trained for up to 450 epochs with a batch size of 8, and the best checkpoint (by validation loss) is retained. All input and output fields are standardized using training-set statistics. 4. Results Before turning to the results, we fix the metrics used throughout. We quantify the accuracy of a predicted field with four complementary measures, all computed over the valid mesh cells of the held-out test set. Writing yi for the SOLPS-ITER value and ŷi for the prediction at cell i (over P N valid cells), the mean absolute error is MAE = N −1 i |ŷi − yi |, and the normalized root-mean-square error is NRMSE = 100 × RMSE/(ymax − ymin ), the RMS error written as a percentage of the field’s range. We complement these magnitude errors with two association measures: the Pearson correlation r (linear agreement) and the Spearman correlation ρ. For the ablation of Sec. 4.3 and the one-dimensional ensemble of Sec. 5 we also report the coefficient of determination X X R2 = 1 − (ŷi − yi )2 / (yi − ȳ)2 , (9)

3.5. Loss Functions and Training The forward model Fθ is trained on the SOLPS-ITER database using a composite loss function that balances pixel-level accuracy, gradient fidelity, and boundary emphasis: L = Lbase + λw Ledge + λg Lgrad . (8) The base loss Lbase is a masked Huber (smooth-L1 ) loss with β = 0.05, computed only over valid mesh cells indicated by the binary mask; per-channel weights wc allow prioritization of specific output fields, and we use wc = 1.0 for Te , Ti , wc = 1.2 for ne , and wc = 1.5 for ua . The edge loss Ledge is the same Huber loss reweighted by a boundary proximity map w(r), computed as a Gaussian-decayed distance from the mask boundary (σ = 3 pixels); this emphasizes accuracy near the separatrix and target plates where gradients are steepest. Finally, Lgrad is a Sobelfilter-based spatial gradient penalty that encourages the model to reproduce sharp spatial features (recycling fronts, temperature pedestals). The gradient is computed on both predicted and target fields, and the loss is the masked L1 difference of the resulting gradient maps. We set λg = 0.2

i

i

where ȳ is the mean of the SOLPS-ITER values. Separately, we use robustness to refer to the reliability of the pipeline rather than its pointwise error. At the database level it is the completion rate, the fraction of SOLPS-ITER runs that reach steady state (Sec. 2); at the model level it is the sensitivity of the prediction to input perturbations, E ∥Fθ (c) − Fθ (c + σϵ)∥2 in normalized output space 6

(Table 3), for which a smaller value indicates a smoother, more robust surrogate.

learned pseudo-inverse from scratch. Figure 4(A) and Table 3 report the held-out test mean squared error and the average cyclical R2 score as λcyc is increased, following the ablation protocol of Anirudh et al. [33]. As λcyc increases, the forward and inverse networks become markedly more cyclically self-consistent: the average cyclical R2 increases from 0.59 at λcyc = 0 to 0.99 at λcyc = 0.5. Crucially, this improved self-consistency does not come at the expense of forward accuracy: the held-out test MSE remains comparable to or below the forward-only baseline (λcyc = 0) across the entire sweep. The prediction’s sensitivity to input perturbations (Table 3) remains bounded and comparable to the unregularized baseline across the range explored. The effect is clearest at the level of individual control parameters. Figure 4(B) compares the per-parameter recovery R2 of the forward-only baseline with that of the consistency-regularized model (λcyc = 0.5). Without the cycle term, the trained pair is self-consistent for some control parameters but not others: the recovered R2 for the core particle source Γcore and the cross-field diffusivity D⊥ is low, whereas the input power Ptot and the gas-puff rate ΓD2 are already well recovered. Training with the cycle term makes the pair self-consistent across all five parameters, each above 0.97. We emphasize that this improved recovery is a property of the jointly trained forward/inverse pair rather than of the forward model in isolation: it does not by itself imply that the forward model has become more invertible. We use it only as a self-supervised consistency signal and adopt λcyc = 0.5, which is both the most self-consistent and the most accurate on held-out data.

4.1. Forward-model accuracy Figure 3 compares SOLPS-ITER ground truth with UNet predictions for the four plasma-state fields (Te , Ti , ne , ua ) in a held-out test case. The model accurately reproduces the two-dimensional structure of electron and ion temperatures, including the steep gradients near the separatrix and the characteristic decay into the scrapeoff layer. Electron density is well captured across the full domain, with the core–SOL contrast and divertor compression faithfully represented. The parallel velocity field ua shows correct flow patterns from the outer midplane toward the divertor targets. Quantitatively, the forward model achieves a global Pearson correlation exceeding 0.95 for all four plasma fields: r = 0.994 for Te , r = 0.992 for Ti , r = 0.987 for ne , and r = 0.975 for ua (Table 2). Mean absolute errors are 12.2 eV for Te and 16.7 eV for Ti , which represent small fractions of the dynamic range of these fields across the database. The density field achieves an MAE of 2.6 × 1018 m−3 for ne , while the parallel velocity MAE is 941 m/s. Residual errors are predominantly localized near the separatrix and in the private flux region where spatial gradients are steepest. Table 2 summarizes the forward-model performance across all four output channels on the held-out test set. The four plasma state fields (Te , Ti , ne , ua ) all achieve Pearson correlations above 0.95 and NRMSE below 2.6%, with per-sample mean Pearson values exceeding 0.98. The parallel velocity ua is the hardest of the four to reproduce (r = 0.975, NRMSE 2.54%): its sign changes along the flow, and it carries sharper spatial structure than the temperatures and density, so residual errors concentrate near the separatrix and the divertor targets.

4.4. Inverse parameter recovery We evaluate the inverse model on held-out test cases using the optimization-based inference procedure (Sec. 3.2). To initialize the optimization we use the learned pseudoinverse Gψ of Sec. 3.4: its convolutional encoder maps the observed fields, with the geometry mask supplied as an auxiliary channel, to a parameter estimate ĉ0 = Gψ (y∗ , m) that provides the warm start; no ground-truth parameters are used. Starting from this estimate, the optimizer refines the parameters by minimizing the forward-model discrepancy (Eq. 5) through the frozen Fθ , with Nr = 5 random restarts per case and cosine learning-rate annealing over 1200 Adam steps. The cycle-consistency metric Lcycle (Eq. 6) is then evaluated by passing the recovered parameters through Fθ and comparing the reconstructed fields against the original targets. This procedure provides two complementary diagnostics: (i) the accuracy of parameter recovery (how close ĉ is to the true c∗ ), and (ii) the cycle reconstruction quality (how well Fθ (ĉ, m) matches y∗ ). These need not agree: the inverse problem in edge-plasma modeling can be ill-posed, since the forward mapping is locally insensitive to certain parameter combinations, so good cycle reconstruction does not on its own guarantee accurate recovery. Reporting both diagnostics lets us separate the two and determine

4.2. Model size and inference speed The conditional U-Net has ∼4.3×106 trainable parameters with a base filter width of 48, occupying 16.5 MB in single precision. Each forward pass maps the singlechannel mask input, with the five scalar parameters injected via FiLM, to four output channels at 36 × 96 resolution. On a single CPU core, inference takes 16 ms per sample; on a GPU, sub-millisecond throughput is achievable in batched mode. With GPU inference, these latencies are compatible with real-time plasma control loops (1–10 ms cycle times) [11, 12, 34], digital-twin frameworks [13, 14], and Monte Carlo uncertainty quantification over the input space. 4.3. Invertibility gains from cycle consistency To test whether cyclical consistency improves the forward surrogate, and not merely the inverse, we sweep the regularization weight λcyc in the joint objective (Eq. 7) and, for each value, retrain the forward model with the 7

Figure 3: Comparison of SOLPS-ITER ground truth (left columns), U-Net predictions (middle columns), and absolute errors (right columns) for the four plasma-state fields on a held-out test case.Top to bottom: electron temperature Te (eV), ion temperature Ti (eV), electron density ne (m−3 ), and parallel velocity ua (m/s). The model captures the large-scale structure and gradients across the SOL and divertor.

8

Table 2: Forward surrogate performance on the held-out test set (152 runs). MAE is reported in physical units, and NRMSE (%) is the root-mean-square error normalized by the range of each field (ymax − ymin ) over the valid test cells. Pearson r and Spearman ρ are global correlations computed over all valid points across all test samples. All metrics are defined at the start of Sec. 4.

Field

Unit

MAE

NRMSE (%)

Te Ti ne ua

eV eV m−3 m·s−1

12.24 16.71 2.61 × 1018 940.9

1.26 1.54 0.62 2.54

Pearson r

Spearman ρ

0.994 0.992 0.987 0.975

0.996 0.987 0.964 0.976

Figure 4: Cycle-consistency ablation. (A) Average cyclical R2 (left axis, higher is better) and held-out test MSE (right axis, lower is better) as the cycle weight λcyc is increased. Self-consistency increases from 0.59 at the forward-only baseline (λcyc = 0) to 0.99. At the same time, the test MSE stays at or below the baseline (dotted line) throughout, so self-consistency improves at no cost to forward accuracy. (B) Perparameter self-consistency of the trained forward/inverse pair for the forward-only baseline (λcyc = 0, orange) versus the cycle-consistent model (λcyc = 0.5, blue). The cycle term brings every control parameter above 0.97, including Γcore and D⊥ , which are otherwise weakly recovered by the pair. Higher R2 is better; lower MSE is better. Each λcyc is trained from scratch, so the intermediate points carry run-to-run variation; the baseline (λcyc = 0) and the adopted λcyc = 0.5 are the comparison of interest. Exact sweep values are given in Table 3.

whether a clean field reconstruction reflects a genuinely well-identified parameter set. Figure 5 shows the true versus recovered values for the five control parameters on held-out test cases. All five are recovered well, with Pearson correlations at or above 0.97: Ptot (r = 0.99, ρ = 0.99), ΓD2 (r = 0.99, ρ = 0.99), D⊥ (r = 0.99, ρ = 0.97), χi (r = 0.99, ρ = 0.97), and the core fueling Γcore (r = 0.97, ρ = 0.95). The core fueling remains

the most scattered of the five, consistent with the weaker sensitivity of the downstream SOL and divertor fields to it. Still, it is well identified: the inverse optimization works through the frozen forward model. Hence, a parameter to which that model is only weakly sensitive is the hardest to pin down. Reporting both parameter-recovery accuracy and cycle-reconstruction quality is what distinguishes a genuinely well-identified parameter from a degenerate one. 5. Uncertainty-Quantified Profile Ensemble for Control and Active Learning

Table 3: Cycle-consistency ablation sweep (760-run database; 608 training / 152 test runs, trained from scratch at each λcyc ). As the cycle weight increases, the average cyclical R2 increases toward unity while the held-out test MSE stays comparable to or below the forward-only baseline (λcyc = 0). Sensitivity is the mean squared change of the normalized prediction under Gaussian input perturbations (σ = 0.1); a smaller value indicates a more robust surrogate. Best values in bold. λcyc Test MSE Cyclical R2 Sens. (×10−3 )

0 (baseline) 0.01 0.05 0.1 0.5

0.0288 0.0272 0.0274 0.0287 0.0266

0.585 0.760 0.895 0.975 0.989

Alongside the two-dimensional U-Net, the framework’s second surrogate is a standalone query-by-committee ensemble of multi-task regression networks that predicts onedimensional profiles of electron temperature Te and density ne at the outboard midplane and both divertor targets directly from the five control parameters, as a function of the normalized poloidal flux ψN . We restrict the targets to Te and ne deliberately: these are the independent state variables measured directly by Thomson scattering and are the quantities most relevant to edge diagnostics and real-time control. Derived quantities such as the divertor heat flux are not given a separate learned head, since

2.14 3.93 3.97 5.10 3.51 9

Figure 5: Inverse parameter recovery on held-out test cases. Each panel shows the true versus recovered value for one control parameter, obtained by optimizing through the frozen forward model with multiple random restarts. Pearson and Spearman correlation coefficients are annotated per parameter.

with 10% held out for early-stopping validation and a further held-out fold used to score it. Members that do not reach R2 ≥ 0.90 on every output channel of that fold are discarded. We retain nensemble = 5 trained models. Because ne spans several decades and the target Te ranges from below 1 to above 200 eV, both are trained on their log10 values and mapped back to physical units for evaluation. To assess the trustworthiness of the model prediction, we threshold the committee disagreement with a quality score σi si = k , (10) S where σi is the ensemble standard deviation at point i and k = 2. The global scale S = σ̄ Ēcal /Ēens rescales the raw committee spread to the magnitude of the actual prediction error: σ̄ is the mean ensemble standard deviation, Ēcal the mean absolute error of the ensemble-mean prediction, and Ēens the mean ensemble spread, all computed on the held-out validation fold. If si ≥ 1, the point is flagged for further simulation. We train each ensemble on 12160 points (608 SOLPSITER runs) and evaluate it on 3040 held-out test points (152 runs), using the same 80/20 train/test split as the 2D surrogate so that the two models are directly comparable. The performance on the held-out test set is shown in Fig. 6. The ensemble reproduces the held-out Te and ne profiles with high fidelity, with the per-panel coefficient of determination annotated in Fig. 6. Points are colored by the quality flag si ; those selected for SOLPS-ITER verifi-

they are constitutive functions of Te and ne and predicting them independently would risk thermodynamic inconsistency; when needed, they can instead be evaluated from the predicted profiles with uncertainty propagated through the ensemble. Where the U-Net provides physically complete fields for design and inverse analysis, this ensemble targets a different deployment mode: fast, pointwise profile evaluation with uncertainty estimates, as required for real-time control and for active-learning selection of the most informative new SOLPS-ITER runs [22, 24]. SOLPS-ITER simulations are large (many cells, many time steps), so we favor parametric models whose evaluation cost does not grow with the size of the training set. We also prefer methods that admit a simple, effective UQ scheme. We therefore use an ensemble of neural networks and take the ensemble variance as a proxy for model uncertainty. This is a practical example of the query-bycommittee (QBC) approach [35], an active-learning strategy that proposes new training points where a committee of models disagrees the most. We have used this approach in previous work [15, 16]. Training minimizes mean-squared error loss, X (t) (t)
2 L = ŷi − yi i,t

using Adam [30] in mini-batches of size 256, with early stopping (patience 20 epochs, learning rate halving) and a cap of 400 epochs. We form the ensemble by bootstrapping: each member is trained on a different random split, 10

cation (si > 1) are drawn as open circles, highlighting the regions where the network is least confident. Among the three locations, the inner divertor target is reproduced least accurately, with the lowest Te coefficient of determination (R2 ≈ 0.89, versus ≳ 0.94 at the upstream and outer-target locations) and by far the largest fraction of flagged points (∼ 13%, against ∼ 3% at the outer target). This is consistent with edge physics rather than a shortcoming of the network: the inner leg detaches at lower upstream density than the outer, so across the database it spans both attached and detached states, and near the strike point the target Te becomes effectively bimodal and steeply varying, a harder mapping to learn from the control parameters alone. The committee makes this visible rather than hiding it: the inner target carries the largest ensemble disagreement. It is therefore the top priority for active-learning acquisition of new SOLPS-ITER runs.

through the differentiable inverse procedure, replacing iterative manual SOLPS-ITER fitting [36]; and (iii) realtime or near-real-time edge-plasma state estimation for control applications [11, 12], conditioned on a limited set of diagnostic measurements. Comparable digital-twin paradigms have been demonstrated for particle accelerator control [37, 34]. Generating the SOLPS-ITER database is computationally demanding: each run solves a stiff, coupled nonlinear system whose convergence trajectory is highly sensitive to the initial state. The k-d tree restart strategy (Sec. 2) functions as a discrete form of numerical continuation [38], where the solution at one parameter-space point predicts the solution at a nearby point. Classical continuation traces a single curve through parameter space; here, the k-d tree generalizes this concept to a multidimensional Latin hypercube design [39] by selecting, for each new sample, the nearest previously converged neighbor in (log N ) time [26, 40]. The database grows as a branching process: starting from a single seed, each newly 6. Discussion converged case extends the frontier of reliable initial condiOur framework builds on the manifold-and-cycle-consistencytions, progressively reducing the parameter-space distance to subsequent targets. Beyond this branching construcparadigm for surrogate modeling of tokamak edge plastion, the warm-start strategy improves robustness: in a mas. Anirudh et al. [33] enforce manifold consistency with prior cold-started ensemble, roughly 30% of runs initiated a Wasserstein autoencoder that learns a low-dimensional from generic initial conditions diverged outright, and a latent representation; we instead use a conditional U-Net further fraction never reached steady state [22], whereas whose skip connections preserve fine spatial detail. Both the branching approach achieved a completion rate above share the core insight that cycle consistency regularizes 95% across the full five-dimensional design, yielding the the inverse model, and in both, the reconstructed outputs 762-run database used here. This infrastructure is essenare constrained to lie on the learned physics manifold. The tial for any surrogate relying on a large, uniformly samtokamak setting adds its own difficulties: the mesh geomepled training set; without it, gaps in the database would try is irregular, which is what forces the masking described directly translate into blind spots in the learned model. above, and the fields carry hundreds of sharp recycling Several limitations remain, and we state them plainly. fronts that do not admit a smooth, low-dimensional strucFirst, the current database covers only a single DIII-D ture. lower-single-null configuration with deuterium-only fuelTable 2 shows that all four fields are reproduced acing and spatially uniform transport coefficients; extendcurately (Pearson correlation > 0.95), with a mild hiering it to other machines, geometries, plasma mixtures, archy among them: the temperatures and density, which and spatially varying transport is necessary for validation vary smoothly over the scale of the scrape-off layer (SOL) and, perhaps, for transfer learning. Second, the surrogate width, are captured most efficiently by the convolutional has not yet been validated against experimental DIII-D receptive field, whereas the parallel velocity ua carries sharper measurements, a critical step before deployment. Third, spatial structure and a sign change along the flow and is physics constraints beyond admissibility, such as global ratherefore the hardest to fit. That same spatial structure is diated power and Z consistency, and the conservation of eff what motivates the convolutional inductive bias: coherent particles, momentum, and energy, are not yet enforced; large-scale flow patterns like ua are reproduced far more imposing them is a natural way to sharpen the physical accurately by a model that shares information between fidelity of the surrogate. neighboring cells than by a pointwise regressor. This supFinally, although the current inverse procedure is warmports the choice of a deterministic, parameter-conditioned started by the learned pseudo-inverse Gψ (Sec. 3.4), the reconvolutional operator for the data-scarce regime studied finement still requires gradient-based optimization through here, where generative alternatives would have to learn a the frozen forward model; end-to-end joint training with full distribution from only a few hundred samples. cycle consistency could further improve convergence speed The millisecond-scale inference described in Sec. 4.2 enand inverse accuracy. ables several practical applications: (i) rapid exploration of the five-dimensional parameter space for sensitivity analysis and scenario optimization [15, 3]; (ii) transport-coefficient and parameter inference from diagnostic measurements

11

Outer target

18.9 18.6

1.0

R²=0.968

20.4

NN log10 ne [m−3]

NN log10 ne [m−3]

NN log10 ne [m−3]

19.2

19.8 19.2 18.6

20.8

20.0

0.8

19.2

18.4

18.3 18.3

18.6

18.9

19.2

19.5

18.0 18.0

SOLPS log10 ne [m−3]

18.6

19.2

19.8

20.4

300 150

102

2

NN Te [eV]

NN Te [eV]

450

101

0

150

300

450

SOLPS Te [eV]

600

20.0

20.8

R²=0.894

0.4

101

0.2

100

100

0

19.2

SOLPS log10 ne [m−3]

R²=0.942 10

0.6 18.4

SOLPS log10 ne [m−3]

600 R²=0.937

NN Te [eV]

Inner target

R²=0.955

Quality Flag si

Upstream 19.5 R²=0.963

100

101

102

SOLPS Te [eV]

100

101

102

SOLPS Te [eV]

Figure 6: Evaluation of the one-dimensional profile ensemble. Columns correspond to the three profile locations (upstream/outboard midplane, outer divertor target, and inner divertor target); rows correspond to the two predicted quantities, electron density ne and electron temperature Te , on the independent test set (3040 points per location; 12160 training points). Target Te is shown on a logarithmic axis; the per-panel coefficient of determination is annotated in each panel. All points are colored by the quality flag si ; open circles mark the highest-uncertainty points (si > 1) and are flagged for additional SOLPS-ITER simulation.

12

netic neutral step, is developed separately.

7. Conclusions

A cycle-consistent neural surrogate framework is introCRediT authorship contribution statement duced for modeling two-dimensional SOLPS-ITER edgeplasma fields. The framework incorporates a k-d tree warmA. Diaw: Conceptualization, Methodology, Software, start strategy that selects nearest-neighbor restarts from Formal analysis, Investigation, Data curation, Visualizaa database of previously converged SOLPS-ITER simulation, Writing – original draft. S. De Pascuale: Contions, thereby reducing both simulation failure rates and ceptualization, Methodology. J.-S. Park: Conceptualwall-clock time: across the five-dimensional parameter scan, ization, Methodology. I. Paradela Perez: Conceptualthe warm start reached a completion rate above 95%, yieldization, Methodology. J.D. Lore: Software, Resources, ing 762 converged runs, whereas roughly 30% of comparaMethodology, Funding acquisition. S. Dasbach: Concepble cold-started runs have been reported to fail outright, tualization, Methodology. with more never reaching steady state [22]. A parameter-conditioned U-Net maps scalar control parameters directly to the two-dimensional plasma state Code and Data availability (Te , Ti , ne , ua ), and achieves Pearson correlation coeffiThe code is available at https://github.com/abdoudiaw/ cients exceeding 0.95 for all four plasma state variables. solpex. The sampled data is archived on figshare at https: Cycle-consistency regularization is shown to yield a self//doi.org/10.6084/m9.figshare.32048490. consistent forward/inverse pair at negligible cost to forward accuracy, supplying a self-supervised consistency signal that requires no ground-truth labels. Finally, an inAcknowledgments verse inference procedure coupled with cycle-consistency regularization enables parameter recovery from target fields. This work was supported by the U.S. Department of It ensures that reconstructed outputs lie on the learned Energy (DOE), Office of Science, Office of Fusion Enphysics manifold, providing a self-supervised reliability met- ergy Sciences. This research used resources of the Oak ric that requires no ground-truth labels at inference time. Ridge Leadership Computing Facility at Oak Ridge NaFor held-out cases, the inverse recovers all five control pational Laboratory, which is supported by DOE under Conrameters with Pearson correlations of 0.97 or higher, intract DE-AC05-00OR22725. This research used resources cluding the cross-field transport coefficients D⊥ and χi of the Oak Ridge National Laboratory Research Cloud. (both r ≈ 0.99); since these are typically set by hand-tuned The results are obtained with the help of the EIRENE SOLPS-ITER fits to diagnostics, the differentiable inverse package (see www.eirene.de) including the related code, offers an automated route toward transport-coefficient indata and tools [18]. DIFFER is part of the institutes orference. A companion query-by-committee profile ensemganisation of NWO. ble supplies predictive uncertainty estimates whose committee disagreement concentrates on the physically hardDeclaration of generative AI and AI-assisted techest regime, the near-strike-point inner divertor closest to nologies in the writing process detachment (flagging roughly 13% of inner-target points against 3% at the outer target), and thereby nominates During the preparation of this work, the authors used the most informative new SOLPS-ITER runs for activeClaude (Anthropic) to improve the manuscript’s readabillearning acquisition. ity and language. After using this tool, the authors reLooking ahead, several directions are being pursued. viewed and edited the content as needed and took full reFirst, validation against experimental DIII-D measurements sponsibility for the published article. from discharge 174310 is underway, which will test the surrogate’s ability to generalize beyond the simulation database. Appendix A. k-d Tree Warm-Start Second, extending the inverse procedure to operate on sparse one-dimensional diagnostic signals (e.g., Thomson Construction. We build a balanced k-d tree over scattering profiles) would directly replace the manual transportall LHS design points {xi } by recursively splitting the coefficient fitting workflow of interpretive boundary moddata along coordinate directions. This structure supports eling [36]. Third, incorporating additional physics connearest-neighbor queries in O(log N ) time on average for straints (global radiated power Prad and effective charge d ≪ N , where d = 5 is the number of input parameters Zeff ) and extending to spatially varying transport coeffiand N is the number of design points. Given a target input cients will improve the surrogate’s physical fidelity and apt⋆ , we define a weighted Euclidean distance plicability. Finally, extending the approach to other magnetic configurations and machines (e.g., ITER, SPARC, 5 X 2 MAST-U) will test its generality. The companion neutrald2w (xi , t⋆ ) = wj (xi,j − t⋆,j ) , (A.1) source model, which maps the plasma state predicted here j=1 to EIRENE source terms for inline replacement of the ki13

as a noisy relaxation process whose fluctuations encode the residual coupling between the plasma and neutral subsystems. −1 Given a time series {qk }N k=0 with sample mean q̄, we compute the normalized autocorrelation function NX −ℓ−1

C(ℓ) =

(qk − q̄) (qk+ℓ − q̄)

k=0 N −1 X

, (qk − q̄)

ℓ = 0, 1, 2, . . . ,

2

k=0

(B.1) so that C(0) = 1 and |C(ℓ)| ≤ 1 for ℓ > 0. The integral correlation time is then estimated as ! ℓmax 1 X + C(ℓ) , (B.2) τint ≈ ∆t 2 ℓ=1

where the sum is truncated at the first non-positive value of C(ℓ) or at a prescribed maximum lag [46]. This quantity plays the same role as the correlation time in MD block averaging [45]: it sets the minimum window length over which successive block means can be treated as approximately independent. Using τint , we define a tail region consisting of the last Ttail = Ntail τint timesteps of the series and split it into two non-overlapping windows of equal length L = Nwin τint , denoted A and B. The drift between the two window means, ∆q̄ = q̄B − q̄A , is compared against its expected statistical fluctuation, r √ 2 τint σ∆q̄ ≈ 2 σq , (B.3) L where σq is the standard deviation within each window. A quantity is declared to be in a steady state when

Figure A.7: The k-d tree warm-start strategy applied during database generation, projected onto the (Ptot , ΓD2 ) plane. Grey points are the existing converged database; colored markers are new samples, each initialized from its k-d tree nearest neighbor among the already-converged runs. Arrows link each new sample to the neighbor it was warm-started from, and marker color encodes the batch index, so the branching growth of the database frontier is visible as sampling progresses.

where wj is a user-defined weight for the j-th parameter (here wj = 1 for all j). For each target, we query the tree for the nearest already-converged neighbor and initialize the SOLPS-ITER state from that solution rather than from a generic seed (Fig. A.7). The database thus grows as a branching process. Each newly converged case extends the frontier of reliable initial conditions, a discrete form of numerical continuation [38] generalized to a multidimensional Latin hypercube design [39].

|∆q̄| < εrel |q̄B |

Appendix B. Database Convergence Detection

and

|∆q̄| < zmax , σ∆q̄

(B.4)

with εrel = 0.01 and zmax = 2. The first condition ensures that the mean has not drifted by more than 1%; the second ensures that the observed drift is statistically consistent with equilibrium fluctuations at the estimated correlation time. A simulation is accepted as converged only when all four diagnostics simultaneously satisfy Eq. (B.4). Figure B.8 illustrates the procedure for a representative run, showing the time series, the two comparison windows, and the corresponding autocorrelation functions. This approach gives a consistent steady-state check inline while SOLPS-ITER is running.

Determining whether a SOLPS-ITER run has reached steady state is nontrivial. Density, momentum, energy, and radiation fields relax toward equilibrium at different rates, and the iterative coupling between B2.5 and EIRENE introduces correlated fluctuations that persist over many timesteps. We address this using autocorrelation and blockaveraging techniques from equilibrium sampling in molecular dynamics [41, 42]. The same machinery underpins Green-Kubo transport-coefficient calculations in molecular dynamics simulations [43, 44]; here, we repurpose it to estimate the correlation time of SOLPS-ITER diagnostics and construct a statistically grounded convergence criterion [45, 46]. We monitor four scalar diagnostics extracted at each SOLPS-ITER timestep: the electron and ion temperatures OMP OMP at the outboard-midplane separatrix (Te,sep , Ti,sep ), the OMP electron density at the same location (ne,sep ), and the total number of particles (⟨tmne⟩). Each diagnostic is treated

Appendix C. One-Dimensional Profile Cycle Consistency To complement the two-dimensional cycle-consistency results of Sec. 4.4, we illustrate the same forward/inverse round-trip at the level of the one-dimensional profiles extracted from the conditional U-Net. For a representative 14

Figure B.8: Steady-state diagnostics for a representative SOLPS-ITER run. Top row: time series of the four monitored quantities with windows A and B used for the drift test. Bottom row: normalized autocorrelation functions with the estimated integral correlation time τint (dashed line).

held-out test case (Fig. C.9), the U-Net reproduces the SOLPS-ITER profiles with reasonable fidelity, recovering the correct amplitudes and characteristic profile shapes, both when evaluated at the true control parameters (forward) and when using the parameters recovered by the inverse procedure (inverse cycle). We emphasize that this forward/inverse comparison is a property of the U-Net surrogate of Sec. 3; it is distinct from the query-by-committee profile ensemble of Sec. 5, which performs no parameter inversion.

[5] R. Schneider, X. Bonnin, K. Borrass, D. P. Coster, H. Kastelewicz, D. Reiter, V. A. Rozhansky, B. J. Braams, Plasma edge physics with B2-Eirene, Contributions to Plasma Physics 46 (1-2) (2006) 3–191. doi:10.1002/ctpp.200610001. [6] R. A. Pitts, X. Bonnin, F. Escourbiac, T. Hirai, V. Komarov, A. S. Kukushkin, A. Loarte, A. Martin, M. Merola, R. Mitteau, Divertor and scrapeoff layer plasma physics for reactor-grade fusion energy devices, Nuclear Materials and Energy 20 (2019) 100696. doi:10.1016/j.nme.2019.100696.

References

[7] D. Reiter, M. Baelmans, P. Börner, EIRENE – a Monte Carlo linear transport solver, Tech. rep., Forschungszentrum Jülich, http://www.eirene.de (2005).

[1] P. Wigley, P. Everitt, A. van den Hengel, J. Bastian, M. Sooriyabandara, G. McDonald, K. Hardman, C. Quinlivan, P. Manju, C. Kuhn, I. Peterson, A. Luiten, J. Hope, N. Robins, M. Hush, Fast machine-learning online optimization of ultracold-atom experiments, Scientific Reports 6 (25890) (2016).

[8] J. Kates-Harbeck, A. Svyatkovskiy, W. Tang, Machine learning for disruption warnings on Alcator CMod, DIII-D, and EAST, Nature 568 (7753) (2019) 526–531. doi:10.1038/s41586-019-1116-4.

[2] A. Scheinker, S. Gessner, Adaptive method for electron bunch profile prediction, Physical Review Accelerators and Beams 18 (102801) (2015).

[9] D. P. Kingma, M. Welling, An introduction to variational autoencoders, Foundations and Trends® in Machine Learning 12 (4) (2019) 307–392. doi:10. 1561/2200000056. URL http://dx.doi.org/10.1561/2200000056

[3] M. Noack, K. Yager, M. Fukuto, G. Doerk, L. Ruipeng, J. Sethian, A kriging-based approach to autonomous experimentation with applications to xray scattering, Scientific Reports 9 (11809) (2019).

[10] O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation (2015). arXiv:1505.04597. URL https://arxiv.org/abs/1505.04597

[4] A. Diaw, M. McKerns, I. Sagert, L. G. Stanton, M. S. Murillo, Efficient learning of accurate surrogates for simulations of complex systems, Nature Machine Intelligence 6 (2024) 568–577. doi:10.1038/ s42256-024-00839-1.

[11] J. Degrave, F. Felici, J. Buchli, M. Neunert, B. Tracey, F. Carpanese, T. Ewalds, R. Hafner, 15

Figure C.9: Predicted electron profiles as a function of normalized poloidal flux ψN for a representative held-out test case. Solid lines: SOLPS-ITER ground truth. Dashed lines (“NN forward”): the forward surrogate Fθ evaluated at the true control parameters, isolating forward-model accuracy. Dotted lines (“NN inverse cycle”): the cycle reconstruction, in which the control parameters are first recovered from the SOLPS-ITER fields by the inverse procedure G⋆ and then propagated back through the forward model. Panels show (a) upstream (outboard-midplane) ne , (b) upstream Te , (c) outer-target Te , and (d) inner-target Te . The inverse-cycle profiles remain close to both the forward prediction and SOLPS-ITER, demonstrating cycle consistency: parameters recovered by the inverse procedure reproduce the original fields.

A. Abdolmaleki, D. de Las Casas, et al., Magnetic control of tokamak plasmas through deep reinforcement learning, Nature 602 (2022) 414–419. doi: 10.1038/s41586-021-04301-9.

[14] K. Willcox, B. Segundo, The role of computational science in digital twins, Nature Computational Science 4 (2024) 147–149. doi:10.1038/ s43588-024-00609-4.

[12] J. Seo, S. Kim, A. Jalalvand, R. Conlin, A. Rothstein, J. Abbate, K. Erickson, J. Wai, R. Shousha, E. Kolemen, Avoiding fusion plasma tearing instability with deep reinforcement learning, Nature 626 (2024) 746– 751. doi:10.1038/s41586-024-07024-9.

[15] A. Diaw, K. Barros, J. Haack, C. Junghans, B. Keenan, Y. W. Li, D. Livescu, N. Lubbers, M. McKerns, R. S. Pavel, D. Rosenberger, I. Sagert, T. C. Germann, Multiscale simulation of plasma flows using active learning, Phys. Rev. E 102 (2) (2020) 023310. doi:10.1103/PhysRevE.102.023310.

[13] W. Tang, E. Feibush, G. Dong, N. Borthwick, AImachine learning-enabled tokamak digital twin, arXiv preprint (2024). arXiv:2409.03112.

[16] A. Diaw, M. McKerns, I. Sagert, L. G. Stanton, M. S. Murillo, Efficient Learning of Accurate Surrogates for Simulations of Complex Systems, arXiv 16

e-prints (2022) arXiv:2207.12855arXiv:2207.12855, doi:10.48550/arXiv.2207.12855.

[25] B. Zhu, M. Zhau, X. Xu, A. Gupta, K. B. Kwon, X. Ma, D. Eldon, Latent space mapping: Revolutionizing predictive models for divertor plasma detachment control, Physics of Plasmas 32 (062508) (2025) 1–19. doi:http://dx.doi.org/10.1063/5. 0267930.

[17] S. Wiesen, D. Reiter, V. Kotov, M. Baelmans, W. Dekeyser, A. S. Kukushkin, S. W. Lisgo, R. A. Pitts, V. Rozhansky, G. Saibene, I. Veselova, S. Voskoboynikov, The new solps-iter code package, Journal of Nuclear Materials 463 (2015) 480– 484. doi:http://dx.doi.org/10.1016/j.jnucmat. 2014.10.012.

[26] J. Bentley, Multidimensional binary search trees used for associative searching, Communications of the ACM 18 (9) (1975) 509–517. doi:http://dx.doi. org/10.1145/361002.361007.

[18] D. Reiter, M. Baelmans, P. Borner, The eirene and b2-eirene codes, Fusion Science and Technology 47 (2) (2005) 172–186. doi:http://dx.doi.org/10. 13182/FST47-172.

[27] A. Krizhevsky, I. Sutskever, G. E. Hinton, Imagenet classification with deep convolutional neural networks, Advances in neural information processing systems 25 (2012).

[19] J. D. Lore, S. De Pascuale, P. Laiu, B. Russo, J. S. Park, J. M. Park, S. L. Brunton, J. N. Kutz, A. A. Kaptanoglu, Time-dependent solps-iter simulations of the the tokamak plasma boundary for model predictive control using sindy, Nuclear Fusion 63 (046015) (2023) 1–12. doi:http://dx.doi.org/ 10.1088/1741-4326/acbe0e.

[28] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, L. D. Jackel, Backpropagation applied to handwritten zip code recognition, Neural computation 1 (4) (1989) 541–551. [29] E. Perez, F. Strub, H. de Vries, V. Dumoulin, A. Courville, FiLM: Visual reasoning with a general conditioning layer, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 32, 2018, pp. 3942–3948. doi:10.1609/aaai.v32i1.11671. URL https://arxiv.org/abs/1709.07871

[20] J. S. Park, J. D. Lore, M. Reinke, A. Q. Kuang, S. De Pascuale, A. Creely, Full time-dependent solps-iter simulation of the sparc tokamak: actuator design for particle and divertor condition control, Nuclear Fusion 64 (4) (2024) 1–14. doi:http://dx.doi.org/ 10.1088/1741-4326/ad4f9d.

[30] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization (2017). arXiv:1412.6980. URL https://arxiv.org/abs/1412.6980

[21] A. Lasa, J. S. Park, J. Lore, S. Blondel, D. E. Bernholdt, J. M. Canik, M. Cianciosa, J. Coburn, D. Curreli, Exploreing the efect of elm and code-coupling frequencies on plasma and material modeling of dynamic recycling in divertors, Nuclear Fusion 64 (7) (2024) 1– 13. doi:http://dx.doi.org/10.1088/1741-4326/ ad4894.

[31] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala, Pytorch: An imperative style, high-performance deep learning library (2019). arXiv:1912.01703. URL https://arxiv.org/abs/1912.01703

[22] S. Dasbach, S. Wiesen, Towards fast surrogate models for interpolation of tokamak edge plasmas, Nuclear Materials and Energy 34 (101396) (2023) 1–6. doi: http://dx.doi.org/10.1016/j.nme.2023.101396.

[32] B. Kustowski, J. A. Gaffney, B. K. Spears, G. J. Anderson, J. J. Thiagarajan, R. Anirudh, Erratum to “Transfer Learning as a Tool for Reducing Simulation Bias: Application to Inertial Confinement Fusion” [Jan 20 46-53], IEEE Transactions on Plasma Science 48 (9) (2020) 3275–3275. doi:10.1109/TPS. 2020.3012229.

[23] S.Wiesen, S. Dasbach, A. Kit, A. E. Jaervinen, A. Gillgren, A. Ho, A. Panera, D. Reiser, M. Brenzke, Y. Poels, E. Westerhof, V. Menkovski, G. F. Derks, P. Strand, Data-driven models in fusion exhaust: Ai methods and perspectives, Nuclear Fusion 64 (086046) (2024) 1–10. doi:http://dx.doi.org/ 10.1088/1741-4326/ad5a1d.

[33] R. Anirudh, J. J. Thiagarajan, P.-T. Bremer, B. K. Spears, Improved surrogates in inertial confinement fusion with manifold and cycle consistencies, Proceedings of the National Academy of Sciences 117 (18) (2020) 9741–9746. doi:10.1073/pnas.1916634117.

[24] B. Zhu, M. Zhao, H. Bhatia, X. Xu, P. Bremer, W. Meyer, N. Li, T. Rognlien, Data-driven model for divertor plasma detachment prediction, Journal of Plasma Physics 88 (895880504) (2022) 1–23. doi: http://dx.doi.org/10.1017/S002237782200085X.

[34] J. St. John, C. Herwig, D. Kafkes, J. Mitrevski, W. A. Pellico, G. N. Perdue, A. Quintero-Parra, B. A. Schupbach, K. Seiya, N. Tran, M. Schram, 17

J. M. Duarte, Y. Huang, R. Keller, Real-time artificial intelligence for accelerator control: A study at the Fermilab Booster, Physical Review Accelerators and Beams 24 (2021) 104601. doi:10.1103/ PhysRevAccelBeams.24.104601.

[44] C. Ticknor, J. D. Kress, L. A. Collins, J. Clérouin, P. Arnault, A. Decoster, Transport properties of an asymmetric mixture in the dense plasma regime, Phys. Rev. E 93 (6) (2016) 063208. doi:10.1103/ PhysRevE.93.063208.

[35] H. S. Seung, M. Opper, H. Sompolinsky, Query by committee, in: Proceedings of the Fifth Annual Workshop on Computational Learning Theory, COLT ’92, Association for Computing Machinery, New York, NY, USA, 1992, p. 287–294. doi:10.1145/130385. 130417. URL https://doi.org/10.1145/130385.130417

[45] H. Flyvbjerg, H. G. Petersen, Error estimates on averages of correlated data, The Journal of Chemical Physics 91 (1989) 461–466. doi:10.1063/1.457480. [46] A. D. Sokal, Monte Carlo methods in statistical mechanics: Foundations and new algorithms, Lecture Notes in Physics 490 (1997) 131–192.

[36] R. S. Wilcox, M. W. Shafer, J. D. Lore, J. M. Canik, S. R. Haskey, C. J. Lasnier, A. L. Moser, T. H. Osborne, H. Q. Wang, Challenges and approaches to interpretive modeling of boundary plasma and neutral transport in a closed, pumped divertor, Nuclear Fusion 66 (2026) 016023. doi:10.1088/1741-4326/ ae1bb9. [37] R. Roussel, A. Edelen, C. Mayes, D. Ratner, J. P. Gonzalez-Aguilera, S. Kim, E. Wisniewski, J. Power, Phase space reconstruction from accelerator beam measurements using neural networks and differentiable simulations, Physical Review Letters 130 (2023) 145001. doi:10.1103/PhysRevLett.130.145001. [38] E. L. Allgower, K. Georg, Introduction to Numerical Continuation Methods, Vol. 45 of Classics in Applied Mathematics, SIAM, 2003. doi:10.1137/1. 9780898719154. [39] M. D. McKay, R. J. Beckman, W. J. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21 (2) (1979) 239–245. doi:10.1080/00401706.1979.10489755. [40] J. H. Friedman, J. L. Bentley, R. A. Finkel, An algorithm for finding best matches in logarithmic expected time, ACM Transactions on Mathematical Software 3 (3) (1977) 209–226. doi:10.1145/ 355744.355745. [41] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, 2nd Edition, Oxford University Press, 2017. doi:10.1093/oso/9780198803195.001.0001. [42] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd Edition, Academic Press, 2002. [43] A. Diaw, M. S. Murillo, Generalized hydrodynamics model for strongly coupled plasmas, Phys. Rev. E 92 (1) (2015) 013107. doi:10.1103/PhysRevE.92. 013107.

18