Paper deep dive
Cycle-Consistent and Uncertainty-Aware Neural Surrogates for Tokamak Edge Plasmas
Abdourahmane Diaw, Sebastian De Pascuale, Jae-Sun Park, Ivan Paradela Perez, Jeremy D. Lore, Stefan Dasbach
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 7/24/2026, 3:10:14 AM
Summary
The paper introduces a cycle-consistent neural surrogate for predicting tokamak edge plasma states, combining a conditional U-Net forward model with an optimization-based inverse method. This approach enables fast, accurate prediction of 2D plasma fields from control parameters and allows for the recovery of input parameters from observed fields via cycle consistency. Additionally, an ensemble of multilayer perceptrons provides uncertainty estimates for 1D profiles, facilitating real-time control and active learning. The model achieves high accuracy and speed, significantly outperforming traditional SOLPS-ITER simulations.
Entities (8)
Relation Signals (5)
Conditional U-Net → usedfor → Edge Plasma Prediction
confidence 97% · The forward model maps five control parameters to two-dimensional plasma-state fields
SOLPS-ITER → isslowerthan → Neural Surrogate
confidence 96% · five to six orders of magnitude faster than SOLPS-ITER
Cycle Consistency → enables → Inverse Parameter Recovery
confidence 95% · The inverse method enforces consistency between forward and inverse predictions... enables recovery of the core fueling rate
Multilayer Perceptrons → provides → uncertainty estimates
confidence 94% · An ensemble of multilayer perceptrons also predicts electron temperature and density profiles... with uncertainty estimates
DIII-D → sourceof → Training Data
confidence 92% · The training dataset was generated using SOLPS-ITER simulations of a deuterium plasma in DIII-D
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The boundary and divertor plasma govern how a tokamak exhausts power and particles, setting heat fluxes, target conditions, and the onset of detachment. Predicting these quantities is essential for operating current and future devices, but edge simulations that resolve them are too slow for parameter scans, optimization, or real-time control. Machine-learning surrogates offer a fast alternative, yet most are forward-only: they cannot recover input parameters from observations or assess the reliability of their predictions. We introduce a cycle-consistent neural surrogate for edge plasmas, combining a conditional U-Net forward model with an optimization-based inverse method built on the frozen forward network. The forward model maps five control parameters to two-dimensional plasma-state fields on the SOLPS-ITER mesh; the inverse method enforces consistency between forward and inverse predictions, a self-supervised quality check needing no ground-truth labels at inference. An ensemble of multilayer perceptrons also predicts electron temperature and density profiles at the outboard midplane and divertor targets, with uncertainty estimates that flag 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 fields. Cycle-consistency regularization raises the average cyclical $R^2$ from 0.59 to 0.99 without degrading forward accuracy and enables recovery of the core fueling rate; all five control parameters are recovered with Pearson $r\ge0.97$. A $k$-d tree warm start yields a database completion rate above 95%, versus roughly 30% outright failures when cold-started. With about $4\times10^6$ parameters, the model produces full 2D predictions in milliseconds, five to six orders of magnitude faster than SOLPS-ITER, enabling real-time control, parameter scans, uncertainty analysis, and digital twins.
Tags
Links
- Source: https://arxiv.org/abs/2607.21407v1
- Canonical: https://arxiv.org/abs/2607.21407v1
Trouble viewing inline? Open PDF directly →
Full Text
74,616 characters extracted from source content.
Expand or collapse full text
Cycle-Consistent and Uncertainty-Aware Neural Surrogates for Tokamak Edge Plasmas A. Diaw diawa@ornl.gov S. De Pascuale J.-S. Park I. Paradela Perez J.D. Lore S. Dasbach Fusion Energy Division, Oak Ridge National Laboratory, Oak Ridge, TN 37932, USA DIFFER - Dutch Institute for Fundamental Energy Research, De Zaale 20, 5612 AJ Eindhoven, the Netherlands Abstract 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 ground-truth 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%2.6\% and Pearson correlations above 0.950.95 for all plasma-state fields. Adding cycle-consistency regularization raises the average cyclical R2R^2 from 0.590.59 to 0.990.99 without reducing forward accuracy, and it allows recovery of the core fueling rate Γcore _core, which is hard to determine with forward-only training. The inverse method recovers all five control parameters with Pearson r≥0.97r≥ 0.97. Using a k-d tree warm start to build the database yields a completion rate above 95%95\%, whereas in a comparable cold-started ensemble roughly 30%30\% of runs failed outright and a further fraction never reached steady state. With about 4×1064× 10^6 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 †journal: Nuclear Materials and Energynoticenoticefootnotetext: 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). 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 [42, 33, 23, 9]. 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 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 [34, 27]. For SOL physics, the state-of-the-art tool is SOLPS-ITER, which couples a two-dimensional Braginskii fluid solver (B2.5) to a kinetic Monte Carlo neutral transport model (EIRENE) on a triangulated vessel mesh [28]. 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 fixed-point 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 [14]. 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θF_θ that maps five scalar control parameters to the two-dimensional plasma state (TeT_e, TiT_i, nen_e, uau_a). 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θ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 U-Net. 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 [16, 30]. 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 [6, 35] and digital-twin workflows [39, 44]. Alongside this two-dimensional model, we train a query-by-committee ensemble of multilayer perceptrons for the one-dimensional TeT_e and nen_e profiles at the outboard midplane and the divertor targets. We restrict it to TeT_e and nen_e as an initial demonstration of the method, and because these are the quantities measured directly by Thomson scattering, facilitating future validation. The ensemble supplies predictive uncertainty estimates for closed-loop control, and its committee variance flags the under-sampled regimes where new SOLPS-ITER runs are most needed, a query-by-committee strategy we have used in earlier work [7, 8]. 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 warm-start 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 uncertainty-quantified 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 [41, 29] simulations of a deuterium plasma in DIII-D, employing the same lower-single-null configuration as used in Ref. [21]. This setup was developed using DIII-D discharge 174310 at t=3500t=3500\,ms, has been demonstrated to have excellent numerical stability, and has been used to assess divertor conditions and neutral dynamics [21, 24, 19]. 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 [5, 32, 45, 46]. 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: =(Ptot,Γcore,ΓD2,D⟂,χi),c=(P_tot, _core, _D_2,D_ , _i), which sets the total power across the separatrix, the target core flux, the deuterium puff rate, and the uniform cross-field 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 [5], about 29%29\% of runs (11981198 of 40964096) diverged outright, and a further fraction of the surviving runs never reached steady state within the allotted runtime. Table 1: Control parameters and their sampling ranges for the SOLPS-ITER database. PtotP_tot Γcore _core ΓD2 _D_2 D⟂D_ χi _i (MW) (at/s) (at/s) (m2/s) (m2/s) min 2 1×10201× 10^20 1.3×10201.3× 10^20 0.1 0.1 max 16 7.5×10207.5× 10^20 5×10215× 10^21 2.0 2.0 To mitigate these issues, we adopt a branching k-d tree strategy [4] that both partitions the input space and warm-starts each simulation from the nearest already-converged design point, propagating reliable initial conditions outward from a single seed run. This raises the completion rate above 95%95\%, compared with the roughly 70%70\% completion of the cold-started ensemble above. The construction and distance metric are detailed in Appendix A. The detailed convergence detection procedure, including autocorrelation analysis and block averaging, is described in Appendix B. The completed database comprises 762762 converged runs, of which two were discarded during quality screening, leaving 760760 for surrogate modeling. We split these 80/2080/20 into a training/validation set of 608608 runs and a held-out test set of 152152 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. 3 Surrogate Model Design The surrogate framework follows the manifold and cycle-consistency paradigm. At inference, it runs as a three-stage cycle (Fig. 1): a forward model, an optimization-based inverse, and a cycle-consistency self-check. All three stages share a single trained forward model FθF_θ and a learned pseudo-inverse Gψ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 U-Net FθF_θ is trained on the SOLPS-ITER database to map scalar control parameters =(Ptot,Γcore,ΓD2,D⟂,χi)c=(P_tot,\, _core,\, _D_2,\,D_ ,\, _i) and a spatial binary mask m to the four-channel two-dimensional plasma state. Once trained, FθF_θ defines the learned physics manifold: every output ^=Fθ(,m) y=F_θ(c,m) is, by construction, a physically plausible edge-plasma state, (2) an inverse model: given a target field ∗y^*, an optimization-based inverse procedure G⋆G recovers the control parameters c by minimizing the discrepancy ‖Fθ(^,m)−∗‖\|F_θ( c,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θF_θ to obtain a reconstructed field ^cycle=Fθ(G⋆(∗),m) y_cycle=F_θ\! (G (y^*),m ). The cycle-consistency loss ‖^cycle−∗‖\| y_cycle-y^*\| measures whether the round-trip ∗→^→^cycley^*\!→\! c\!→\! y_cycle preserves the target, providing a self-supervised quality metric that does not require ground-truth parameter labels at test time. Each stage is described in detail in the following subsections. Figure 1: Surrogate model architecture: the three-stage inference cycle. (a) Stage 1 (forward model): a conditional U-Net maps a single-channel geometry mask, with the five scalar control parameters injected via FiLM, to the four-channel plasma state (Te,Ti,ne,uaT_e,T_i,n_e,u_a); the fields 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ψ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θF_θ and pseudo-inverse GψG_ψ used here are produced beforehand by the training-time regularizer of Fig. 2. 3.1 Conditional U-Net Forward Model The forward surrogate FθF_θ maps scalar control parameters c and the spatial binary mask m to the multi-channel plasma state y = = Fθ(,m), F_θ(c,\,m), (1) y = = Te,Ti,ne,ua, \T_e,\,T_i,\,n_e,\,u_a\, (2) on the SOLPS (R,Z)(R,Z) grid, where TeT_e and TiT_i are the electron and ion temperatures, nen_e is the electron density, and uau_a is the parallel velocity of the deuterium ions. Although the SOLPS state arrays are dense in index space (i,j)(i,j), the physical plasma domain is irregular in (R,Z)(R,Z). It occupies only a subset of the rectangular tensor used by convolutional neural networks (CNNs) [17, 20]. Standard U-Net 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) [26]: 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 single-channel 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 log-transformed before being centered and scaled to unit variance (with a small offset, 10−210^-2 for Te,TiT_e,T_i and 1016m−310^16\,m^-3 for nen_e, to regularize near-zero values). The parallel velocity uau_a 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 [30] with three resolution levels. At each level, an input feature map is processed by two 3×33× 3 convolutions, each followed by group normalization and a SiLU activation: ′ = \;=\; σ(GN(∗)), σ\! (GN(W*Y) ), (3) σ(s) σ(s) = = s⋅sigmoid(s). s·sigmoid(s). (4) Downsampling is performed using 2×22× 2 max pooling with a stride of 2; at each step, the spatial size halves while the number of channels doubles. The bottleneck applies the same double-convolution block at 8×8× base channels. Decoder. Each up-step uses a 2×22× 2 transposed convolution to increase spatial resolution and reduce channels, concatenates the corresponding encoder feature map via a skip connection, and then applies two 3×33× 3 convolutions with group normalization and SiLU. After the last decoder block, a 1×11× 1 convolution maps the base channels to Cout=4C_out=4 output channels corresponding to the four plasma state fields. Because the targets are continuous, we use a linear activation at the final layer. Additionally, because ne≡nin_e≡ n_i in the plasma, we omit the ion density channel. For thermodynamic consistency, we enforce non-negativity of the temperature and density outputs (Te,Ti,ne≥0T_e,T_i,n_e≥ 0) by clamping negative values during post-processing. 3.2 Inverse Model The inverse problem, in which we infer control parameters c from target plasma fields ∗y^*, is ill-posed: many 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θF_θ. We denote this optimization-based inverse G⋆G ; unlike FθF_θ and the pseudo-inverse GψG_ψ (Sec. 3.4), it has no trainable parameters of its own: it recovers c 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 c that minimize the discrepancy between the forward prediction and the target: c = = argmin′ℒinv(′), _c \;L_inv(c ), ℒinv(′) _inv(c ) = = ∑c,i,jwcmij(Fθ(′,m)c,i,j−yc,i,j∗)2∑c,i,jwcmij _c,i,jw_c\,m_ij (F_θ(c ,m)_c,i,j-y^*_c,i,j )^\!2 _c,i,jw_c\,m_ij (5) +λreg‖′‖2, + _reg\|c \|^2, where wcw_c are the per-channel weights and λreg _reg is an L2L_2 regularization coefficient that prevents the recovered parameters from drifting to unphysical extremes. Optimization is performed using Adam [15] with learning rate 10−210^-2 for 12001200 steps through the frozen forward model, with gradient-based parameter updates computed via automatic differentiation [25]. The optimization is warm-started from the pseudo-inverse estimate ^0=Gψ(∗,m) c_0=G_ψ(y^*,m) (Sec. 3.4), which is computed from the observed fields alone. To mitigate local minima, we employ Nr=5N_r=5 random restarts about this estimate, perturbed with Gaussian noise (σnoise=0.2 _noise=0.2 in standardized space), and retain the solution with the lowest ℒinvL_inv; no ground-truth parameters are used at any point. 3.3 Cycle Consistency The cycle-consistency constraint links the forward and inverse models by requiring that parameters recovered from a target field, when fed back through the forward model, reproduce that target field. Formally, for a ground-truth field ∗y^* with mask m: ℒcycle=‖Fθ(G⋆(∗),m)−∗‖,L_cycle= \|F_θ\! (G (y^*),\,m )-y^* \|, (6) where G⋆(∗)G (y^*) denotes the inverse-recovered parameters (Eq. 5) and ∥⋅∥\|·\| is a masked norm over valid mesh cells. This constraint serves several purposes. First, it regularizes the inverse mapping: even when the forward model’s parameter-to-field mapping is locally flat (so that many parameter vectors produce similar fields), the cycle loss penalizes inverse solutions whose forward reconstructions deviate from the target. Second, it enforces manifold consistency [18, 3]: the reconstructed fields ^cycle=Fθ(^,m) y_cycle=F_θ( c,m) are guaranteed to lie on the learned physics manifold of the forward model, ensuring physically plausible outputs. Third, the cycle loss provides a self-supervised quality metric for the inverse procedure that does not require ground-truth parameter labels at test time. 3.4 Cycle Consistency as a Training Regularizer Figure 2: Training-time cycle-consistency regularizer. The forward model FθF_θ and the learned pseudo-inverse GψG_ψ are trained jointly; the cycle penalty λcyc‖−Gψ(Fθ(,m),m)‖ _cyc\|c-G_ψ(F_θ(c,m),m)\| back-propagates into both networks, so consistency reshapes the forward surrogate. Setting λcyc=0 _cyc=0 recovers the forward-only baseline; the weight is swept in the ablation of Sec. 4.3. 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 [3], we introduce a learned pseudo-inverse Gψ:↦G_ψ:y and train it jointly with the forward model under minθ,ψρ(Fθ(,m),)⏟forward+ρ(Gψ(,m),)⏟inverse+λcyc‖−Gψ(Fθ(,m),m)‖⏟cycle, split _θ,ψ\;& ρ\! (F_θ(c,m),y )_forward+ ρ\! (G_ψ(y,m),c )_inverse\\ &+ _cyc\, \|c-G_ψ\! (F_θ(c,m),m ) \|_cycle, split (7) where ρ is a masked smooth-L1L_1 discrepancy. Unlike Eq. 6, the cycle term here back-propagates into the forward parameters θ, so consistency can reshape the surrogate. The weight λcyc _cyc multiplies only the cycle term; the forward and inverse reconstruction terms are kept at unit weight, so λcyc=0 _cyc=0 recovers the forward-only baseline. The pseudo-inverse Gψ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 2424 to 4848 to 9696; each stage applies two 3×33× 3 convolutions, each followed by group normalization and a SiLU activation, and then a 2×22× 2 max-pooling that halves the spatial resolution. A global average pooling collapses the final feature map to a 9696-dimensional vector, which a two-layer perceptron (hidden width 128128) maps to the five control parameters. The 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 R2R^2 score analogous to the metric of Anirudh et al. [3]: for each control parameter, we sweep it linearly across its range (holding the others fixed), push the resulting parameter vectors through FθF_θ and then Gψ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. 3.5 Loss Functions and Training The forward model FθF_θ is trained on the SOLPS-ITER database using a composite loss function that balances pixel-level accuracy, gradient fidelity, and boundary emphasis: ℒ=ℒbase+λwℒedge+λgℒgrad.L=L_base+ _wL_edge+ _gL_grad. (8) The base loss ℒbaseL_base is a masked Huber (smooth-L1L_1) loss with β=0.05β=0.05, computed only over valid mesh cells indicated by the binary mask; per-channel weights wcw_c allow prioritization of specific output fields, and we use wc=1.0w_c=1.0 for Te,TiT_e,T_i, wc=1.2w_c=1.2 for nen_e, and wc=1.5w_c=1.5 for uau_a. The edge loss ℒedgeL_edge is the same Huber loss reweighted by a boundary proximity map w()w(r), computed as a Gaussian-decayed distance from the mask boundary (σ=3σ=3 pixels); this emphasizes accuracy near the separatrix and target plates where gradients are steepest. Finally, ℒgradL_grad is a Sobel-filter-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 L1L_1 difference of the resulting gradient maps. We set λg=0.2 _g=0.2 and apply a linear warmup from epoch 20 to 60 to stabilize early training. Training uses the Adam optimizer [15] with initial learning rate 3×10−43× 10^-4, ReduceLROnPlateau scheduling (factor 0.50.5, patience 5 epochs), mixed-precision (AMP) on GPU, gradient clipping at norm 1.01.0, and early stopping with patience of 80 epochs. The model is trained for up to 450450 epochs with a batch size of 88, 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 yiy_i for the SOLPS-ITER value and y^i y_i for the prediction at cell i (over N valid cells), the mean absolute error is MAE=N−1∑i|y^i−yi|MAE=N^-1 _i| y_i-y_i|, and the normalized root-mean-square error is NRMSE=100×RMSE/(ymax−ymin)NRMSE=100×RMSE/(y_ -y_ ), 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 R2=1−∑i(y^i−yi)2/∑i(yi−y¯)2, R^2=1- _i( y_i-y_i)^2/ _i(y_i- y)^2, (9) where y¯ y 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, ‖Fθ()−Fθ(+σϵ)‖2E\,\|F_θ(c)-F_θ(c+σ ε)\|^2 in normalized output space (Table 3), for which a smaller value indicates a smoother, more robust surrogate. 4.1 Forward-model accuracy Figure 3 compares SOLPS-ITER ground truth with U-Net predictions for the four plasma-state fields (TeT_e, TiT_i, nen_e, uau_a) 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 scrape-off layer. Electron density is well captured across the full domain, with the core–SOL contrast and divertor compression faithfully represented. The parallel velocity field uau_a shows correct flow patterns from the outer midplane toward the divertor targets. 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 TeT_e (eV), ion temperature TiT_i (eV), electron density nen_e (m-3), and parallel velocity uau_a (m/s). The model captures the large-scale structure and gradients across the SOL and divertor. Quantitatively, the forward model achieves a global Pearson correlation exceeding 0.950.95 for all four plasma fields: r=0.994r=0.994 for TeT_e, r=0.992r=0.992 for TiT_i, r=0.987r=0.987 for nen_e, and r=0.975r=0.975 for uau_a (Table 2). Mean absolute errors are 12.212.2 eV for TeT_e and 16.716.7 eV for TiT_i, which represent small fractions of the dynamic range of these fields across the database. The density field achieves an MAE of 2.6×10182.6× 10^18 m-3 for nen_e, while the parallel velocity MAE is 941941 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. 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−yminy_ -y_ ) 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 (%) Pearson r Spearman ρ TeT_e eV 12.24 1.26 0.994 0.996 TiT_i eV 16.71 1.54 0.992 0.987 nen_e m-3 2.61×10182.61× 10^18 0.62 0.987 0.964 uau_a m⋅·s-1 940.9 2.54 0.975 0.976 The four plasma state fields (TeT_e, TiT_i, nen_e, uau_a) all achieve Pearson correlations above 0.950.95 and NRMSE below 2.6%2.6\%, with per-sample mean Pearson values exceeding 0.980.98. The parallel velocity uau_a is the hardest of the four to reproduce (r=0.975r=0.975, NRMSE 2.54%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.2 Model size and inference speed The conditional U-Net has ∼4.3×106 4.3× 10^6 trainable parameters with a base filter width of 48, occupying 16.516.5 MB in single precision. Each forward pass maps the single-channel mask input, with the five scalar parameters injected via FiLM, to four output channels at 36×9636× 96 resolution. On a single CPU core, inference takes 1616 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 (11–1010 ms cycle times) [6, 35, 38], digital-twin frameworks [39, 44], 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 _cyc in the joint objective (Eq. 7) and, for each value, retrain the forward model with the learned pseudo-inverse from scratch. Figure 4(A) and Table 3 report the held-out test mean squared error and the average cyclical R2R^2 score as λcyc _cyc is increased, following the ablation protocol of Anirudh et al. [3]. As λcyc _cyc increases, the forward and inverse networks become markedly more cyclically self-consistent: the average cyclical R2R^2 increases from 0.590.59 at λcyc=0 _cyc=0 to 0.990.99 at λcyc=0.5 _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 _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 R2R^2 of the forward-only baseline with that of the consistency-regularized model (λcyc=0.5 _cyc=0.5). Without the cycle term, the trained pair is self-consistent for some control parameters but not others: the recovered R2R^2 for the core particle source Γcore _core and the cross-field diffusivity D⟂D_ is low, whereas the input power PtotP_tot and the gas-puff rate ΓD2 _D_2 are already well recovered. Training with the cycle term makes the pair self-consistent across all five parameters, each above 0.970.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 _cyc=0.5, which is both the most self-consistent and the most accurate on held-out data. Figure 4: Cycle-consistency ablation. (A) Average cyclical R2R^2 (left axis, higher is better) and held-out test MSE (right axis, lower is better) as the cycle weight λcyc _cyc is increased. Self-consistency increases from 0.590.59 at the forward-only baseline (λcyc=0 _cyc=0) to 0.990.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) Per-parameter self-consistency of the trained forward/inverse pair for the forward-only baseline (λcyc=0 _cyc=0, orange) versus the cycle-consistent model (λcyc=0.5 _cyc=0.5, blue). The cycle term brings every control parameter above 0.970.97, including Γcore _core and D⟂D_ , which are otherwise weakly recovered by the pair. Higher R2R^2 is better; lower MSE is better. Each λcyc _cyc is trained from scratch, so the intermediate points carry run-to-run variation; the baseline (λcyc=0 _cyc=0) and the adopted λcyc=0.5 _cyc=0.5 are the comparison of interest. Exact sweep values are given in Table 3. Table 3: Cycle-consistency ablation sweep (760-run database; 608 training / 152 test runs, trained from scratch at each λcyc _cyc). As the cycle weight increases, the average cyclical R2R^2 increases toward unity while the held-out test MSE stays comparable to or below the forward-only baseline (λcyc=0 _cyc=0). Sensitivity is the mean squared change of the normalized prediction under Gaussian input perturbations (σ=0.1σ=0.1); a smaller value indicates a more robust surrogate. Best values in bold. λcyc _cyc Test MSE Cyclical R2R^2 Sens. (×10−3× 10^-3) 0 (baseline) 0.02880.0288 0.5850.585 2.142.14 0.010.01 0.02720.0272 0.7600.760 3.933.93 0.050.05 0.02740.0274 0.8950.895 3.973.97 0.10.1 0.02870.0287 0.9750.975 5.105.10 0.50.5 0.02660.0266 0.9890.989 3.513.51 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 pseudo-inverse Gψ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ψ(∗,m) c_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θF_θ, with Nr=5N_r=5 random restarts per case and cosine learning-rate annealing over 1200 Adam steps. The cycle-consistency metric ℒcycleL_cycle (Eq. 6) is then evaluated by passing the recovered parameters through FθF_θ and comparing the reconstructed fields against the original targets. This procedure provides two complementary diagnostics: (i) the accuracy of parameter recovery (how close c is to the true ∗c^*), and (i) the cycle reconstruction quality (how well Fθ(^,m)F_θ( c,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 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.970.97: PtotP_tot (r=0.99r=0.99, ρ=0.99ρ=0.99), ΓD2 _D_2 (r=0.99r=0.99, ρ=0.99ρ=0.99), D⟂D_ (r=0.99r=0.99, ρ=0.97ρ=0.97), χi _i (r=0.99r=0.99, ρ=0.97ρ=0.97), and the core fueling Γcore _core (r=0.97r=0.97, ρ=0.95ρ=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. 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. 5 Uncertainty-Quantified Profile Ensemble for Control and Active Learning 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 one-dimensional profiles of electron temperature TeT_e and density nen_e at the outboard midplane and both divertor targets directly from the five control parameters, as a function of the normalized poloidal flux ψN _N. We restrict the targets to TeT_e and nen_e 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 they are constitutive functions of TeT_e and nen_e 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 [5, 45]. 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-by-committee (QBC) approach [36], 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 [7, 8]. Training minimizes mean-squared error loss, ℒ=∑i,t(y^i(t)−yi(t))2L\;=\; _i,t\, ( y^(t)_i-y^(t)_i )^2 using Adam [15] 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, with 10%10\% held out for early-stopping validation and a further held-out fold used to score it. Members that do not reach R2≥0.90R^2≥ 0.90 on every output channel of that fold are discarded. We retain nensemble=5n_ ensemble=5 trained models. Because nen_e spans several decades and the target TeT_e ranges from below 11 to above 200200 eV, both are trained on their log10 _10 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 si=kσiS,s_i\;=\;k\, _iS, (10) where σi _i is the ensemble standard deviation at point i and k=2k=2. The global scale S=σ¯E¯cal/E¯ensS= σ\, E_ cal/ E_ ens rescales the raw committee spread to the magnitude of the actual prediction error: σ¯ σ is the mean ensemble standard deviation, E¯cal E_ cal the mean absolute error of the ensemble-mean prediction, and E¯ens E_ ens the mean ensemble spread, all computed on the held-out validation fold. If si≥1s_i≥ 1, the point is flagged for further simulation. We train each ensemble on 1216012160 points (608608 SOLPS-ITER runs) and evaluate it on 30403040 held-out test points (152152 runs), using the same 80/2080/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 TeT_e and nen_e profiles with high fidelity, with the per-panel coefficient of determination annotated in Fig. 6. Points are colored by the quality flag sis_i; those selected for SOLPS-ITER verification (si>1s_i>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 TeT_e coefficient of determination (R2≈0.89R^2≈ 0.89, versus ≳0.94 0.94 at the upstream and outer-target locations) and by far the largest fraction of flagged points (∼13% 13\%, against ∼3% 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 TeT_e 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. 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 nen_e and electron temperature TeT_e, on the independent test set (30403040 points per location; 1216012160 training points). Target TeT_e 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 sis_i; open circles mark the highest-uncertainty points (si>1s_i>1) and are flagged for additional SOLPS-ITER simulation. 6 Discussion Our framework builds on the manifold-and-cycle-consistency paradigm for surrogate modeling of tokamak edge plasmas. Anirudh et al. [3] enforce manifold consistency with a Wasserstein autoencoder that learns a low-dimensional latent representation; we instead use a conditional U-Net whose skip connections preserve fine spatial detail. Both share the core insight that cycle consistency regularizes the inverse model, and in both, the reconstructed outputs are constrained to lie on the learned physics manifold. The tokamak setting adds its own difficulties: the mesh geometry is irregular, which is what forces the masking described above, and the fields carry hundreds of sharp recycling fronts that do not admit a smooth, low-dimensional structure. Table 2 shows that all four fields are reproduced accurately (Pearson correlation >0.95>0.95), with a mild hierarchy among them: the temperatures and density, which vary smoothly over the scale of the scrape-off layer (SOL) width, are captured most efficiently by the convolutional receptive field, whereas the parallel velocity uau_a carries sharper spatial structure and a sign change along the flow and is therefore the hardest to fit. That same spatial structure is what motivates the convolutional inductive bias: coherent large-scale flow patterns like uau_a are reproduced far more accurately by a model that shares information between neighboring cells than by a pointwise regressor. This supports the choice of a deterministic, parameter-conditioned convolutional operator for the data-scarce regime studied here, where generative alternatives would have to learn a full distribution from only a few hundred samples. The millisecond-scale inference described in Sec. 4.2 enables several practical applications: (i) rapid exploration of the five-dimensional parameter space for sensitivity analysis and scenario optimization [7, 23]; (i) transport-coefficient and parameter inference from diagnostic measurements through the differentiable inverse procedure, replacing iterative manual SOLPS-ITER fitting [43]; and (i) real-time or near-real-time edge-plasma state estimation for control applications [6, 35], conditioned on a limited set of diagnostic measurements. Comparable digital-twin paradigms have been demonstrated for particle accelerator control [31, 38]. 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 [2], 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 multi-dimensional Latin hypercube design [22] by selecting, for each new sample, the nearest previously converged neighbor in (logN)( N) time [4, 13]. The database grows as a branching process: starting from a single seed, each newly converged case extends the frontier of reliable initial conditions, progressively reducing the parameter-space distance to subsequent targets. Beyond this branching construction, the warm-start strategy improves robustness: in a prior cold-started ensemble, roughly 30%30\% of runs initiated from generic initial conditions diverged outright, and a further fraction never reached steady state [5], whereas the branching approach achieved a completion rate above 95%95\% across the full five-dimensional design, yielding the 762762-run database used here. This infrastructure is essential for any surrogate relying on a large, uniformly sampled training set; without it, gaps in the database would directly translate into blind spots in the learned model. Several limitations remain, and we state them plainly. First, the current database covers only a single DIII-D lower-single-null configuration with deuterium-only fueling and spatially uniform transport coefficients; extending it to other machines, geometries, plasma mixtures, and spatially varying transport is necessary for validation and, perhaps, for transfer learning. Second, the surrogate has not yet been validated against experimental DIII-D measurements, a critical step before deployment. Third, physics constraints beyond admissibility, such as global radiated power and ZeffZ_eff consistency, and the conservation of particles, momentum, and energy, are not yet enforced; imposing them is a natural way to sharpen the physical fidelity of the surrogate. Finally, although the current inverse procedure is warm-started by the learned pseudo-inverse GψG_ψ (Sec. 3.4), the refinement still requires gradient-based optimization through the frozen forward model; end-to-end joint training with cycle consistency could further improve convergence speed and inverse accuracy. 7 Conclusions A cycle-consistent neural surrogate framework is introduced for modeling two-dimensional SOLPS-ITER edge-plasma fields. The framework incorporates a k-d tree warm-start strategy that selects nearest-neighbor restarts from a database of previously converged SOLPS-ITER simulations, thereby reducing both simulation failure rates and wall-clock time: across the five-dimensional parameter scan, the warm start reached a completion rate above 95%95\%, yielding 762762 converged runs, whereas roughly 30%30\% of comparable cold-started runs have been reported to fail outright, with more never reaching steady state [5]. A parameter-conditioned U-Net maps scalar control parameters directly to the two-dimensional plasma state (TeT_e, TiT_i, nen_e, uau_a), and achieves Pearson correlation coefficients exceeding 0.950.95 for all four plasma state variables. Cycle-consistency regularization is shown to yield a self-consistent forward/inverse pair at negligible cost to forward accuracy, supplying a self-supervised consistency signal that requires no ground-truth labels. Finally, an inverse inference procedure coupled with cycle-consistency regularization enables parameter recovery from target fields. It ensures that reconstructed outputs lie on the learned physics manifold, providing a self-supervised reliability metric that requires no ground-truth labels at inference time. For held-out cases, the inverse recovers all five control parameters with Pearson correlations of 0.970.97 or higher, including the cross-field transport coefficients D⟂D_ and χi _i (both r≈0.99r≈ 0.99); since these are typically set by hand-tuned SOLPS-ITER fits to diagnostics, the differentiable inverse offers an automated route toward transport-coefficient inference. A companion query-by-committee profile ensemble supplies predictive uncertainty estimates whose committee disagreement concentrates on the physically hardest regime, the near-strike-point inner divertor closest to detachment (flagging roughly 13%13\% of inner-target points against 3%3\% at the outer target), and thereby nominates the most informative new SOLPS-ITER runs for active-learning acquisition. Looking ahead, several directions are being pursued. First, validation against experimental DIII-D measurements from discharge 174310 is underway, which will test the surrogate’s ability to generalize beyond the simulation database. Second, extending the inverse procedure to operate on sparse one-dimensional diagnostic signals (e.g., Thomson scattering profiles) would directly replace the manual transport-coefficient fitting workflow of interpretive boundary modeling [43]. Third, incorporating additional physics constraints (global radiated power PradP_rad and effective charge ZeffZ_eff) and extending to spatially varying transport coefficients will improve the surrogate’s physical fidelity and applicability. Finally, extending the approach to other magnetic configurations and machines (e.g., ITER, SPARC, MAST-U) will test its generality. The companion neutral-source model, which maps the plasma state predicted here to EIRENE source terms for inline replacement of the kinetic neutral step, is developed separately. CRediT authorship contribution statement A. Diaw: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data curation, Visualization, Writing – original draft. S. De Pascuale: Conceptualization, Methodology. J.-S. Park: Conceptualization, Methodology. I. Paradela Perez: Conceptualization, Methodology. J.D. Lore: Software, Resources, Methodology, Funding acquisition. S. Dasbach: Conceptualization, Methodology. Code and Data availability The code is available at https://github.com/abdoudiaw/solpex. The sampled data is archived on figshare at https://doi.org/10.6084/m9.figshare.32048490. Acknowledgments This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Fusion Energy Sciences. This research used resources of the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory, which is supported by DOE under Contract DE-AC05-00OR22725. This research used resources of the Oak Ridge National Laboratory Research Cloud. The results are obtained with the help of the EIRENE package (see w.eirene.de) including the related code, data and tools [29]. DIFFER is part of the institutes organisation of NWO. Declaration of generative AI and AI-assisted technologies in the writing process During the preparation of this work, the authors used Claude (Anthropic) to improve the manuscript’s readability and language. After using this tool, the authors reviewed and edited the content as needed and took full responsibility for the published article. Appendix A k-d Tree Warm-Start Construction. We build a balanced k-d tree over all LHS design points i\x_i\ by recursively splitting the data along coordinate directions. This structure supports nearest-neighbor queries in (logN)O( N) time on average for d≪Nd N, where d=5d=5 is the number of input parameters and N is the number of design points. Given a target input ⋆t_ , we define a weighted Euclidean distance dw2(i,⋆)=∑j=15wj(xi,j−t⋆,j)2,d_w^2(x_i,t_ )= _j=1^5w_j (x_i,j-t_ ,j )^2, (11) where wjw_j is a user-defined weight for the j-th parameter (here wj=1w_j=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. 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 [2] generalized to a multi-dimensional Latin hypercube design [22]. Figure 7: The k-d tree warm-start strategy applied during database generation, projected onto the (Ptot,ΓD2)(P_tot,\, _D_2) 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. Appendix B Database Convergence Detection 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 block-averaging techniques from equilibrium sampling in molecular dynamics [1, 12]. The same machinery underpins Green-Kubo transport-coefficient calculations in molecular dynamics simulations [10, 40]; here, we repurpose it to estimate the correlation time of SOLPS-ITER diagnostics and construct a statistically grounded convergence criterion [11, 37]. We monitor four scalar diagnostics extracted at each SOLPS-ITER timestep: the electron and ion temperatures at the outboard-midplane separatrix (Te,sepOMPT_e,sep^OMP, Ti,sepOMPT_i,sep^OMP), the electron density at the same location (ne,sepOMPn_e,sep^OMP), and the total number of particles (⟨tmne⟩ tmne ). Each diagnostic is treated as a noisy relaxation process whose fluctuations encode the residual coupling between the plasma and neutral subsystems. Given a time series qkk=0N−1\q_k\_k=0^N-1 with sample mean q¯ q, we compute the normalized autocorrelation function C(ℓ)=∑k=0N−ℓ−1(qk−q¯)(qk+ℓ−q¯)∑k=0N−1(qk−q¯)2,ℓ=0,1,2,…,C( )= _k=0^N- -1(q_k- q)\,(q_k+ - q) _k=0^N-1(q_k- q)^2, =0,1,2,…, (12) so that C(0)=1C(0)=1 and |C(ℓ)|≤1|C( )|≤ 1 for ℓ>0 >0. The integral correlation time is then estimated as τint≈Δt(12+∑ℓ=1ℓmaxC(ℓ)), _int≈ t ( 12+ _ =1 _ C( ) ), (13) where the sum is truncated at the first non-positive value of C(ℓ)C( ) or at a prescribed maximum lag [37]. This quantity plays the same role as the correlation time in MD block averaging [11]: it sets the minimum window length over which successive block means can be treated as approximately independent. Using τint _int, we define a tail region consisting of the last Ttail=NtailτintT_tail=N_tail\, _int timesteps of the series and split it into two non-overlapping windows of equal length L=NwinτintL=N_win\, _int, denoted A and B. The drift between the two window means, Δq¯=q¯B−q¯A q= q_B- q_A, is compared against its expected statistical fluctuation, σΔq¯≈22τintLσq, _ q≈ 2\, 2\, _intL\; _q, (14) where σq _q is the standard deviation within each window. A quantity is declared to be in a steady state when |Δq¯||q¯B|<εreland|Δq¯|σΔq¯<zmax, | q|| q_B|< _rel | q| _ q<z_ , (15) with εrel=0.01 _rel=0.01 and zmax=2z_ =2. The first condition ensures that the mean has not drifted by more than 1%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. (15). Figure 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. Figure 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 _int (dashed line). 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 held-out test case (Fig. 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. Figure 9: Predicted electron profiles as a function of normalized poloidal flux ψN _N for a representative held-out test case. Solid lines: SOLPS-ITER ground truth. Dashed lines (“N forward”): the forward surrogate FθF_θ evaluated at the true control parameters, isolating forward-model accuracy. Dotted lines (“N inverse cycle”): the cycle reconstruction, in which the control parameters are first recovered from the SOLPS-ITER fields by the inverse procedure G⋆G and then propagated back through the forward model. Panels show (a) upstream (outboard-midplane) nen_e, (b) upstream TeT_e, (c) outer-target TeT_e, and (d) inner-target TeT_e. 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. References [1] M. P. Allen and D. J. Tildesley (2017) Computer simulation of liquids. 2nd edition, Oxford University Press. External Links: Document Cited by: Appendix B. [2] E. L. Allgower and K. Georg (2003) Introduction to numerical continuation methods. Classics in Applied Mathematics, Vol. 45, SIAM. External Links: Document Cited by: Appendix A, §6. [3] R. Anirudh, J. J. Thiagarajan, P. Bremer, and B. K. Spears (2020) Improved surrogates in inertial confinement fusion with manifold and cycle consistencies. Proceedings of the National Academy of Sciences 117 (18), p. 9741–9746. External Links: Document Cited by: §3.3, §3.4, §3.4, §4.3, §6. [4] J. Bentley (1975) Multidimensional binary search trees used for associative searching. Communications of the ACM 18 (9), p. 509–517. External Links: Document Cited by: §2, §6. [5] S. Dasbach and S. Wiesen (2023) Towards fast surrogate models for interpolation of tokamak edge plasmas. Nuclear Materials and Energy 34 (101396), p. 1–6. External Links: Document Cited by: §2, §2, §5, §6, §7. [6] J. Degrave, F. Felici, J. Buchli, M. Neunert, B. Tracey, F. Carpanese, T. Ewalds, R. Hafner, A. Abdolmaleki, D. de Las Casas, et al. (2022) Magnetic control of tokamak plasmas through deep reinforcement learning. Nature 602, p. 414–419. External Links: Document Cited by: §1, §4.2, §6. [7] 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, and T. C. Germann (2020-08) Multiscale simulation of plasma flows using active learning. Phys. Rev. E 102 (2), p. 023310. External Links: Document Cited by: §1, §5, §6. [8] A. Diaw, M. McKerns, I. Sagert, L. G. Stanton, and M. S. Murillo (2022-07) Efficient Learning of Accurate Surrogates for Simulations of Complex Systems. arXiv e-prints, p. arXiv:2207.12855. External Links: Document, 2207.12855 Cited by: §1, §5. [9] A. Diaw, M. McKerns, I. Sagert, L. G. Stanton, and M. S. Murillo (2024) Efficient learning of accurate surrogates for simulations of complex systems. Nature Machine Intelligence 6, p. 568–577. External Links: Document Cited by: §1. [10] A. Diaw and M. S. Murillo (2015-07) Generalized hydrodynamics model for strongly coupled plasmas. Phys. Rev. E 92 (1), p. 013107. External Links: Document Cited by: Appendix B. [11] H. Flyvbjerg and H. G. Petersen (1989) Error estimates on averages of correlated data. The Journal of Chemical Physics 91, p. 461–466. External Links: Document Cited by: Appendix B, Appendix B. [12] D. Frenkel and B. Smit (2002) Understanding molecular simulation: from algorithms to applications. 2nd edition, Academic Press. Cited by: Appendix B. [13] J. H. Friedman, J. L. Bentley, and R. A. Finkel (1977) An algorithm for finding best matches in logarithmic expected time. ACM Transactions on Mathematical Software 3 (3), p. 209–226. External Links: Document Cited by: §6. [14] J. Kates-Harbeck, A. Svyatkovskiy, and W. Tang (2019) Machine learning for disruption warnings on Alcator C-Mod, DIII-D, and EAST. Nature 568 (7753), p. 526–531. External Links: Document Cited by: §1. [15] D. P. Kingma and J. Ba (2017) Adam: a method for stochastic optimization. External Links: 1412.6980, Link Cited by: §3.2, §3.5, §5. [16] D. P. Kingma and M. Welling (2019-11) An introduction to variational autoencoders. Foundations and Trends® in Machine Learning 12 (4), p. 307–392. External Links: ISSN 1935-8245, Link, Document Cited by: §1. [17] A. Krizhevsky, I. Sutskever, and G. E. Hinton (2012) Imagenet classification with deep convolutional neural networks. Advances in neural information processing systems 25. Cited by: §3.1. [18] B. Kustowski, J. A. Gaffney, B. K. Spears, G. J. Anderson, J. J. Thiagarajan, and R. Anirudh (2020-09) 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), p. 3275–3275. External Links: Document Cited by: §3.3. [19] A. Lasa, J. -S. Park, J. Lore, S. Blondel, D. E. Bernholdt, J. M. Canik, M. Cianciosa, J. Coburn, and D. Curreli (2024) Exploreing the efect of elm and code-coupling frequencies on plasma and material modeling of dynamic recycling in divertors. Nuclear Fusion 64 (7), p. 1–13. External Links: Document Cited by: §2. [20] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel (1989) Backpropagation applied to handwritten zip code recognition. Neural computation 1 (4), p. 541–551. Cited by: §3.1. [21] J. D. Lore, S. De Pascuale, P. Laiu, B. Russo, J. -S. Park, J. M. Park, S. L. Brunton, J. N. Kutz, and A. A. Kaptanoglu (2023) Time-dependent solps-iter simulations of the the tokamak plasma boundary for model predictive control using sindy. Nuclear Fusion 63 (046015), p. 1–12. External Links: Document Cited by: §2. [22] M. D. McKay, R. J. Beckman, and W. J. Conover (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), p. 239–245. External Links: Document Cited by: Appendix A, §6. [23] M.M. Noack, K.G. Yager, M. Fukuto, G.S. Doerk, L. Ruipeng, and J.A. Sethian (2019) A kriging-based approach to autonomous experimentation with applications to x-ray scattering. Scientific Reports 9 (11809). Cited by: §1, §6. [24] J. -S. Park, J. D. Lore, M. Reinke, A. Q. Kuang, S. De Pascuale, and A. Creely (2024) Full time-dependent solps-iter simulation of the sparc tokamak: actuator design for particle and divertor condition control. Nuclear Fusion 64 (4), p. 1–14. External Links: Document Cited by: §2. [25] 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, and S. Chintala (2019) PyTorch: an imperative style, high-performance deep learning library. External Links: 1912.01703, Link Cited by: §3.2. [26] E. Perez, F. Strub, H. de Vries, V. Dumoulin, and A. Courville (2018) FiLM: visual reasoning with a general conditioning layer. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 32, p. 3942–3948. External Links: Document, Link Cited by: §3.1. [27] R. A. Pitts, X. Bonnin, F. Escourbiac, T. Hirai, V. Komarov, A. S. Kukushkin, A. Loarte, A. Martin, M. Merola, and R. Mitteau (2019) Divertor and scrape-off layer plasma physics for reactor-grade fusion energy devices. Nuclear Materials and Energy 20, p. 100696. External Links: Document Cited by: §1. [28] D. Reiter, M. Baelmans, and P. Börner (2005) EIRENE – a Monte Carlo linear transport solver. Technical report Forschungszentrum Jülich. Note: http://w.eirene.de Cited by: §1. [29] D. Reiter, M. Baelmans, and P. Borner (2005) The eirene and b2-eirene codes. Fusion Science and Technology 47 (2), p. 172–186. External Links: Document Cited by: §2, Acknowledgments. [30] O. Ronneberger, P. Fischer, and T. Brox (2015) U-net: convolutional networks for biomedical image segmentation. External Links: 1505.04597, Link Cited by: §1, §3.1. [31] R. Roussel, A. Edelen, C. Mayes, D. Ratner, J. P. Gonzalez-Aguilera, S. Kim, E. Wisniewski, and J. Power (2023) Phase space reconstruction from accelerator beam measurements using neural networks and differentiable simulations. Physical Review Letters 130, p. 145001. External Links: Document Cited by: §6. [32] 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, and P. Strand (2024) Data-driven models in fusion exhaust: ai methods and perspectives. Nuclear Fusion 64 (086046), p. 1–10. External Links: Document Cited by: §2. [33] A. Scheinker and S. Gessner (2015) Adaptive method for electron bunch profile prediction. Physical Review Accelerators and Beams 18 (102801). Cited by: §1. [34] R. Schneider, X. Bonnin, K. Borrass, D. P. Coster, H. Kastelewicz, D. Reiter, V. A. Rozhansky, and B. J. Braams (2006) Plasma edge physics with B2-Eirene. Contributions to Plasma Physics 46 (1-2), p. 3–191. External Links: Document Cited by: §1. [35] J. Seo, S. Kim, A. Jalalvand, R. Conlin, A. Rothstein, J. Abbate, K. Erickson, J. Wai, R. Shousha, and E. Kolemen (2024) Avoiding fusion plasma tearing instability with deep reinforcement learning. Nature 626, p. 746–751. External Links: Document Cited by: §1, §4.2, §6. [36] H. S. Seung, M. Opper, and H. Sompolinsky (1992) Query by committee. In Proceedings of the Fifth Annual Workshop on Computational Learning Theory, COLT ’92, New York, NY, USA, p. 287–294. External Links: ISBN 089791497X, Link, Document Cited by: §5. [37] A. D. Sokal (1997) Monte Carlo methods in statistical mechanics: foundations and new algorithms. Lecture Notes in Physics 490, p. 131–192. Cited by: Appendix B, Appendix B. [38] 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, J. M. Duarte, Y. Huang, and R. Keller (2021) Real-time artificial intelligence for accelerator control: a study at the Fermilab Booster. Physical Review Accelerators and Beams 24, p. 104601. External Links: Document Cited by: §4.2, §6. [39] W. Tang, E. Feibush, G. Dong, and N. Borthwick (2024) AI-machine learning-enabled tokamak digital twin. arXiv preprint. External Links: 2409.03112 Cited by: §1, §4.2. [40] C. Ticknor, J. D. Kress, L. A. Collins, J. Clérouin, P. Arnault, and A. Decoster (2016-06) Transport properties of an asymmetric mixture in the dense plasma regime. Phys. Rev. E 93 (6), p. 063208. External Links: Document Cited by: Appendix B. [41] 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, and S. Voskoboynikov (2015) The new solps-iter code package. Journal of Nuclear Materials 463 (), p. 480–484. External Links: Document Cited by: §2. [42] P.B. Wigley, P.J. Everitt, A. van den Hengel, J.W. Bastian, M.A. Sooriyabandara, G.D. McDonald, K.S. Hardman, C.D. Quinlivan, P. Manju, C.C.N. Kuhn, I.R. Peterson, A.N. Luiten, J.J. Hope, N.P. Robins, and M.R. Hush (2016) Fast machine-learning online optimization of ultra-cold-atom experiments. Scientific Reports 6 (25890). Cited by: §1. [43] R. S. Wilcox, M. W. Shafer, J. D. Lore, J. M. Canik, S. R. Haskey, C. J. Lasnier, A. L. Moser, T. H. Osborne, and H. Q. Wang (2026) Challenges and approaches to interpretive modeling of boundary plasma and neutral transport in a closed, pumped divertor. Nuclear Fusion 66, p. 016023. External Links: Document Cited by: §6, §7. [44] K. Willcox and B. Segundo (2024) The role of computational science in digital twins. Nature Computational Science 4, p. 147–149. External Links: Document Cited by: §1, §4.2. [45] B. Zhu, M. Zhao, H. Bhatia, X. Xu, P. Bremer, W. Meyer, N. Li, and T. Rognlien (2022) Data-driven model for divertor plasma detachment prediction. Journal of Plasma Physics 88 (895880504), p. 1–23. External Links: Document Cited by: §2, §5. [46] B. Zhu, M. Zhau, X. Xu, A. Gupta, K. -B. Kwon, X. Ma, and D. Eldon (2025) Latent space mapping: revolutionizing predictive models for divertor plasma detachment control. Physics of Plasmas 32 (062508), p. 1–19. External Links: Document Cited by: §2.