Building Confidence in Simulation: Applications of EasyVVUQ

Validation, verification, and uncertainty quantification (VVUQ) of simulation workflows are essential for building trust in simulation results, and their increased use in decision‐making processes. The EasyVVUQ Python library is designed to facilitate implementation of advanced VVUQ techniques in new or existing workflows, with a particular focus on high‐performance computing, middleware agnosticism, and multiscale modeling. Here, the application of EasyVVUQ to five very diverse application areas is demonstrated: materials properties, ocean circulation modeling, fusion reactors, forced human migration, and urban air quality prediction.


Introduction
In order for the results of computational science to become widely accepted components of decision making processes, such as in medicine and industry, it is essential that we quantify the trust one can have in the model in question. Confidence can only be gained by ensuring not only that simulation codes are solving the correct governing equations (validation), but they are solving them correctly (verification) and we have a comprehensive estimate of the uncertainties in the result uncertainty DOI: 10.1002/adts.201900246 quantification). [1,2] Collectively the processes involved in evaluating our level of trust in the results obtained from models are known as VVUQ. While the need for rigorous model assessment is widely acknowledged, it is far from being universally implemented within the scientific literature. The reasons for this are wide ranging, but include lack of specialist knowledge of VVUQ techniques and, until recently, the difficulty in obtaining sufficient computational power to perform the necessary sampling in large scale simulations. We have recently developed EasyVVUQ, [3] a package designed to help leverage recent advances in the scale of computational resources to make state of the art VVUQ algorithms available and accessible to a wide range of computational scientists. EasyVVUQ is a component of the VECMA open source toolkit (http://www.vecma-toolkit.eu), which provides tools to facilitate the use of VVUQ techniques in multiscale, multiphysics applications. [4] In order to enable straightforward computations of EasyVVUQ scenarios on HPC resources, the tool has been designed to work with a variety of middleware technologies, such as FabSim3 [5] or www.advancedsciencenews.com www.advtheorysimul.com QCG. [6] The integration with pilot job mechanisms, in particular with QCG-PilotJob [7] and Dask JobQueue, allowed us to bypass limitations of regular queuing systems related to the scheduling of workloads composed of a very large number of relatively short tasks.
In this paper, we showcase the use of EasyVVUQ in a selection of applications chosen to have highly varied computational and VVUQ requirements. The examples come from a wide range of domains; materials science, climatology, fusion physics, forced population displacement, and environmental science. All of the examples come from active research projects and were chosen to highlight the range of capabilities of EasyVVUQ: 1. Materials-a simple parameter sweep performed using a computationally expensive molecular dynamics simulation; 2. Ocean circulation-estimation of Sobol sensitivity indices using stochastic collocation in a differential equation based model; 3. Fusion-estimation of Sobol sensitivity indices using the polynomial chaos expansion in a multiscale simulation workflow; 4. Forced migration-estimation of Sobol sensitivity indices in an agent based model; 5. Environmental-estimating uncertainties using stochastic collocation in a model forecasting urban air quality.

EasyVVUQ
EasyVVUQ is a Python library, developed within the VECMA project, designed to simplify the implementation of VVUQ work-flows in new or existing applications. The library is designed around a breakdown of such workflows into four distinct stages; sampling, simulation execution, result collation (or aggregation), and analysis. In the sampling stage, the uncertainty on the inputs of the model are defined, for instance, by specifying independent probability density functions p( i ) for each model parameter i . This leads to a sampling plan, that is, a collection of points in the input space where the model must be executed. This execution stage is deemed beyond the remit of the package (it can be handled for instance by Dask JobQueue, FabSim3, [5] QCG-PilotJob [8] , RADICAL Cybertools, [9] etc.) but EasyVVUQ does provide some functionality to address it. The final collation and analysis stages, which are handled by EasyVVUQ, deal with post processing the simulation outcomes into mean predictions, uncertainty estimates, and sensitivity measures. A common object, the Campaign, contains information on the application being analyzed alongside the runs mandated by the sampling algorithm being employed, and is used to transfer information between each stage. All applications outlined below share a similar Campaign creation step, up until the point where a specific sampler and input uncertainties are selected. This general procedure consists of creating an EasyVVUQ Campaign object, defining the parameter space and code outputs, and selecting an encoder, decoder, and collation element. The following code can be used as a generic template for all applications we consider (up to sampler selection), where variables indicated by < ⋅ > have to be replaced with application-specific values.
www.advancedsciencenews.com www.advtheorysimul.com Most such variables are self explanatory, hence we only highlight • "<path_to_input_template>": This is the path to the input template of a particular application. Essentially, this is just the standard input file that the application uses, except that the value of uncertain variables (those in params), must be flagged by a delimiter ($ in this case, e.g., param_1=$param_1), such that they will change their values for each sample. • "<input_filename>": This is the file name that will be given to each realisation of the input template.
Note also that the parameter space definition (in Listing 1) has optional specification of the type and minimum/maximum allowed values. EasyVVUQ's Cerberus dependency uses this information to apply verification of input variables such as type, range, and conditional checks. [10] EasyVVUQ additionally provides version checking for the library (and each of the component VVUQ elements) so that the user is made aware when a given element they have been using in the past may now have a new algorithm/behavior. This functionality, along with detailed logging of element application and "fail-early" checks, is intended to aid the user in verifying that a VVUQ workflow is doing what was intended.
For more information on the various EasyVVUQ elements, we refer to the software release publication. [3] The following sections detail the use of EasyVVUQ as applied to a variety of different application domains. All can be considered to have gone through the Campaign creation process as described above, hence we will not repeat this setup code. Only the assignment of input distributions, sampler selection, and post-processing will be described for each example application. Relevant information, code, and output data for the following example applications may be found in Supporting Information. [11]

Goals
We intend, through the five following example sections, to demonstrate how EasyVVUQ can be used to augment existing applications with VVUQ features or capabilities, notably: A) In a non-intrusive manner (all solvers may be used as "black boxes", with no changes to their internals). This applies to all five example applications. B) Favouring consistency and interoperability between approaches (a particular UQ approach may be painlessly swapped for another due to EasyVVUQ's standard interface for VVUQ elements). In this work, we demonstrate a basic parameter sweep (Section 5), stochastic collocation (Sections 6, 8, and 9), and polynomial chaos expansion (Section 7), showing a similar pattern of application. C) Combining VVUQ elements together into single elements (to create complex) behavior easily using small, existing parts). This is demonstrated with Encoders in the Fusion (Section 7) and UrbanAir (Section 9) applications. D) Allowing execution of generated runs in any order, using any desired middleware (of particular importance to HPC applications, where job submission and execution patterns are key to performance and highly dependent on the computing resources). This design principle is demonstrated by the mix of execution methods used in the five example applications, ranging from local or manual execution through to dynamic pilot job schedulers.
The focus is not on the scientific results of each section, but on the consistency of the approach when applied to different techniques, for different solvers from different scientific domains. EasyVVUQ seeks to abstract out both the underlying model (with its application specific inputs and execution needs) and the implementation of VVUQ (particularly UQ) algorithms. These algorithms, which may be custom implementations in EasyVVUQ or sourced from existing libraries such as chaospy or SAlib, all interact via standardized interfaces, such that the user should not have to worry about the provenance of the underlying implementation, but rather about connecting the operations together (or swapping them for others).

Background in UQ Techniques
This section gives a brief background in the various UQ techniques which are used in the example applications.

Stochastic Collocation
Once an input distribution is defined, the output quantities of interest (QoIs) become random variables. The stochastic collocation (SC) method creates a polynomial approximation of a quantity of interest q in the stochastic space ∈ ℝ d via the following expansion: Here, the stochastic space is the space of the uncertain code input parameters, for which independent, user-specified, probability density functions (pdfs) must be provided: i ∼ p( i ), i = 1, ⋯ , d. Furthermore, q j are the code samples which are computed on a structured multi-dimensional grid, and N p is the total number of collocation points, that is, the total number of code evaluations. The samples q j are interpolated to an arbitrary point within the stochastic space by means of Lagrange interpolation polynomials L( ). For interpolation in multiple dimensions (d > 1), L( ) is built as a tensor product of 1D Lagrange polynomials. The SC method, and similarly the polynomial chaos expansion (PCE) method (described briefly in the next section), are well-known and we refer to ref. [12] for more details on these techniques. Suffice it to say that the tensor product construction yields an exponential increase in N p with the number of uncertain variables d and the chosen polynomial order, an example of the familiar "curse of dimensionality." However, for moderate values of d, the SC and PCE methods can display exponential convergence with N p , thereby outperforming Monte Carlo sampling. [12] There are three main uses of the SC expansion (1). First, the N p code samples q j can be used to estimate the first two moments of q in the stochastic space, giving a mean prediction and an estimate of the output uncertainty due to the prescribed distributions on the inputs. Second, (1) acts as a computationally inexpensive surrogate model for the code. Using the Lagrange polynomials, the code samples q j (evaluated at a specific parameter values j ), can be interpolated to an unsampled location . Finally, the SC expansion is amenable to variance-based global sensitivity analysis. Estimates of the well-known Sobol sensitivity indices can be obtained from Equation (1) as a post processing step, (which is outlined in Section 4.3).

Polynomial Chaos
The PCE method is an expansion technique that is closely related to SC method presented in Section 4.1. Whereas in SC we build Lagrange interpolation functions for known coefficients, in PCE we estimate coefficients for known orthogonal polynomial basis functions. Here, we can approximate the quantity of interest q with the following expansion: In this equation, , c j , and N p are the uncertain parameter, expansion coefficients, and number of expansion factors, [13,14] respectively. The polynomials P j are chosen such that they are orthogonal to the input distributions, which differ from the SC expansion in Equation (1).
To compute the c j coefficients, two variants have been implemented: spectral projection and linear regression. In the spectral projection variant, we project the response against each basis function (composed of the polynomials set (P j )) and we exploit their orthogonality properties to extract each coefficient. In the linear regression variant (also known as point collocation), we use a least squares method that minimizes a normed difference between the PC expansion and the output for a set of samples; the coefficients c j are then the solution of the resulted linear system. [12] By using the PCE method, like the SC method, we can obtain the statistical moments (mean, standard deviation, variance, and (100 − x) th percentile) of the quantities of interest, and we can also provide a global sensitivity analysis in the form of Sobol indices (which is outlined in the next section).

Sobol Indices
Sobol indices are variance-based sensitivity measures of a function q( ) with respect to its inputs ∈ ℝ d . [15] As in the case of the SC method, an independent probability density function p( i ) is assigned to each input j , which makes this a global method. Local sensitivity methods, on the other hand, measure the sensitivity of q at some point 0 in the domain, and are uninformative away from this point. Another advantage of global methods is that they can capture the sensitivity due to higher-order interactions (several parameters changing at once).
Sobol indices are derived from the analysis of variance (ANOVA) decomposition of q( ). This decomposes q into a sum of basis functions of increasing input dimension, which in long forms reads as: A more concise notation is www.advancedsciencenews.com

www.advtheorysimul.com
where u is a multi-index and  is the power set of  := {1, 2, ⋯ , d}. Let us define u as u := { i |∀i ∈ u}, that is, the set of all inputs with an index in u. Furthermore, u ′ is the complement of u, that is, u ∪ u ′ :=  and u ∩ u ′ = ∅.
In the ANOVA decomposition, the basis functions q u satisfy the following properties: that is, they have zero mean and are orthogonal when integrated over the distributions. These properties hold when the basis functions are defined as It is perhaps more clear to write this in terms of conditional expectations: Hence, q ∅ represents the mean of q( ), and the q i basis functions represent the effect of varying a single parameter i , minus the mean. Basis functions such as q ij capture the effect of changing i and j simultaneously, minus all lower-order interactions, etc.
Therefore, the variances of these basis functions are the sensitivity measures we aim to approximate. Since the q u have a zero mean, these are defined as Using the orthogonality property of the basis functions, (8) can be rewritten as Expression (9) allows us to compute all D u in increasing order, if we can compute the first integral on the right-hand side. The authors of ref. [16] developed a method to approximate this integral using the SC expansion (1) for q( ), and similar techniques exist for the PCE method. [17] It is out of the scope of the current paper to go into detail, and we refer the interested reader to refs. [16,17] for the mathematical details. Essentially, once all the code samples are obtained, the Sobol indices, which are defined as can be approximated in a post-processing step. Here, D := ar[q] = ∑ u⊆ D u . [15] Note that all D u are positive, and that the sum of all possible S u equals 1. Each D u measures the amount of variance in the output q that can be attributed to the parameter combination indexed by u.

Application Outline
Molecular dynamics (MD) simulations are often used to investigate the properties of materials, [18] including as part of multiscale material prediction applications. [19] Here, we take a well understood soft-matter system and study how calculations of its Young's modulus (stiffness) using MD can vary with the system size and starting configuration.
The system under consideration (Figure 1) is an epoxy resinepoxy tetraglycidyl methylene dianiline (TGMDA) cured with polyetheramine (PEA) in a 1:1 ratio. Epoxies are thermosetting polymers. Small reactant monomer molecules have several reactive sites which create strong covalent bonds between several other molecules, forming a dense network of crosslinks. The resulting polymer network is very strong and epoxies are widely used in manufacturing, in the aerospace industry, as adhesives, and as multipurpose insulators.
MD simulations in the condensed phase are almost always periodic, which means that we only simulate a comparatively small simulation cell to approximate the bulk properties. The size of this simulation cell has many implications for computational cost and, more importantly, the scientific results it furnishes. Finite size effects, self-interaction across periodic boundary boundaries, and thermal fluctuations in small systems can all affect the simulation's outcome. We measure the Young's modulus (YM) of an epoxy system by measuring the pressure exerted along one axis before and after a small strain. www.advancedsciencenews.com www.advtheorysimul.com

VVUQ Algorithm
Since the instantaneous pressure of a molecular dynamics simulation can fluctuate by several GPa, it is necessary to average this value over a long sample period to measure the change in pressure due to an applied strain. The YM could also be affected by starting velocities of the atoms in the system, and the configuration of the epoxy network. To measure the system size dependence of all of these potential sources of variance, we design an EasyVVUQ Campaign that will take samples across each variable. Then, bootstrap analyses will measure their effect on the YM. A closer look at the variance due to each variable will show which is the most significant.

Execution Pattern
This application makes use of the BasicSweep sampler in EasyVVUQ, which recursively carries out a parameter sweep across the range of allowed values specified for each input variable. The system size is limited by computational cost at the high end, and the system stability at the lower. In this case, we know these approximate limits beforehand, so choose the specific range we want to sample using this method. The sampler is set up like this: In the above, we create a sampling element to sample across 6 simulation sizes, build 10 epoxy networks (structures) at each size, then measure the YM for each network starting from 10 different snapshots. These numbers are somewhat arbitrary, and more parameters could be swept depending on availability of computational resources.
Building the epoxy networks is done with an in-house developed script, [20] used in ref. [21] Simulating the epoxy network is accomplished using LAMMPS. [22] The execution of the system building procedure and measurement simulations are submitted on a remote computing resource. The "restart campaign" functionality of EasyVVUQ is required here, as the sampling and analysis stages were performed in separate python scripts. This ability to restart a campaign from a different script is useful in cases where, for example, the runs are expected to take a long time on a remote computing resource, and the user cannot or does not wish to have an EasyVVUQ script running locally, waiting for such jobs to finish.
Each replica generates three values for the YM, measured by separately straining along each principle axis of the polymer simulation. So that we can treat each of these values as equivalent measurements, we change the results pandas DataFrame format to have one row for each value; all YM values are in one column which makes some analysis more straightforward. This mode may be set using

Results and Analysis
The system can be characterized by simply employing a boostrap analysis of the campaign.
The results of the above analysis are shown in Figure 2, along with a histograms of all measured YMs associated with each box size. We can clearly see that the average YM is independent of simulation size, above 4 nm. There is approximately a 20% increase in YM for a system size of 3 nm. We can safely say that for this system the characteristic length is therefore less than 4 nm. We would like to know if the structure of an epoxy network has a significant bearing on the YM of a system, that is, if there is a large variation in the expected YM given an epoxy network. We approach this with the law of total variance where "#" is used to denote a specific network of cross-links. We can calculate the first and third terms of this law (the total variance in YM, and the expected variance given a specific structure) by some straightforward manipulation of the campaign results DataFrame.
Detailed results are shown in Supporting Information; [11] however, the analysis shows that Var[E(YM|#)] ≪ E[Var(YM|#)]. Therefore, the epoxy network structure has no significant effect on the YM. The variance is due to the inefficient sampling of MD. We studied low strains (0.5%) because epoxies are often brittle above these strains, but simulating further (into plastic deformation) could resolve the dependence on the network structure. Chaotic dynamical systems may manifest a pathology of IEEE floating point arithmetic which was hitherto unknown, [23] providing a potentially interesting overlap between uncertainty quantification and verification in affected systems.

Application Outline
In this section, we consider the forced-dissipative vorticity equations for 2D incompressible flow (as described in Verkley et al. [24] ), used as a simplified study for the general circulation in the oceans. The governing equations are Here, is the vertical component of the vorticity, defined from the curl of the velocity field V as := e 3 ⋅ ∇ × V, where e 3 := (0, 0, 1) T . The stream function Ψ relates to the horizontal velocity components by the well-known relations u = − Ψ∕ y and v = Ψ∕ x. The non-linear advection term is defined as This system generates flow fields such as those shown in Figure 3, which depicts a snapshot of the vorticity . www.advancedsciencenews.com www.advtheorysimul.com As in ref. [24], the forcing term is chosen as the single Fourier mode F = 2 3∕2 cos(5x) cos(5y). The system is fully periodic in the x and y directions over a period of 2 L, where L is a user-specified length scale, chosen as the Earth's radius (L = 6.371 × 10 6 [m]). The inverse of the earth's angular velocity Ω −1 is chosen as a time scale, where Ω = 7.292 × 10 −5 [s −1 ]. Thus, a simulation time period of a single "day" can now be expressed as 24 × 60 2 × Ω ≈ 6.3 non-dimensional time units. Given these chosen length and time scales, we non-dimensionalize (12) and solve by using a spectral method with the second-order accurate AB/BDI2 time-stepping scheme. [25] The viscosity and the forcing term coefficient are tunable parameters, and are typically set to a value such that the build up of grid-scale noise at the smallest resolved scale is prevented. In our example code, their values are computed such that a Fourier mode at this scale is exponentially damped with a user-specified e-folding time scale, that is, a time scale over which a decay of 63 % occurs (1 − e −1 ). This leads to the following expressions for and : Here, K is the highest resolved wave number in our spectral method, which is fixed at 85. More important for our current discussion are 1 and 2 , that is, the aforementioned damping time scales (expressed in days), which we treat as uncertain. We use EasyVVUQ to estimate the effect of this uncertainty on certain measures derived from the solution of (12). Our focus will be on the (time-dependent) energy E and enstrophy Z, defined as Specifically, we are interested in the time-averaged statistical moments of the energy E and enstrophy Z; for example, our quantities of interest q take the form of and likewise for the enstrophy. The integration interval [T 0 , T] will be defined later.

VVUQ Algorithm
For this particular problem, we will use the stochastic collocation method, as outlined in Section 4.1. In addition to the statistical moments of the aforementioned energy and enstrophy, the Sobol sensitivity indices of our damping time scales will serve as our QoIs as well. Specific implementation details are given next.

Execution Pattern
EasyVVUQ is designed to work with the Chaospy library, [26] for the specification of the input distributions. We will assume the following distributions for the uncertain decay times associated with and : That is, we assume that the viscous term ( ∇ 2 ) in Equation (12) has a uniformly distributed uncertain decay time at the smallest retained scale between 1 and 5 days, whereas our forcing term is damped somewhere between 85 and 95 days. We then select the stochastic collocation sampler via By selecting a polynomial order of 6, a seven-point quadrature rule for each uncertain dimension is created. Hence, since we have two uncertain variables, we obtain a tensor grid of 49 points in the stochastic space, see Figure 4a. At each point, we have to evaluate the code solving (12). Instead of directly creating a full tensor product of the seven-point 1D quadrature rule, we can also construct a sparse grid (see Figure 4b), which uses a linear combination of tensor products of quadrature rules of different orders. [12] By using carefully chosen 1D quadrature rules, many points in the different tensor products will coincide, leading to a more efficient sampling plan in high dimensions. To switch to a sparse grid, one might use Here, quadrature_rule="C" denotes the use of 1D Clenshaw-Curtis quadrature rules, which are a common choice in sparse grid constructions. Furthermore, growth=True selects an exponential growth rule, which ensures that the Clenshaw-Curtis rules are "nested" such that a quadrature rule of the next order contains all points of the previous order, leading to the aforementioned more efficient sampling plan in high dimensions. However, since we just have two uncertain variables here, we will use the full tensor product construction. Depending on the spatial resolution of the computational grid (in our case, we employ a 2D grid of 256 × 256 points), the cost of sampling Equation (12) at all collocation points j can be high. Moreover, since we are interested in the time-averaged statistics as in Equation (16), we must run each sample until convergence in these statistics can be safely demonstrated. We use FabSim3 [5] to facilitate the execution of these samples in parallel on the Eagle supercomputer at the Poznan Supercomputing and Networking Center (PSNC). A FabSim3 plugin "FabUQCampaign" has been created to execute the ensemble run of EasyVVUQ samples on a remote resource, with minimal change in the code that is executed on the localhost. For a tutorial on the setup of FabUQCampaign, see ref. [27].
The fab module is a wrapper around FabSim3 command-line instructions, such that these can be executed from within Python. Furthermore, machine specifies the name of the remote HPC resource (the PSNC Eagle cluster in our case), and campaign_dir is the directory containing the EasyVVUQ campaign. Finally, "ocean" is the name of the script which executes a single run of our model (12); see the tutorial [27] for more details. Once the ensemble run has completed, the results can be retrieved through: If one wishes to run a (small) local ensemble for testing or debugging purposes, specifying machine="localhost" will make sure that everything is executed locally. Note that FabSim3 is not the only available execution interface between EasyVVUQ and HPC clusters. EasyVVUQ-QCGPilotJob is a lightweight integration code that simplifies usage of EasyVVUQ with a QCG-PilotJob execution engine; see ref. [7] for a tutorial.

Results and Analysis
In our example, d = 2 ( and ), and our quantities of interest are time-averaged moments of Equation (16) of the energy and enstrophy. For each sample of the ensemble run, Equation (12) is simulated for 11 years, and the last 10 years are used to compute the time-averaged E and Z moments. To perform the postprocessing analysis of these samples, an SCAnalysis object is created: The results dictionary contains the statistical moments and the Sobol indices of the quantities of interest, the latter of which are given below for this particular case: A value close to one means that this variable, or combination of variables, explains most of the variance in the selected output. Clearly, is the more influential parameter for both the timeaveraged energy E and enstrophy Z. However, for the corresponding standard deviations (stdevs), does play an important role in the second-order Sobol index, indicating a significant interaction between and for these QoI.
In order to use the SC expansion as a surrogate model for the real code, we can draw random samples from Equation (1) via Here, xi is an array containing a random sample from the input distributions of and . The surrogate is far cheaper than the original model, such that we can use it to evaluate the output probability density function via a kernel-density estimate (KDE). Figure 5 shows the KDE of E, evaluated using 5 × 10 4 samples from (1).

Application Outline
Thermonuclear fusion is potentially a solution to the provision of base load electricity, which is carbon free and not subject to geopolitical problems. Understanding the mechanisms of heat and particle transport in hot fusion plasma is one of the keys to obtain a cost-efficient reaction in the fusion devices. Our present understanding of the problem is that turbulence at small scales is responsible for much of this transport, but the profiles of temperature and density evolve over much larger scales.
A wide standardization effort toward integrated modeling [28] for fusion plasmas has allowed us to build modular applications in the form of a workflow. The code-to-code coupling is done via standardized data-objects [29] (referred to hereafter as CPO files), while specific parameters are stored in XML. This setup allows users to swap codes with others of different complexity. Based on this effort, a multi-scale application is developed to study the turbulence effects on plasma transport at larger scales. [30] However, much remains to be done on the validation of such simulations as well as on the control of their uncertainties. In this work, we present an early validation pattern we uncovered by extracting and comparing experimental and UQ simulation output distributions. In our application, these uncertainties originate from applied heating sources (extrinsic) and/or from the noisy, chaotic nature of the turbulence (intrinsic). We focus here on quantifying extrinsic uncertainties for the heating source as well as boundary conditions for both electron and ion temperatures. The heat source (energy per unit time) for each species is a Gaussian function with respect to the radial (or toroidal flux) coordinate tor , and it is characterized by its amplitude, width, and position. The boundary conditions refer to the initial temperatures at the plasma edge for both species; the edge is positioned at the maximum tor value, or at normalized tor = 1.0.

VVUQ Algorithm
The EasyVVUQ library provides both quasi-Monte Carlo (QMC) and polynomial Chaos expansion (PCE) methods (described in Section 4.2) that we can select from to conduct UQ and SA in the multi-scale fusion workflow. [31] In the work we present here, the PCE method was selected because it can carry out the calculations much faster than the QMC method. However, this is only valid if the number of uncertain parameters remain relatively low.

Execution Pattern
Similar to the ocean circulation example, we specify the input distributions using chaospy through EasyVVUQ. In addition, to fully benefit from the standardized interface for each code within our multi-scale workflow, we extended the EasyVVUQ base encoder with a new domain specific CPOEncoder (for boundary conditions of electron and ion temperature profiles) and a generic XMLEncoder (for electron and ion heating sources approximated by the amplitude, position and width of a Gaussian function). These format-bound encoders allow us to update real data and parameter files without having to create a template, which in turn gives us more flexibility. Since we are interested in uncertainties driven by both the heating sources and the boundary conditions for electrons and ions temperatures, we need to combine these two encoders with the MultiEncoder provided by EasyVVUQ. Therefore, the encoder creation from listing 1 is modified with the following snippet of code: www.advancedsciencenews.com www.advtheorysimul.com • common_dir is a folder that contains all required input files.
• uncertain_params is a python dictionary specified by the user, and it contains the list of parameters with their probability distributions types followed by the chaospy glossary.
In addition, the new encoders have a specific function that provides two dictionaries containing the names and types of all parameters to be varied and their corresponding distributions.
We set up the PCE sampler using a polynomial of order 4 to ensure good accuracy: The output of the application is composed of several CPO format files, so the same kind of modification is done for the creation of the decoder, which uses our domain-specific CPODecoder.
Finally, to generate all samples needed for the analysis, we can either call the function ExecuteLocalexecute provided directly by EasyVVUQ or resort to a wrapper enabling the execution using the QCG-PilotJob mechanism. [7,8]

Results and Analysis
To perform a post-processing analysis on the generated samples, we use the PCEAnalysis object from EasyVVUQ. For the results, as in the ocean circulation example, we use analysis_results, the output dictionary of the campaign object's get_last_analysis() method, in which the statistical moments, and the Sobol indices of the Quantities of Interest are stored.
The current version of the fusion workflow uses an analytical turbulence code, with four uncertain parameters (amplitude, width, position of heating source, and boundary condition). We assumed each of these parameters has a normal distribution in the range of ±20% around its original value, and as the number of samples is determined by the uncertain parameter number and polynomial degree in the PCE method, the number of runs required for this example is 1296.
The uncertainty quantification of the fusion workflow is shown in Figures 6 and 7. The quantities of interest are the electron and ion temperature profiles, spanning from the radial position of plasma core ( tor = 0) to the edge (normalized tor = 1.0). 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 and, as expected, boundary condition parameter at the edge region. The parameter width has no direct effect on the variance of the two quantities, so according to ref. [32], this parameter can be neglected and then the number of samples can be reduced while keeping the same variance behavior.
In addition to uncertainty quantification, the fusion application performs validation on the simulation results by comparing the distribution of the QoI to the distribution from experimental measurements (first results are shown in the Figure 8). Specifically, we create the ValidationSimilarity object to determine the similarity between two distribution functions: www.advancedsciencenews.com www.advtheorysimul.com   For the similarity measure, we use the Jensen-Shannon distance (JSD), which is a symmetrized and smoothed version of the Kullback-Leibler divergence. [33,34] It is defined by Here, Ns is the number of samples, P and Q are defined as two discrete probability distributions, and M = 1 2 (P + Q). As presented in Figure 8, Jensen-Shannon distance takes values in the range [0, 1]. The values closer to 0 indicate a smaller "distance" between the two distributions and therefore a stronger similarity. Two other measures based on the Hellinger and Wasserstein distances [33] are also available in EasyVVUQ. These measures were also tested on the current example, and they give equivalent results as the Jensen-Shannon distances.

Application Outline
Forecasting forced displacement is of considerable importance since 70.8 million people are today being forcibly displaced worldwide, a record level. [35] It is also challenging as many forced population data sets are small and incomplete, and data sources have too little information. [36] Nevertheless, forced population predictions are essential to save the lives of such migrants, to investigate the effects of policy decisions and to help complete incomplete data collections on forced population movements.
Through the use of computational approach, namely the FLEE agent-based simulation code, we predict the distribution of forced population arrivals to potential destinations as governments and NGOs can efficiently allocate humanitarian resources and provide protection to vulnerable people. [37] We represent forcibly displaced people as individual agents, combining simple rulesets for individuals to allow complex movement patterns to emerge and validate simulation results against real data. We are also able to systematically explore the possible impact of policy decisions using the FabSim3-based FabFlee toolkit while accounting for the sensitivity to a subset of parameters and assumptions in the model, such as the probability of migrants making specific moves. In Figure 9, we present a simulation instance of the Mali conflict, in which a number of insurgent groups began a fight for the independence of the Azawad region resulting in an increasing number of forcibly displaced people since January 16, 2012. [36,37]

Figure 9.
Overview of geographic network model for Mali, which includes conflict zones (red circles), camps (dark green circles), forwarding hub (light green circle), and other major towns (yellow circles) interconnected with straight-lines that represent roads and their length in kilometers with adjacent blue numbers.

VVUQ Algorithm
FabFlee uses the EasyVVUQ library to facilitate VVUQ for simulation analysis. It allows us to automate parameter exploration analysis and explore essential one-at-a-time input uncertainty quantification. Importantly, uncertainty quantification and sensitivity analysis are required in multiscale migration studies to understand in what regimes and scenarios our simulation approach performs well. FabSim3, EasyVVUQ, QCG-PilotJob, and other QCG components can be combined in a variety of ways, enabling users to combine their added values while retaining a limited deployment footprint. As previously mentioned, EasyVVUQ can use FabSim3 to facilitate automated execution. Users can convert their EasyVVUQ campaigns to FabSim3 ensembles using a one-line command, and the FabSim3 output is ordered such that it can be directly moved to EasyVVUQ for further decoding and analysis.

Execution Pattern
We use similar approach as described in the ocean circulation example for sensitivity analysis of forced migration application. In particular, we analyze the probability of parameters when agents move from their current location to a different one on a given day. These probabilities depend on the type of locations, namely conflict zone, camp, or other location, where agents reside. [37] We adjust these parameters to understand the importance of our assumptions in regard to the validation results. In Figure 10, we provide the overall workflow of forced displacement application for sensitivity analysis.
To provide the input distributions, for instance, we specify a uniformly distributed move chance probabilities for camps between 0.0001 and 1.0, as well as for conflict locations between 0.1 and 1.0 as illustrated below.
Then, we set up the stochastic collocation (SC) sampler using www.advancedsciencenews.com www.advtheorysimul.com  where a polynomial order is 3 in this instance. In turn, it creates a four-point quadrature rule for each move chance parameter (see Figure 11). EasyVVUQ encodes the generated samples to FLEE input definitions for specific conflict simulations and submits all ensemble runs for execution using FabSim3 to Eagle machine where QCG-PilotJobs schedule submitted ensemble runs and pre-reserved resources.

Results and Analysis
We apply the same SC-based Sobol index method [16] as in the ocean circulation example above to the Mali conflict, and obtain the results illustrated in Table 1 for two parameters. We draw our own distinction on the Sobol indices by accepting parameters with values below 0.05 while identifying parameters with higher values as sensitive to output results. The camp_move_chance parameter is more sensitive in our model compared to the other parameter, namely conflict_move_chance, since camps are primary destination locations for forcibly displaced people fleeing from conflict locations. We also find that our models are not sensitive to the combination of these parameters.

Application Outline
The UrbanAir application concerns the modeling and forecasting of the concentration and dispersion of pollutants. It is a 3D multiscale model that combines a numerical weather prediction (NWP) model, running at larger scale (e.g., mesoscale), with a city-scale geophysical flow solver for accurate prediction of contaminant transportation through the street corridors, over buildings and obstacles. The NWP model is based on the community Weather Research and Forecasting (WRF) model, [38] while the city-scale problem is solved using the EULAG model. [39] EULAG is a numerical solver for all-scale geophysical flows, with many proven scenarios, for example, flows around buildings [40] with comparison against wind tunnel experiments. [41] The coupling between WRF and EULAG model has been evaluated in ref. [42]. Typically, an emergency response situation requires fast and accurate tools. However, the use of more complex and expensive models is dictated by the need for accurate prediction of peak concentrations and plume temporal evolution.
With increased model resolution, small-scale flow characteristics are becoming more essential for prediction, and general urban parameterization coming from the NWP model is not enough. The WRF output is used as the initial and lateral boundary conditions for the EULAG simulation, along with terrain data (terrain elevation, road network, buildings shapes, and height) and emission data. Figure 12 presents the general workflow of the application. The IMB approach is used in EULAG to explicitly resolve complex building structures, accounting for different urban aerodynamic features, such as channeling, vertical mixing, and street-level flow. The pollutant dispersion is simulated using passive tracer equations.
The NWP model may be supplemented with an additional chemistry module, to simulate chemical transportation and mixing over larger scales. [43,44] In order to accurately simulate at small scales (grid resolutions up to 1 metre), HPC resources are required. EULAG is proven to scale up to thousands of CPU cores to support such resolution and to decrease overall time-to-solution. [45] The key problem in providing accurate forecasts is the lack of complete, well-known emission sources. Contaminants-such as NO, NO 2 , PM2.5 (particulate matter under 2.5 µm in size), and PM10 in particular-are emitted by point sources (e.g., industrial chimneys), line sources (e.g., road transportation), and area Adv. Theory Simul. 2020, 3, 1900246 www.advancedsciencenews.com www.advtheorysimul.com sources (e.g., heat appliances). The uncertainty comes from unknown emission details. Taking road transportation as an example, there is a set of parameters that need to be estimated: these include the ratio of cars using gasoline to diesel fuel, fuel usage, emission index, percentage of cars that cold-started, and so on. Through the use of computational ensemble simulations, we can address these issues using statistical data, such as by combining the number of cars passing a given road section within 1 h with previously estimated parameter values.

VVUQ Algorithm
In order to assess the influence of unknowns in the emission sources, we have designed an EasyVVUQ campaign that samples across each of the input variable. It allows us to assess input uncertainty quantification and sensitivity analysis, though we concentrate on the former at the moment. The uncertainty may additionally stem from weather boundary conditions, heights of buildings, etc. To facilitate uncertainty quantification for this computationally demanding application, The QCG-PilotJob is used to choreograph the execution of the ensembles.

Execution Pattern
Currently, we focus on quantifying uncertainty coming from parameters related to NO 2 emission attributed to road transportation. The simulations require input data regarding, for example, NO 2 index from gasoline engines, fuel usage, density, and ratio of gas to diesel cars. The input distribution is specified using chaospy via EasyVVUQ. Here, we focus only on parameters related to petrol-powered vehicle, while a similar setup is needed for diesel vehicle analysis.
Next we set up different samplers for different input parameters we want to be sampled.
A custom encoder is used, EmisEncoder, whose goal is to use the values of the sampled parameters as components to calculate the correct value of road transportation emissions. www.advancedsciencenews.com www.advtheorysimul.com We setup a stochastic collocation sampler and use the multiencoder for our campaign.
To facilitate running ensembles, each of which requires hundreds of cores, we use an integrator between QCG-PilotJob and EasyVVUQ called EasyVVUQ-QCGPJ. [7,8]

Results and Analysis
In this simulation example, a 2 × 2 × 2 metre grid resolution has been used, and the same resolution has been applied to the output results, which contain NO 2 concentration for each given point in 3D space. The output is then transformed into x*y columns with z NO 2 values, that is, for each point in 2Dspace, there is a list of NO 2 concentration at different heights. Such data is then processed and analyzed using the SCAnalsysis object from EasyVVUQ.
The goal of the analysis is to provide us with the mean concentration (and associated uncertainty) for the whole domain at different heights, and to study how the final result may vary due to the incomplete emissions data. While at present the analysis is performed for a given point in 2D space for different heights above street level, future analyses will concern the entire 3D space.
The uncertainty quantification of the UrbanAir workflow for an arbitrary point in 2D-space is shown in Figure 13. Since the NO 2 concentrations is attributed to road transportation, it tends to decrease with increasing height above road level. Note that the interpolation of NO 2 concentration in between every 2 m is here due only to the plotting software. The standard deviation indicates how much uncertainty of input parameters (currently only four are taken into account) is reflected in the air quality predictions. In the forthcoming work, a sensitivity analysis will be conducted to select the most important uncertainty input parameters.

Discussion and Conclusions
In this work, we have applied EasyVVUQ to five diverse application areas, in order to extract information on sensitivity or Figure 13. Emissions of NO 2 at different heights above street-level from road transportation, with the mean (red line) and standard deviation (blue region) calculated using the EasyVVUQ campaign. uncertainty in these pre-existing models, without the need for intrusive modifications to the code. EasyVVUQ provides the tools necessary for computational scientists to add state of the art VVUQ algorithms to their simulation workflows without modifying the underlying codebase.
The library is intentionally execution-method agnostic, providing the base VVUQ workflow elements to allow for different execution patterns (such as Pilot Jobs) facilitated by any choice of middleware solutions. The agnosticism to choice of middleware (including using no middleware at all), and restartability of the workflow, provide the flexibility necessary for EasyVVUQ to be applied to many workflows in the HPC domain. For example, the Fusion application above uses the PSNC Pilot Job Manager to manage job execution, whereas the Ocean Circulation and Migration applications rely on FabSim3. Execution of the materials application, meanwhile, is handled manually by the user. Other middleware solutions may be used, such as RADICAL Cybertools, [9] Dask JobQueue, [46] or cloud submission tools.
The encoding and decoding steps of a standard EasyVVUQ script ensure that application-specific information is abstracted from the rest of the VVUQ workflow. This keeps the UQ algorithms in the sampling elements entirely generic. As such, multiple sampling elements may be chained or combined into more complex sampling elements (such as via use of the MultiSampler element). Complex encoding may also be achieved through combining multiple encoders into a single MultiEncoder element.
This generic approach is intended to accommodate switching to different UQ methods at no development cost to the user, allowing users to easily try out a variety of UQ approaches. It is intended that many more UQ algorithms will be integrated into this framework over time.