Uncertainty Quantification for Multiscale Fusion Plasma Simulations with VECMA Toolkit

Abstract

Within VECMAtk platform we perform Uncertainty Quantification (UQ) for multiscale fusion plasmas simulations. The goal of VECMAtk is to enable modular and automated tools for a wide range of applications to archive robust and actionable results. Our aim in the current paper is to incorporate suitable features to build UQ workflow over the existing fusion codes and to tackle simulations on high performance parallel computers.

1 Introduction

Thermonuclear fusion can potentially provide a carbon free and safe solution to the provision of base load electricity. Heat and particle transport will play the key roles in determining the size and efficiency of future fusion power reactor, so understanding their mechanisms is an important milestone towards the realization of fusion-based electricity. Our present understanding is that turbulence at small space and time scales is responsible for much of this transport, but the profiles of temperature and density evolve over much larger space and time scales. To study these phenomena, multiscale and multiphysics applications have been developed throughout the years. Notable efforts in Europe lean toward implementing such applications by following a workflow approach, in which several single-scale (or single-physics) codes are coupled together [1, 2].

When one element in such a workflow is a turbulence code, whose outputs are inherently noisy due to the chaotic nature of the turbulence (intrinsic uncertainties), or if input data such as external sources contain uncertainties (extrinsic), these need to be propagated through different models of the workflow. The goal, therefore, is to produce temperature and density profiles, along with their confidence intervals. This information will then allow for an improved validation of the simulation results against the experimental measurements that come with error bars.

Uncertainty quantification in the existing fusion workflow is implemented with the VECMA open source toolkit¹, which provides tools to facilitate the verification, validation and uncertainty quantification (VVUQ) process in multiscale, multiphysics applications [3]. So far we have integrated two elements from this toolkit: the EasyVVUQ library [4] to incorporate uncertainty quantification (UQ) and sensitivity analysis (SA) in a non-intrusive manner, and the QCG Pilot-Job middleware [5] to execute efficiently the large number of tasks or samples prepared by EasyVVUQ.

In this work we focus on the different sampling methods used to quantify the uncertainties, their results and performance for two different black box use cases: on the entire workflow for extrinsic uncertainties coming from the heating source terms and boundary conditions for temperature profiles (on electron and ions), and on the turbulence model for uncertain ions and electron temperatures values and gradients. This paper is organized as follows. After a brief overview on the background theory for uncertainty quantification, we present the various sampling methods used in our application that are also available in the EasyVVUQ. We then describe our two use cases together with UQ and SA results for the method that demonstrates the best performance. Finally, we conclude with future works related to more intrusive approaches.

2 Theory Overview

We consider a numerical model $\mathcal {M}$ that has a number d of uncertain input parameters $Q_1, Q_2, ..., Q_d$, and gives the output Y:

$$\begin{aligned} Y = \mathcal {M}(x , Q). \end{aligned}$$

(1)

Here, Q is the vector $[Q_1, Q_2, ..., Q_d]$ and x is a space variable. The model can also depends on other variables (e.g. time) and on additional but invariant input parameters. In this study we focus on the uncertain parameters $({Q_i})_{i=1}^d$ and how the uncertainty is propagated through the model $\mathcal {M}$. That way we can quantify the uncertainty in the output Y and identify the relative contribution of each $Q_i$ into this uncertain result.

First, we assume that the uncertain parameters are statistically independent components, and each one has an univariate probability density function $f_{{Q}_i}$, hence the vector Q can be described by the joint probability density function:

$$\begin{aligned} f_{Q} = \prod _{i=1}^d f_{{Q}_i}. \end{aligned}$$

(2)

Thereafter, we seek relevant properties of the output Y that characterize the uncertainty propagation and where possible, we build an approximation of the unknown output distribution $f_Y$.

2.1 Statistical Moments

The principal aim of UQ is to describe the output distribution $f_Y$ through useful statistics [6, 7]. Several metrics are based on statistical moments, and among these the most commonly used are the mean $\mathbb {E}$ and the variance $\mathbb {V}$:

$$\begin{aligned} \mathbb {E}[Y] = \int _{\varOmega _Y} yf_Y(y)dy, \end{aligned}$$

(3)

$$\begin{aligned} \mathbb {V}[Y] = \int _{\varOmega _Y} (y-\mathbb {E}[Y])^2f_Y(y)dy. \end{aligned}$$

(4)

Here, ${\varOmega _Y}$ is the output space.

The mean indicates the expected value of the model output Y and the variance measures how much the output varies around the mean. In general, instead of using the variance, we take its square root to obtain the standard deviation $\sigma$:

$$\begin{aligned} \sigma [Y]= \sqrt{\mathbb {V}[Y]}. \end{aligned}$$

(5)

The advantage of using $\sigma$ instead of the variance is that, like the mean, $\sigma$ has the same units as Y and that makes it easier to interpret.

Another useful measure for UQ is percentiles which can be used to define confidence intervals [6].

2.2 Sensitivity Analysis

The principal aim of SA is to evaluate the contribution of different inputs parameters $(Q_i)_{i=1}^d$ to the uncertainty in output Y. Among the various sensitivity metrics that exist, we choose in this study the variance-based sensitivity analysis by computing the Sobol indices [8, 9]. It is the most commonly used metric as it allows to measure sensitivity across the whole input space and can cope with nonlinear outputs.

There are several orders of Sobol indices. Here, we select two of the most-used indices in UQ: the first-order Sobol index $S_i$ and the total order index $S_{Ti}$. They measure the direct effect of the parameter $Q_i$ on the variance of Y and the interactions between $Q_i$ and other parameters, respectively:

$$\begin{aligned} S_i = \frac{\mathbb {V}[\mathbb {E}[Y|Q_i]]}{\mathbb {V}[Y]}, \end{aligned}$$

(6)

$$\begin{aligned} S_{Ti} = 1- \frac{\mathbb {V}[\mathbb {E}[Y|Q_{-i}]]}{\mathbb {V}[Y]}. \end{aligned}$$

(7)

We denote by $\mathbb {E}[Y|Q_i]$ the conditional mean, which represents the expected value of Y for a fixed value of $Q_i$, and by $Q_{-i}$ all uncertain parameters except $Q_i$.

The first-order Sobol index, also known as the main effect index, cannot exceed one, and if it is close to zero we can claim that its input parameter have a small or no effect on the output. The higher the index value is, the more influence the associated parameter has on the resulting uncertainty. For the total Sobol indices, their sum is equal to or greater than one, and it is only equal to one if there is no interaction between the parameters [10].

3 Sampling Approaches

Sampling is the usual computational method used to estimate uncertainty measures, statistical moments and sensitivity as they have been described in the previous section. We can choose between several widely used sampling methods for a non-intrusive (or black box) approach: Quasi-Monte Carlo (QMC) [11] or Polynomial Chaos Expansion (PCE) [12]. This choice can be motivated by the numerical cost of such estimation, which is proportional to the number of samples and mainly depends on the number of uncertain input parameters.

3.1 Quasi-Monte Carlo

Monte Carlo is one of the most widespread method due to its simplicity and ease of implementation. The basic idea of the method is a purely random selection of sample evaluations without assumptions about the model. The number of samples is independent of the number of uncertain parameters, but a very high number is required in order to obtain reliable statistics, thus UQ estimations can rapidly become computationally expensive.

QMC methods are different approaches that aim at improving Monte Carlo method. They are based on variance-reduction techniques to reduce the number of model evaluations needed [11, 13]. In this work, instead of random selection of samples from $f_{Q}$, we use Sobol sequences which have a low-discrepancy sequence [8]. Using $N_s$ samples to evaluate the model $\mathcal {M}$, which gives the output results $Y = [Y_1 , Y_2, ..., Y_{N_s}]$, the mean and the variance can approximated by:

$$\begin{aligned} \mathbb {E}[Y] \approx \frac{1}{N_s} \sum _{i=1}^{N_s} Y_i, \end{aligned}$$

(8)

$$\begin{aligned} \mathbb {V}[Y] \approx \frac{1}{N_s-1} \sum _{i=1}^{N_s} (Y_i-\mathbb {E}[Y])^2. \end{aligned}$$

(9)

For the SA, we use Saltelli’s procedure [9, 14]. It is also based on Sobol sequences and to obtain a full set of main and total sensitivity indices, it is sufficient to take:

$$\begin{aligned} N_s=\frac{(d+2)N}{2}, \end{aligned}$$

(10)

where N is the number of samples required to get a given accuracy.

3.2 Polynomial Chaos Expansion

The PCE method is one of most competitive alternatives to (Quasi-)Random methods. The output Y is expanded into a series of orthogonal polynomials $(P_j)_j$ of degree p, which are functions of the input parameters Q and are chosen such that they are orthogonal to the input distributions [7, 15,16,17]:

$$\begin{aligned} Y \approx \widehat{Y}(Q) = \sum _{j=1}^{N_p} c_jP_j(Q). \end{aligned}$$

(11)

Here, $N_p$ is the number of expansion factors, which depends on both the degree p and the number of the uncertain parameters d, and is given by the binomial coefficient:

$$\begin{aligned} N_p = \left( {\begin{array}{c}d+p\\ d\end{array}}\right) ; \end{aligned}$$

(12)

And $(c_j)_j$ are the expansion coefficients, which can be estimated using tensor product quadrature or linear regression. These two methods that will be detailed later in this section.

By exploiting the properties of the orthogonal polynomials, the expectation value and the variance for the output model are:

$$\begin{aligned} \mathbb {E}[Y] \approx c_1, \end{aligned}$$

(13)

$$\begin{aligned} \mathbb {V}[Y] \approx \sum _{j=2}^{N_p}\gamma _{j} c_j \end{aligned}$$

(14)

Here, $\gamma _{j}$ is a normalization factor and it is defined as: $\gamma _{j} = \mathbb {E}[P_{j}^2(Q)]$.

For the sensitivity analysis, the first and total-order Sobol indices can also be computed using Sobol’s variance decomposition and by exploiting the polynomial chaos expansion properties [16, 18].

The estimation of the expansion coefficients is the most important step in the PCE construction. As mentioned before, we use two of the most common methods: spectral projection with tensor product quadrature and linear regression.

Spectral Projection. In this PCE variant, each expansion coefficient $c_j$ is estimated by projecting the solution onto the output space $\varOmega _Y$. This is analogous to the calculation of Fourier coefficients in the approximation of periodic function. By exploiting the orthogonality of the polynomial, Eq. 11 yields:

$$\begin{aligned} c_j = \frac{1}{H_j} \int _{\varOmega _Y} Y(q) P_j(q)f_{Q}(q)dq, \end{aligned}$$

(15)

where $H_j$ is a normalization factor associated with the polynomial $P_j$.

To compute the integral from Eq. 15, we use quadrature generators with weights $(\omega _k)_k$ and nodes $(q_k)_k$, and then we evaluate the outputs at those nodes. As the integrands are smooth polynomials, $c_j$ can be approximated with a good accuracy by the expression:

$$\begin{aligned} c_j \approx \frac{1}{H_j} \sum _{k=1}^{N_s}Y(q_k)P_j(q_k)\omega _k. \end{aligned}$$

(16)

The number of samples $N_s$ in this case is determined by the quadrature rule. For our applications, we use tensored quadrature with Gaussian schemes and in this case $N_s$ is:

$$\begin{aligned} N_s = (p+2)^d. \end{aligned}$$

(17)

As the required number of samples increases exponentially with the number of uncertain parameters, this method is beneficial only for a relatively small number of parameters.

Regression. This variant, also known as point collocation, uses a single linear least square approximation of the form:

$$\begin{aligned} \left[ \begin{array}{cccc} P_1(q_1) &{} P_2(q_1) &{} \ldots &{} P_{N_p}(q_1) \\ P_1(q_2) &{} P_2(q_2) &{} \ldots &{} P_{N_p}(q_2) \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ P_1(q_{N_s})&{} P_2(q_{N_s}) &{} \ldots &{} P_{N_p}(q_{N_s}) \\ \end{array}\right] \begin{bmatrix} c_{1} \\ c_{2} \\ \vdots \\ c_{N_p} \end{bmatrix} {}&= \begin{bmatrix} Y(q_1) \\ Y(q_2) \\ \vdots \\ Y(q_{N_s}) \end{bmatrix} \end{aligned}$$

(18)

The set of linear equations above results from the expansion of Eq. 11 at a set of collocation nodes $(q_k)_{k=1}^{N_s}$. To get stable nodes, we use pseudo-random samples generated from the input distribution $f_Q$ with specific rules as indicated in Chaospy [19]. To avoid an under-determined system, we should choose a number of samples such that $N_s \ge N_p$. According to [20], it is recommended that we use at least $2N_p$ samples. Therefore, we set:

$$\begin{aligned} N_s = 2\frac{(p+d)!}{p!d!}. \end{aligned}$$

(19)

In this case, we obtain a good least square approximation which becomes significantly more affordable than the full tensor product quadrature used in the spectral projection.

3.3 Implementation and Choice

All these sampling methods and variants are available, among others, in the EasyVVUQ² library; it provides a simple decomposition of UQ and SA as generic steps (sampling, encoder, decoder, collation, analysis) in a modular workflow approach. Sampling and analysis methods are mainly based on existing Python packages such as Chaospy [19], which provides probability distributions, statistics, quadrature generators and SA for PCE, and SALib [21] for SA with the QMC method. This decomposition allows users to exchange sampling methods easily without impacting the overall structure of the application.

In terms of sample sizes, we can use PCE with spectral projection as long as the number of uncertain parameters is low, typically smaller than 6. For a moderate number of parameters, PCE with regression becomes advantageous, and for a high number of parameters, quasi-Monte Carlo is the only method that is applicable on current computing platforms. In terms of accuracy, for smooth problems, PCE with spectral projection is expected to be more efficient compared with the regression variant or QMC. A more in-depth analysis and comparison between these methods and others can be found in [16, 19, 20].

4 Numerical Applications and Results

To perform UQ and SA on any multiscale application that couples single-scale models as a workflow, various level of intrusiveness can be considered. In a non-intrusive approach, either the entire workflow is treated as a black box, or only one single-scale model is being studied and treated as a black box. In a semi-intrusive approach, each single scale code is placed inside a black box, and UQ and SA are performed on the entire workflow.

In this section, we will present the two use cases for non-intrusive approach with the fusion multiscale application, and discuss their results as well as some performance considerations.

4.1 Multiscale Fusion Workflow as a Black Box

The multiscale fusion workflow³ was created to study the effects of turbulence happening at the micro time and space scales on the evolution of temperature profiles at the macro level of the fusion device. Currently, the workflow combines three main models as shown in Fig. 1: an equilibrium code that describes the plasma geometry, a turbulence code that approximates transport coefficient, and a transport solver that evolves temperature profiles at the macro time scale.

When considering the entire workflow as a black box, we use a simpler analytical model for the turbulence to reduce drastically the computational cost of the simulation. The extrinsic uncertainties we consider are:

• Boundary conditions defining the electron temperature at the plasma edge.

• Simplified heating sources, for which electron heating power Gaussian distributions are characterized by their amplitude, width, and position.

Fig. 1.

Sketch of the targeted fusion workflow diagram

We select the PCE method, described in Sect. 3, and assume that each of the uncertain parameters has a normal distribution in the range of $\pm 20\%$ around its original value. The workflow runs within a black-box for multiple time iterations until the plasma reaches steady state [2, 22]. The quantities of interest are the outputs Te and Ti, the electron and ion temperatures respectively, coming from the macroscopic transport model. The outcome from the EasyVVUQ analysis is plotted in Fig. 2.

Fig. 2.

Descriptive statistics and sensitivity analysis for electron and ion temperatures: on the left hand-side, the expected value and standard deviation are shown with respect to the normalized toroidal flux coordinate $\rho _{tor}$, and on the right hand-side the first-order Sobol indices are plotted for each of the uncertain parameters.

The standard deviation indicates that the ion temperature varies weakly since the uncertainties are carried by the electrons sources. The sensitivity analysis reveals that the variance in the electron and ion temperatures is mainly due to the uncertainty from three parameters: the position and amplitude parameters of the sources at core region of the plasma ($\rho _{tor} = 0$) and, as expected, boundary condition parameter at the edge region ($\rho _{tor} = 1$). The parameter width has no direct effect on the variance of the two quantities, so according to [23], this parameter can be neglected and then the number of samples can be reduced while keeping the same variance behavior.

We did the same simulation using uncertain parameters from heating sources and temperature boundary conditions for ions, and both electrons and ions. In all cases the UQ and SA results are qualitatively the same.

This use case was performed using both PCE variants: spectral projection and regression. They required 1296 and 140 samples, respectively, because we have 4 uncertain parameters and use quadratic polynomials. The results are qualitatively the same, and the differences between output means and variances of each method are at the order of $\mathcal {O}(10^{-6})$.

4.2 3D Gyrofluid Turbulence as a Black Box

In the workflow, given a plasma state and the geometry, the turbulence code provides the transport coefficients by computing the turbulent fluxes of particles and energy. Due to the chaotic nature of the turbulence model, these outputs are inherently noisy. And running such stochastic model is computationally expensive. Therefore, we begin the UQ study by isolating a 3D gyrofluid turbulence model that computes heat and particle fluxes [24] and consider one flux tube. We vary the electron and ion temperatures, Te and Ti, and their respective gradients, $\nabla Te$ and $\nabla Ti$, at the position of the flux tube.

We introduce four uncertain input parameters into the black box. Similar to the previous study, we assume that each parameter has a normal distribution in the range of $\pm 20\%$ around its original value. Finally, we perform the sensitivity analysis to explore how these parameters affect the resulting heat fluxes that determine the transport coefficients. It will also allow us to reduce the number of uncertain inputs that have negligible effect on the output variance [16, 25].

Fig. 3.

First-order Sobol indices electron and ion heat fluxes.

As indicated in Fig. 3, the variance in the electrons and ion heat fluxes are mainly due to the uncertainty in temperature gradients. However, we should be cautious about the impact coming from temperature uncertainties, especially Ti on ion heat flux.

In this case, we also use PCE method with spectral projection and with regression. As in the previous case, we also set the polynomial order to 4 and have 4 uncertain parameters, hence, the number of samples required are 1296 and 140, respectively. The parallel execution of the samples, whether the black box contains serial or parallel codes, is done via QCG PilotJob [5]. This middleware considers each sample as an independent job, but schedule them dynamically through a single batch allocation. For example, using 32 cores for each sample, we notice that the time to run UQ simulation with regression is 10 times less than the spectral projection.

5 Conclusion

In this paper, we presented results from the introduction of uncertainty quantification and sensitivity analysis methods in our multiscale fusion workflow application. We focused so far on two black box use cases: the entire coupled workflow with simple analytical model for turbulence, and singled-out 3D gyrofluid turbulence model. Integration of the non-intrusive methods was done with the EasyVVUQ library and the parallel execution of samples was performed in a single batch allocation with QCG-PilotJob. Both tools, available as part of the VECMA toolkit and in conjunction with the coupled workflow nature of our application, have proven to be simple to use and to allow a large variety of experiments with different methods and models.

A next step will be to use quantified uncertainties in order to improve the validation of our simulation results against experimental measurements. We also want to explore semi-intrusive methods, by coupling black boxes with propagation of uncertainties or by replacing the expensive 3D model with a cheap surrogate, in order to reduce the computational cost for quantifying the uncertainties such as it can be performed routinely in our simulation. Other optimizations (e.g. parallelization of the sampling step with a closer EasyVVUQ–QCG-PilotJob integration) and scalability tests are being performed and will be the subject of a more in-depth performance study.