A general framework for quantifying uncertainty at scale

The paper addresses the challenge of performing uncertainty quantification (UQ) and global sensitivity analysis (SA) for complex, computationally expensive, physics-based simulations. Realistic models often require supercomputers, making brute-force sampling infeasible and limiting predictive capability. The authors propose a framework based on a sensitivity-driven, dimension-adaptive sparse grid interpolation that adaptively explores important input directions and anisotropic interactions, mitigating the curse of dimensionality. The study focuses on a key problem in fusion research—quantifying and predicting turbulent transport in magnetic confinement devices—where simulations are extremely costly and inputs are uncertain. The goal is to enable accurate UQ and SA at scale while producing an efficient surrogate model.

The authors situate their work within decades of advances in model-based and data-centric computing and the established need for UQ and SA in predictive simulations. They cite foundational works on stochastic finite elements, UQ, and global SA, and note the prohibitive costs of ensemble evaluations in large-scale simulations. Reduced-order models can help but often at accuracy and training-cost trade-offs. Sparse grid and adaptive collocation approaches have been explored previously, including applications to linear gyrokinetic simulations, but most standard methods remain infeasible for fully nonlinear, large-scale problems. The proposed method builds on prior sensitivity-driven adaptive sparse approximations and dimension-adaptive sparse grids, aiming to exploit lower intrinsic dimensionality and anisotropic parameter interactions in a general, non-intrusive way applicable across disciplines (e.g., CFD, combustion, climate, materials, geophysics, medicine, epidemiology).

The framework employs a sensitivity-driven, dimension-adaptive sparse grid interpolation to perform UQ and SA efficiently. Uncertain inputs are modeled as independent random variables with a product density π over a product domain X = X1 × ... × Xd. The high-fidelity model f_hi maps θ ∈ R^d to a scalar output of interest. Sparse grid approximations are formed as linear combinations of d-variate products of univariate Lagrange interpolation operators on weighted (L)-Leja points, enabling hierarchical refinement and reuse of evaluations. A multi-index set L ⊂ N^d (admissible/downward closed) defines tensorized subspaces; it is split into old (O) and active (A) sets for sequential, adaptive refinement. In each step, for each candidate subspace, the method computes unnormalized sensitivity indicators by decomposing variance contributions into first-order and interaction components (equivalent to unnormalized Sobol’-type indices via a spectral-projection representation). User-defined tolerances for each index yield an integer sensitivity score per subspace, incremented when indicators exceed tolerances. The subspace with the largest score is refined (ties broken by largest sum of indicators), preferentially exploring important directions and interactions. For moderate d, all 2^d − 1 indices can be used; for larger d, only up to pairwise or triple interactions are considered. The interpolation is implemented stably via barycentric formula using (L)-Leja sequences (hierarchical, interpolatory, accurate, applicable for arbitrary densities); for uniform densities, the first point is the interval midpoint. The interpolation is equivalently expressed in an orthonormal polynomial basis (e.g., Legendre for uniform inputs), allowing straightforward computation of hierarchical spectral surpluses and sensitivity measures that drive adaptivity. The algorithm is non-intrusive, requiring only parameter specification and extraction of outputs from the high-fidelity solver; it runs from laptops to supercomputers. Outputs at termination include mean/variance, sensitivity indices (first-order and selected interactions), and an interpolation-based surrogate of the parameter-to-output map. Application setup: fully nonlinear gyrokinetic turbulence simulations (GENE code) of ETG-driven transport in the near-edge (ρ = 0.95) pedestal region of a DIII-D tokamak discharge. Grid: 256 × 24 × 168 × 32 × 8 ≈ 2.642×10^8 DoF. Computing: Frontera supercomputer (16 nodes, 896 cores), average run time >8000 core-hours (min ~4000, max >14,000). Eight uncertain inputs modeled as uniform RVs with representative experimental bounds: Te and ne varied ±10%; ω_Te, ω_ne, Ti/Te (τ), Zeff, q, and magnetic shear ŝ varied ±20%. Output of interest Q_high is the time-averaged electron heat flux across the flux surface (MW) over a quasi-steady interval, converted to SI units for consistency. Tolerances: a vector of length 2^8−1 = 63 with all entries 10^−4. The method adaptively selects sparse-grid points; 57 high-fidelity simulations were required to meet tolerances.

• The adaptive method completed accurate UQ and SA with only 57 high-fidelity simulations for an 8D problem with ~2.64×10^8 DoF per run. - Mean and variance of the output: E[Q_high] ≈ 0.7530 MW; Var[Q_high] = 0.2571 MW^2. - Sensitivity analysis identified four important parameters: ω_Te (inverse electron temperature scale length), ω_ne (inverse electron density scale length), Te, and the temperature ratio τ = Ti/Te; ω_Te and ω_ne are the most influential. - Interactions are anisotropic; the strongest interaction is between ω_Te and ω_ne; the second strongest involves ω_Te with another key parameter (reported as τ). - Adaptive sampling concentrated points in planes involving ω_Te and ω_ne, with few points for unimportant inputs (e.g., Zeff, magnetic shear), confirming structure exploitation. - Computational savings: a naive 3-point-per-dimension full tensor grid would need 3^8 = 6561 simulations (~53 million core-hours). The adaptive approach used ~460,000 core-hours, a factor ≈115 reduction. - The interpolation surrogate closely matches high-fidelity outputs on 32 independent test samples spanning ~0.1–2.6 MW. Mean squared error MSE = 3.0707×10^−4. - Average surrogate evaluation cost c_SG = 9.4046×10^−3 s, about nine orders of magnitude faster than a high-fidelity run. - Surrogate predictions stabilize by ~36–57 points with negligible variation, indicating small prediction uncertainty for tested cases.

The study demonstrates that sensitivity-driven dimension-adaptive sparse grids can enable robust UQ and SA for complex, fully nonlinear, and computationally intensive simulations that are otherwise infeasible with standard methods. By exploiting lower intrinsic dimensionality and anisotropic interactions, the method concentrates computational effort on influential inputs and their key interactions, thereby overcoming the curse of dimensionality in practice. In the fusion application, it systematically assessed the impact of experimentally uncertain inputs on turbulent electron heat flux and identified the most influential parameters and interactions. The accurate surrogate provides rapid evaluations, enabling exploration of parameter regimes beyond current experimental accessibility, and supporting downstream tasks such as design optimization and multi-fidelity workflows. These results address the central research question by providing accurate statistics and sensitivity information at a fraction of the cost, without intruding on the high-fidelity solver.

The paper introduces and validates a general, non-intrusive framework for large-scale UQ and SA based on sensitivity-driven, dimension-adaptive sparse grid interpolation. In a challenging fusion turbulence case with eight uncertain inputs and costly simulations, the framework achieved accurate statistics and sensitivity indices using only 57 runs, and produced an accurate surrogate model nine orders of magnitude cheaper to evaluate. The approach is broadly applicable to computational science and engineering problems and integrates seamlessly into existing pipelines. Future work includes extending to basis functions with local support (e.g., wavelets) to handle discontinuities or sharp parameter-space gradients, developing stronger theoretical convergence guarantees for general problems, and incorporating data-driven methods to quantify prediction uncertainty while maintaining generality. The surrogate can further support optimization and multi-fidelity studies in fusion and beyond.

• Global interpolation polynomials are not well-suited for problems with discontinuities or sharp gradients in the parameter-to-output map; local-support bases (e.g., wavelets) may be required. - Accuracy of dimension-adaptive algorithms is often assessed empirically; rigorous convergence guarantees are limited to specific settings. - The method assumes independent inputs (product structure of domain and density); while relaxable via transformations (e.g., transport maps), dependence requires additional preprocessing. - As with any surrogate-based UQ/SA, validation is limited by the feasible number of high-fidelity test samples in extremely expensive applications.