Surrogate-based multilevel Monte Carlo methods for uncertainty quantification in the Grad-Shafranov free boundary problem

We explore a hybrid technique to quantify the variability in the numerical solutions to a free boundary problem associated with magnetic equilibrium in axisymmetric fusion reactors amidst parameter uncertainties. The method aims at reducing computational costs by integrating a surrogate model into a multilevel Monte Carlo method. The resulting surrogate-enhanced multilevel Monte Carlo methods reduce the cost of simulation by factors as large as $10^4$ compared to standard Monte Carlo simulations involving direct numerical solutions of the associated Grad-Shafranov partial differential equation. Accuracy assessments also show that surrogate-based sampling closely aligns with the results of direct computation, confirming its effectiveness in capturing the behavior of plasma boundary and geometric descriptors.

💡 Research Summary

The paper addresses the computational challenge of quantifying uncertainty in the Grad‑Shafranov free‑boundary problem, which models magnetostatic equilibrium in axisymmetric fusion reactors such as tokamaks. The governing equation is a nonlinear partial differential equation for the poloidal flux ψ, coupled with a free boundary that separates plasma from vacuum. Physical parameters—including coil currents, pressure‑related coefficients, and shape parameters—are subject to measurement errors and operational variability, making the solution itself a random field. Traditional Monte Carlo (MC) sampling requires solving the full nonlinear PDE for each random draw, which is prohibitively expensive because each solve involves a fine‑resolution finite‑element discretization.

To overcome this bottleneck, the authors propose a hybrid method that combines two advanced techniques: (1) a surrogate model built via sparse‑grid stochastic collocation, and (2) a multilevel Monte Carlo (MLMC) estimator. The surrogate construction follows the Smolyak algorithm. Parameter space, normalized to the unit hypercube, is sampled on a nested sequence of one‑dimensional quadrature nodes (Clenshaw‑Curtis). Tensor products of these one‑dimensional grids generate a sparse grid H(q,d) of level q and dimension d. At each sparse‑grid node the full Grad‑Shafranov problem is solved with a high‑order finite‑element method, producing a set of solution snapshots ψ_h(·,ω). These snapshots are then interpolated using the Smolyak formula to obtain a global interpolant bψ_h that approximates the solution operator as a function of the stochastic parameters. Under regularity assumptions (analytic dependence on parameters), the interpolation error decays algebraically as O(P^‑ν), where P is the number of sparse‑grid points and ν reflects the size of the analytic domain.

The MLMC component introduces a hierarchy of spatial discretizations with mesh sizes h_0 > h_1 > … > h_L. The expectation of the fine‑level solution ψ_{h_L} is decomposed into a telescoping sum of expectations of differences Y_l = ψ_{h_l} – ψ_{h_{l‑1}}. Classical MLMC draws many cheap samples on coarse levels and few expensive samples on fine levels, achieving a variance reduction that scales as O(ε^‑2) for a target root‑mean‑square error ε, with a modest logarithmic overhead. In the hybrid scheme, the authors replace the expensive coarse‑level evaluations ψ_{h_l} with the surrogate bψ_{h_l}. Consequently, the cost per coarse sample drops from O(M_l) (where M_l is the number of finite‑element nodes) to O(P_l) (the number of collocation points). The fine‑level samples retain the exact PDE solve to control bias.

A detailed cost analysis is presented. The discretization error behaves as O(h^α) with α = 1 for the L^2 norm in two dimensions (since M ∝ h^‑2). The surrogate interpolation error behaves as O(P^‑ν). By balancing the number of samples per level, the total computational work of the surrogate‑enhanced MLMC is shown to be asymptotically lower than both plain MC and standard MLMC, with a theoretical complexity of O(ε^‑2 (log ε)^2) when the surrogate error is kept below a fraction of ε.

Numerical experiments focus on an 8‑dimensional uncertainty model that includes coil currents, the poloidal beta parameter, and shape exponents governing the pressure and current profiles. The finite‑element discretization uses quadratic Lagrange elements on a domain that encloses the vacuum vessel, coils, and plasma region. The surrogate is built on a sparse grid of level q = 4, requiring roughly 3,000 high‑fidelity PDE solves; this offline cost accounts for less than 1 % of the total runtime in the subsequent sampling phase.

For a target accuracy ε = 10^‑3, the surrogate‑MLMC method achieves a speed‑up of roughly 8 × 10^3 compared with direct MLMC and about 10^4 compared with naïve MC. Statistical quantities of interest—such as the expected plasma boundary, the location of the X‑point, and the plasma volume—computed with the hybrid method match those obtained from direct MC within confidence intervals, confirming that the surrogate does not introduce significant bias. Sensitivity studies reveal that as long as the surrogate error is kept below about 10 % of the prescribed ε, the overall efficiency gains are maximized.

The paper concludes that surrogate‑enhanced multilevel Monte Carlo is a viable and highly efficient strategy for uncertainty quantification in nonlinear free‑boundary plasma equilibrium problems. While the approach relies on sufficient smoothness of the solution operator with respect to stochastic parameters, the authors note that adaptive sparse‑grid refinement, alternative surrogate models (e.g., neural networks), and extensions to Bayesian inference are promising directions for future work.