A nonintrusive adaptive reduced order modeling approach for a molten salt reactor system

We use a novel nonintrusive adaptive Reduced Order Modeling method to build a reduced model for a molten salt reactor system. Our approach is based on Proper Orthogonal Decomposition combined with locally adaptive sparse grids. Our reduced model captures the effect of 27 model parameters on k eff of the system and the spatial distribution of the neutron ﬂux and salt temperature. The reduced model was tested on 1000 random points. The maximum error in multiplication factor was found to be less than 50 pcm and the maximum L 2 error in the ﬂux and temperature were less than 1%. Using 472 snapshots, the reduced model was able to simulate any point within the deﬁned range faster than the high-ﬁdelity model by a factor of 5 (cid:1) 10 6 . We then employ the reduced model for uncertainty and sensitivity analysis of the selected parameters on k eff and the maximum temperature of the system. (cid:1) 2020 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
Complex systems such as molten salt reactors impose a modeling challenge because of the interaction between multi-physics phenomena (radiation transport, fluid dynamics and heat transfer). Such complex interaction is captured with high-fidelity, coupled models. However, these models are computationally expensive for applications of uncertainty quantification, design optimization, and control, where many repeated evaluations of the model are needed. Reduced Order Modeling (ROM) is an effective tool for such applications. This technique is based on recasting the high fidelity, high dimensional model into a simpler, low dimensional model that captures the prominent dynamics of the system with a controlled level of accuracy. Many ROM approaches can be found in literature (Antoulas et al., 2001). However, amongst studied ROM methods, Proper Orthogonal Decomposition (POD) is the suitable method for parametrized, nonlinear systems (Benner et al., 2015). The POD approach is divided into two main phases: the first is the offline phase, where the reduced order model is constructed by solving the high fidelity model at several points in parameter space to obtain a reduced basis space; the second is the online phase, in which the reduced model is used to replace the high fidelity model in solving the system at any desired point with a reduced computational burden.
POD can be implemented intrusively by projecting the reduced basis onto the system's governing equations or non-intrusively by building a surrogate model for the POD coefficients. Many studies have successfully implemented projection based POD for nuclear applications (Buchan et al., 2013;Sartori et al., 2014;Lorenzi et al., 2017;Manthey et al., 2019;German and Ragusa, 2019).
However, for practical nuclear reactor applications, the intrusive approach is often challenging because these models are usually implemented with legacy codes that prohibit access to the governing equations, or built with coupled codes that renders modifying the governing equations a complicated task. In this case, a nonintrusive approach can be adopted to build a surrogate model for the coefficients of the POD basis. Simple interpolation or splines can be used (Ly and Tran, 2001) or for high-dimensional problems, Radial Basis Function (RBF) is usually employed (Buljak, 2011). Neural networks (Hesthaven and Ubbiali, 2018) and Gaussian regression (Nguyen and Peraire, 2016) have also been studied to build the surrogate model. These approaches rely on standard sampling schemes (Monte Carlo, Latin Hypercube Sampling, tensorized uniform) to generate the snapshots. Such strategies do not take into account the dynamics of the problem and can be expensive

Contents lists available at ScienceDirect
Annals of Nuclear Energy j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / a n u c e n e for problems parametrized on high-dimensional spaces. Audouze et al. (2009) suggested tackling this issue by combining the POD-RBF method with a greedy residual search. In this approach, the residual of the PDE is used as an error estimator by iteratively placing sampling points at locations that minimize the residual until a certain global criterion is achieved. However, this method requires repeated evaluations of the residual, which can be expensive in some solvers (e.g. matrix-free solvers) or unavailable for legacy solvers.
In this work, we propose the use of ROM method that combines the nonintrusive POD approach with the sparse grids technique (Bungartz and Griebel, 2004) to build a reduced model of a fastspectrum molten salt system. Our approach is implemented using a previously developed algorithm (Alsayyari et al., 2019) that uses locally adaptive sparse grids as a sampling strategy for selecting the POD snapshots efficiently. The adaptivity is completely nonintrusive to the governing equations. In addition, the algorithm provides a criterion to terminate the iterations, which can be used as a heuristic estimation for the error in the developed reduced model. In this work, we extend the algorithm to deal with multiple fields of outputs. In addition, we demonstrate how local derivatives can be computed for local sensitivity analysis. The liquid-fueled system under investigation is a simplified system that captures the main characteristics of the Molten Salt Fast Reactor (Allibert et al., 2016). An in-house multi-physics tool (Tiberga et al., 2019), coupling an S N radiation transport code with an incompressible Navier-Stokes solver, was considered as the reference model of the molten salt system. We use the developed adaptive-POD (aPOD) algorithm to construct a ROM for this reference model. We then employ the built reduced model for an uncertainty and sensitivity analysis application to study the effect of the parameters on the maximum temperature and the multiplication factor. The uncertainty and sensitivity analysis was accomplished with extensive random sampling of the reduced model. Such approach is only achievable due to the efficiency provided by the reduced model over the reference model.
The remainder of this paper is organized as follows: the POD method is briefly introduced in Section 2. Section 3 presents the sparse grids approach by introducing the interpolation technique first followed by the method for selecting the sampling points. The aPOD algorithm along with the approach to deal with multiple fields of outputs and computing the local derivatives are presented in Section 4. The model for the molten salt system is given in Section 5. The discussion of the results of constructing the reduced model is in Section 6. The uncertainty and sensitivity analysis is in Section 7. Finally, conclusions are presented in Section 8.

Proper orthogonal decomposition
In a nonintrusive manner, the Proper Orthogonal Decomposition can build a ROM by considering the reference, high fidelity model as a black box mapping a given input to the desired output. Let the reference model f ðy; xÞ be dependent on state y and a vector of input parameters x. We can then find an expansion approximating the model as follows: where c i is the expansion coefficients which depends on the input parameter x and u i ðyÞ is the corresponding basis function. The POD method seeks to find the optimal basis functions u i ðyÞ that minimizes the error in L 2 norm, The basis functions are chosen such that they are orthonormal. Thus, the coefficients c i ðxÞ can be computed as c i ðxÞ ¼< f ðy; xÞ; u i ðyÞ >; where < f ðxÞ; gðxÞ >¼ R f ðxÞgðxÞdx. Assuming that the reference model is discretized (f ðy; xÞ ! fðxÞ), the POD snapshot method finds the solution to the minimization problem using the Singular Value Decomposition (SVD). This approach begins with sampling the reference model at discrete points in parameter space ½x 1 ; x 2 ; . . . ; x p , where p is the number of sampling points. Then, the corresponding outputs ½fðx 1 Þ; . . . ; fðx p Þ can be arranged in a matrix M called the snapshot matrix. Finally, we obtain the basis vectors (also called POD modes) u i as the first r left singular vectors of the SVD on the matrix M, where r is chosen to be less than or equal to the rank of the matrix M. A truncation error can be quantified using the singular values of the SVD ðrÞ if r is chosen to be strictly less than the rank of M, where n is the rank of M. e tr quantifies the error in approximating the solutions contained in the snapshot matrix.

Sparse grids
For an accurate POD reduced model, the snapshots need to cover the entire dynamics of the reference model within the defined range of input parameters. Therefore, selecting an effective sampling strategy is crucial for the success of the reduced model. We propose an algorithm that is based on locally adaptive sparse grids to select the sampling points. The sparse grid algorithm builds a surrogate model for each of the POD coefficients using a Smolyak interpolant. Iteratively, the algorithm identifies a set of important points and samples their neighbouring points in the next iteration (Griebel, 1998). This process is repeated until a global convergence criterion is met. In this section we introduce the methods for the interpolation and the selection of the sampling points.

Interpolation
The Smolyak interpolation is a hierarchical interpolant that can be implemented in an iterative manner such that the accuracy is increased with each iteration (Barthelmann et al., 2000). Different basis functions can be used for the interpolant. We choose piecewise linear functions with equidistant anchor nodes since they are suitable for local adaptivity. The equidistant anchor nodes, x i j , corresponding to level i are defined as (Klimke, 2006) The unidimensional nodes from Eq. (6) can be shown in a tree structure (Fig. 1) where the depth of the tree is assigned a level index i. The algorithm is iterative where at each iteration k, it defines a set of important points Z k . The criterion for selecting the important points is presented in Section 3.2. Once Z k is identified, the interpolant at iteration k for a function (cðxÞ) depending on a d-dimensional input x can be given by with A 0;d ðcÞðxÞ ¼ 0, where m D k is the cardinality of Z k , and H n is the d-variate basis function for the point x n 2 Z k , where x n has support nodes ðx i 1 n;1 ; . . . ; x i d n;d Þ, and i p is the level (tree depth) index for the support node x ip n;p . w k n is called the surplus which is defined as The union of the important points from all iterations up to k are collected in the set Because of the tree structure arrangement of the points, each point in the sparse grid (x ¼ ðx 1 ; . . . ; x d Þ) has ancestry and descendant points. All the descendant points fall within the support of the basis function anchored at that point. The first generation descendants of a point are neighbouring points called forward points. The forward points for n points in the set S ¼ fx q jq ¼ 1; . . . ; ng are defined with an operator WðSÞ as follows: where bðxÞ is a function that returns the parent of a node x from the tree. Likewise, the first generation ancestor points are called backward points and defined with an operator W À1 ðSÞ as follows: Finally, an operator CðSÞ that return all ancestors for the points in S can be defined as 3.2. Selecting the important points The algorithm builds the reduced model in an iterative fashion. At each iteration, we generate a set of trial points to test the model. The model is then updated according to the results of this test. Let the generated trial points be stored in the set T k , where k is the iteration number. The method for generating the trial points will be discussed in Section 4. For any point x q 2 T k , we can define a local error measure k q in the L 2 -norm as follows: where r k is the number of POD modes selected at iteration k. The number of POD modes is selected such that the truncation error (Eq. 4) is below a defined tolerance c tr . Once k q is computed for all points in T k , we can select points with an error above a certain threshold to be stored as candidate points. The candidate points are defined as where c int is an interpolation threshold and f abs is the absolute tolerance, which is introduced to deal with functions of small magnitude. The candidate points indicate the regions in which the model needs to be enriched. To enrich the model, the ancestor points of these candidate points are first considered because ancestors have wider support. If all ancestors of the candidate points were considered important from previous iterations, that point is taken as important because the error at that point ( k q ) is above the desired threshold. This is formulated as follows: On the other hand, if a point x q in iteration k has an error k q above the threshold but has also an ancestor point y i which was not included in the important set in the previous iterations, x q will not be marked important but its ancestor y i will be marked important, because it is possible that the error k q was large due to missing the ancestor which has a wider support, that is Then, the complete set of important points at iteration k is formed by Eqs. (18) and (19)as

Algorithm
Points that are not included in the important set Z k are added to the inactive set I k to be tested in subsequent iterations. The trial set of the next iteration (k þ 1) is generated as where cardð:Þ is the cardinality operator, and l is a greediness parameter which has a value 2 ½0; 1. The trial set (T kþ1 ) is formed by the forward points of Z k . However, some of these forward points are excluded from being evaluated if they have some backward points not considered important in previous iterations. The number of excluded points is tuned with l. For l ¼ 1, all points are tested regardless of their ancestry (the algorithm in this case is more exploratory) whereas the algorithm is more efficient for l ¼ 0 by not testing points that have any backward points not included in The trial set (T kþ1 ) is then used to sample both the reduced model and the reference model to compute the error kþ1 q . Then, the important points (Z kþ1 ) are identified and added to the snapshot matrix. Each update to the snapshot matrix generates a complete new set of POD modes, which requires recomputing the interpolant A k;d ðcÞðxÞ because of its dependence on the POD modes. Specifically, the surpluses (w k q;h ) corresponding to POD mode u h need to be recomputed with each POD update. The surpluses are just the deviations of the interpolant from the true value. Therefore, an easy way to update the surpluses after each iteration is as follows: whereû g is the gth POD mode after updating the snapshot matrix, u h is the hth POD mode before updating the snapshot matrix, w k q;h is the surplus at iteration k corresponding to the point x q 2 X k and POD mode u h , andŵ k q;g is the updated surplus corresponding to x q 2 X k andû g . For further reading regarding the adaptive sparse grids technique and the derivation of Equation 22, see Alsayyari et al. (2019) and the references within. Fig. 2 summarizes the algorithm.

Multiple outputs
To deal with models of multiple outputs, we can build a different ROM model for each output, which entails running the adaptive-POD algorithm separately for each output. With such an approach, managing the output field data is important to prevent multiple costly evaluations of the same point. This can be achieved by storing all output fields for any full model evaluation in a data bank, which the algorithm is directed to access when a point is required more than once in different output field constructions. With this strategy, the separate runs of the algorithm are performed in series rather than parallel in order to avoid full evaluations of the same point. Another approach is to combine the output fields by stacking them into a composite vector which is then treated as a single output in the snapshot matrix. In this approach, only a single ROM is built to represent all outputs. Since the first approach is a straightforward application of the algorithm, in this section, we show how the second approach is implemented.
Let the outputs be represented by f 1 ðxÞ; . . . ; f o ðxÞ where o is the number of output fields. The snapshot matrix is formed by stacking the output fields as We can compute the local error measure (Eq. (16) Different interpolation thresholds and absolute tolerances can be defined for each output. A point x q is admitted to the candidate set (Eq. (17)) if the corresponding error k s;q at any of the output fields (s 2 ½1; . . . ; o) is greater than the defined threshold where c int;s and f abs;s are respectively the interpolation threshold and the absolute tolerance defined for output f s ðxÞ.
The algorithm is terminated when a global criterion is met. We define this criterion to be where f rel;s is the global relative tolerance set for output f s ðxÞ. Note that the multiple-outputs approach can yield a different performance compared to the single-output approach in terms of points selected for evaluations. This is because the POD basis is constructed differently. In the single-output approach, the POD modes are tailored to that output specifically whereas in the multipleoutputs approach the POD modes contain information for all output fields.

Calculation of local sensitivities
To compute local sensitivities, we can find an analytical expression for the derivatives of each output with respect to the inputs. The derivative of the ROM model in Eq. (1) with respect to the gth dimension x g is The ROM model interpolates c i ðxÞ with the operator A k;d ðcÞðxÞ. Using Eqs. (8) and (9), Eq. (27) where the derivative of the unidimensional basis function @ @x a i (dropping the dependence on the dimension g) is computed as @ @x a 1 It is evident that due to the choice of piecewise linear basis functions, our reduced model is non-differentiable at the anchor nodes x i j , which implies that we cannot compute local derivatives at the sampled snapshots, including the nominal point. However, we can compute the local derivatives at two points very close to the nominal values and average them out to have a measure of the local sensitivities at the nominal point.

Molten salt system
In this work, we construct a reduced order model of a simplified system representative of the main characteristics of the Molten Salt Fast Reactor (Allibert et al., 2016): strong coupling between neutronics and thermal-hydraulics, fast spectrum, and transport of precursors. The problem was developed as a benchmark for multi-physics tools dedicated to liquid-fuel fast reactors (Aufiero and Rubiolo, 2018;A. Laureau et al., 2015). Fig. 3 depicts the problem domain: a 2 m side square, 2-dimensional cavity filled with fluoride molten salt at initial temperature of 900 K. The cavity is surrounded by vacuum and insulated; salt cooling is simulated via a heat sink equal to hðT ext À TÞ, where T ext ¼ 900 K and h is a volumetric heat transfer coefficient. Zero-velocity boundary conditions are applied to all walls except the top lid, which moves at v lid ¼ 0:5 ms À1 . The steady-state solution is sought with criticality eigenvalue calculations normalizing the reactor power to P 0 . Fluid properties are con-stant with temperature and uniform in space. Neutronics data are condensed into 6 energy groups and temperature corrected only via density feedback, to avoid the complexities related to Doppler feedback modeling; delayed neutron precursors are divided into 8 families. The flow is laminar and buoyancy effects are modeled via the Boussinesq approximation. Cross sections are corrected according to where T ref ¼ 900 K and qðT ref Þ is the density at which macroscopic cross sections are provided. They correspond to the reference values chosen for the Boussinesq approximation. b th is the thermal expansion coefficient. We refer to Aufiero and Rubiolo (2018) and A. Laureau et al. (2015) for a more detailed description of the problem. An in-house multi-physics tool is used to model the molten salt system. It couples a solver for the incompressible Navier-Stokes equations (DGFlows) with a neutronics code solving the multigroup S N Boltzmann equation coupled with the transport equations for the delayed neutron precursors ðPHANTOM À S N ). Both codes are based on the Discontinuous Galerkin Finite Element method for space discretization. Fig. 4 displays the structure of the multiphysics tool and the data exchanged between the codes. The average temperature on each element (T avg ) is outputted to PHANTOM À S N , which applies the density feedback on cross sections taken from the library at 900 K, according to Equation 32. Then, the neutronics problem is solved taking the velocity field (u) from DGFlows as another input for the delayed neutron precursors equation. Finally, the fission power density (P fiss ) is transferred to the CFD code. The steady state solution is sought by iterating DGFlows and PHANTOM À S N until convergence. More details on the multi-physics tool can be found in Tiberga et al. (2019).
Simulations of the molten salt system were performed choosing a 50 Â 50 uniform structured mesh, with a second-order polynomial discretization for the velocity and a first-order one for all the other quantities. An S 2 discretization was chosen for the angular variable. Fig. 5 shows the steady state fields (velocity magnitude, temperature, and total flux) obtained for the nominal values of the input parameters. The nominal multiplication factor in this configuration is k eff ¼ 0:99295. The upper bounds for each of the six energy groups are shown in Table 1 along with the space averaged flux (U avg ) for each group in the nominal case.

Results
A ROM model was built for the molten salt system by considering 27 input parameters. We assumed a uniform distribution for all paramters. The parameters and the corresponding percentage variation from the nominal values are summarized in Table 2, where P 0 is the initial power, b th is the thermal expansion coefficient, R f ;g is the fission cross section for group g; b i is the delayed neutron fraction for precursors family i; k i is the decay constant for precursors Fig. 3. Simplified molten salt fast system: square cavity domain. It is insulated, surrounded by vacuum, and filled with molten fluoride salt at initial temperature of 900 K. The top lid moves with velocity v lid ¼ 0:5 ms À1 . Fig. 4. Computational scheme of the multi-physics tool representing the highfidelity model. The CFD code, DGFlows, exchanges data with the radiation transport code, PHANTOM À SN, at each iteration due to the coupling between the physics characterizing the molten salt nuclear system. family i, v lid is the lid velocity, m is the viscosity, and h is the heat transfer coefficient. Since we aim at using the reduced model for uncertainty and sensitivity analysis, we assigned a variation of AE10% for parameters with typical experimental uncertainties whereas we vary design parameters (P 0 , v lid and h) by AE20%. Our interest is in the effect of these parameters on the spatial distribution of the total flux UðrÞ, the temperature TðrÞ, and the value of the effective multiplication factor k eff . Therefore, the reference model has 27 inputs and returns a value for the k eff and two field vectors each of length 7500 corresponding to the coefficients of the discontinuous Galerkin expansion for the total flux U and temperature T. In this work, we compare the stacking of the outputs approach described in Subsection 4.1 with the single-output approach. For the multiple-outputs approach, the snapshot matrix for the outputs evaluated at points ½x 1 ; . . . ; x p is computed as The global relative tolerances f rel for U and T were set to be 10 À2 , which means we require the error in the L 2 norm for these fields to be less than 1%. For k eff , we require the error to be less than 50 pcm, so we set f rel for k eff to be 50 Â 10 À5 . The interpolation threshold (c int ) was chosen to be one order of magnitude less than the set relative tolerances. Therefore, c int was 10 À3 for both U and T and was set to be 5 Â 10 À5 for k eff .
We first built a reduced model using a greediness value l ¼ 1.
For the multiple-outputs approach, the algorithm required 4495 reference model evaluation to converge. However, only 142 points were included in the important set. The small number of selected important points is an indication of oversampling. The algorithm was then run again with l ¼ 0. In this case, the algorithm sampled 472 points with 105 important points included in the snapshot matrix, which is a reduction by about a factor of 10 in the number of evaluations compared with the case of l ¼ 1. Each reference model evaluation takes about 1.5 h to run (performed on a Linux cluster using 1 CPU operating at 2.60 GHz). Therefore, this reduction in number of evaluations is massive in computational time.
In order to test the model, 1000 Latin Hypercube Sampling (LHS) points were generated. LHS is a method to generate unbiased random points in higher dimensional spaces by partitioning the hypercube first. Then, drawing one sample from each partition. These generated points were not part of the snapshot matrix. Note that the reduced model was trained only on the important set. The rest of the model evaluations served as trial points but were not included in the snapshot matrix. In machine learning terminology, the important set is the training set and the rest of the evaluations served the function of the validation set (Ripley, 1996). Therefore, the generated 1000 unbiased random points in the test set represent 10 times more testing points than training points. Running the reduced model on the 1000 testing points needed only about one second on a personal computer. Table 3 summarizes the maximum L 2 norm error found for each output. It is evident that all tested points resulted in errors well below the set tolerances. We also compare the results of the single-output approach to the multiple-outputs approach in the same table. While both approaches satisfied the required tolerances, the number of full model evaluations required in the offline stage was different. The single-output approach required fewer evaluations compared to the multiple-outputs approach. This is due to the fact that the POD modes in the single-output approach are tailored to that output field. The algorithm in this case, samples points to construct a specific reduced model satisfying the desired tolerance for that output. In the multiple-outputs approach, on the other hand, the algorithm uses POD modes containing information for all output fields, which require more points to satisfy the desired tolerances for every output fields. However, because the reduced model is enriched with every additional sampling point, the multiple-outputs model has a slightly less error in the online phase compared to the single-output approach. Fig. 6 shows the distribution of the L 2 norm error for the tested 1000 random points for each output in the reduced model of the multiple-outputs approach and l ¼ 0. A comparison between the temperature distributions of the reduced model and the reference full order model at the point that resulted in the maximum error is shown in Fig. 7. The L 2 norm error for this case was 0.2% while the maximum absolute difference locally was 13.9 K, which is about 1% of the maximum local temperature (about 1482.6 K). Both cases of l ¼ 1 and l ¼ 0 converged with 3 iterations (k = 3). To highlight the cost effectiveness of the adaptive approach, for such 27dimensional problem, the classical (non-adaptive) sparse grid approach would require 27829 points after 3 iterations, which is extremely expensive to run. Table 4 summarizes the number of unique nodes per dimension, which was found to be the same for both the single and multiple-outputs approaches. This number is indicative of the linearity/non-linearity of the reference model. During the construction stage, the algorithm captures the degree of linearity of the output of the reference model with respect to each dimension within the defined range. A value of 3 means that the algorithm considered that dimension to be constant because after building a constant interpolant at the root 0:5, the error in the model was found to be within the defined tolerances at the children points {0; 1}. The algorithm then stopped further refinements along that dimension. A value of 5 indicates that the model is piecewise linear in the segments (0,0.5) and (0.5,1) with respect to that dimension because the refinement is stopped after testing the piecewise linear interpolant using the first 3 points {0.5,0,1} at the children {0.25,0.75}. A value higher than 5 indicates that the model is nonlinear along that dimension.
It is evident from the number of unique nodes that the algorithm found the outputs of the model to be constant (within the set tolerances) with respect to b i and k i , which means varying these parameters within the 10% range does not significantly affect the defined outputs. Additionally, the model was found to be piecewise  linear with respect to the power, velocity, thermal expansion coefficient, viscosity, and the fission cross section for the groups 1-4. However, for the lowest energy groups (group 5 and 6), the model was nonlinear. This can be explained by the fact that the flux distributions for all groups were not changing significantly due to the homogeneity of the changes to the system. In addition, the group fluxes were found to have the same order of magnitude as shown in Table 1 for the nominal case. However, the nominal values of the fission cross section for R f ;5 and R f ;6 are higher compared to the other fast groups, which weigh more in the calculation of k eff . By examining the cause for the additional unique points along R f ;5 and R f ;6 , we found that they were triggered purely by k eff and not by U or T. The model was also nonlinear in the heat transfer coefficient. The negligible effect of b i and k i explains the reason for the massive reduction in number of evaluations with the setting l ¼ 0. The algorithm in this case recognized that b i and k i have    Note the change of the colour bar scale in the difference plot (right).
no effect within the defined range and stopped sampling points along these dimensions. Since b i and k i amount to 16 out of the 27 dimensions, the reduction in number of points was massive.

Uncertainty quantification and sensitivity analysis
In this section, we demonstrate the potential of the built ROM model in an application of uncertainty quantification and sensitivity analysis. We study the effect of the selected input parameters on the maximum temperature and the multiplication factor k eff . The resulting ROM can be sampled cheaply at any point within the specified range. The ROM model from the multiple-outputs approach and l ¼ 0 is employed for the study in this section. However, we do not expect differences in the results if any of the other 3 ROM models developed in Section 6 were used instead. We use Latin Hypercube Sampling to sample the reduced model with 100,000 random points. The density histograms approximating the Probability Distribution Function (PDF) are shown in Fig. 8. For comparison, the densities resulting from running the reference model on the 1000 testing points are also shown in the figure. The density histogram shows a distribution close to a normal distribution, which can be explained by the fact that all input parameters are assumed to have uniform distribution and the model is linear or almost linear in these parameters. Therefore, the sum of these uniform distribution approaches the normal distribution. The normal probability plot in Fig. 9 confirms that the distribution is normal within the middle range while the deviation from the normal is seen at the tails of the distribution. The mean of the maximum temperature was found to be at 1336.5 K with standard deviation equal to 61.1 K while the mean of k eff was 0.99229 with standard deviation equal to 0.016.
Local and global sensitivity analyses were also performed using the built ROM. For the local sensitivities, Table 5 presents the averaged derivatives computed from several points within a distance of 10 À14 (measured in the unit hypercube ½0; 1 d ) from the input's nominal values. In order to provide a better comparison of the effect of the parameters, the computed derivatives in the table are normalized by the ratio R 0 =x p;0 , where R 0 is the desired response (maximum temperature or k eff ) computed at the nominal values of the input parameters x p;0 .
The results show that the maximum temperature is mainly affected by the initial power P 0 and the heat transfer coefficient h. This is expected because these two parameters directly control the amount of energy present in the system. Higher initial power increases the amount of energy in the system which directly raises the temperature. The heat transfer coefficient, on the other hand, is negatively correlated with T max because lower h decreases the amount of energy being extracted from the system causing the temperature to rise.
The thermal expansion coefficient is related to the natural convection phenomenon. Forced and natural convection play a competing role in terms of mixing of the salt in the cavity. There are two vortexes in the cavity as shown by the streamlines in Figure 5 (left) for the nominal case. When forced convection increases, the larger vortex grows causing the vortex centre to move towards the cavity centre. In this case, salt in the central region of the cavity would always circulate around the centre where the fission power is maximum. On the other hand, when natural convection increases, the smaller vortex in the bottom left corner becomes larger causing the salt to pass through the centre then transported close to the boundaries of the cavity where the thermal energy is minimum. Hence, in the range of variations considered in this work, natural convection tends to redistribute the heat in the cavity, whereas forced convection has the opposite effect. Higher b th causes natural convection to be more prevalent over forced convection. This causes the temperature to be more uniform. For this reason, b th is negatively correlated with T max . The viscosity, on the other hand, has the opposite effect. Increasing the viscosity reduces the mixing of the liquid, which creates more concentrated hot spots that increase the maximum temperature. The lid velocity is also positively correlated with the maximum temperature Table 4 Number of unique nodes per dimension.
Parameter number of unique nodes Parameter number of unique nodes Fig. 8. Density histograms of the maximum temperature (left) and the multiplication factor k eff (right) by sampling the reduced model with 100,000 points. The distributions of same variables from sampling the reference model with the 1000 testing points are also shown. Note that the histogram is normalized such that the sum of the areas of the bars equals to 1.. because it increases the forced convection. However, this correlation is shown to be weak because the range in which the velocity changes (AE20%) is very small. The fission cross sections have negligible effect on T max . The delayed neutron fractions and the precursors decay constants have zero derivatives because our reduced model assumes them to be constants at any point. The multiplication factor is mainly affected by the fission cross sections as expected. The fission cross sections of the two lowest energy groups are the most important. This is because of their higher weight (higher nominal values compared to the fast groups with similar flux magnitudes) in computing k eff . The thermal expansion coefficient is negatively correlated with k eff because by increasing b th , the liquid is mixed more, which in turn causes more precursors to move from regions of higher importance to regions of lower importance near the boundaries. The initial power is negatively correlated with k eff due to the negative temperature feedback coefficient of the system. For the same reason, the heat transfer coefficient is positively correlated with k eff . The lid velocity and viscosity have negligible effect on the multiplication factor.
For the global sensitivities, we computed the first order Sobol indices using quasi Monte Carlo method with Sobol sequence sampling (Sobol, 2001). We selected the size of our sampling matrices to be 10 5 , which generates 2 matrices each of dimension 10 5 Â 27. The first order Sobol indices were then estimated using the estimators recommended by Saltelli et al. (2010). The computed indices for both T max and k eff are shown in Fig. 10. The Sobol indices show agreement with the conclusions of the local sensitivities. The maximum temperature is predominantly sensitive to P 0 and h while b th and m have a slight effect on T max . The multiplication factor, on the other hand, is mainly sensitive to the fission cross sections with the Fig. 9. Normal probability plots for the maximum temperature (left) and the multiplication factor k eff (right) showing the distribution to be normal within the middle parts but deviating from the normal distribution at the tales.  lowest energy groups having the most importance. Although R f ;5 has a nominal value of about half R f ;6 , the Sobol index of R f ;5 is about 50% higher than R f ;6 . This can be explained by the higher flux magnitude of group 5 compared to group 6 as can be seen from the average flux value reported in Table 1 for the nominal case. P 0 and h have a reduced effect while b th has a minimal effect on k eff . The agreement between the local and global sensitivities show that the system is only weakly nonlinear. Additionally, the sum of the computed first order Sobol indices was found to be very close to one, which indicates that second and higher order interactions between the parameters are almost negligible. This confirms the weak nonlinearity of the model.
In total, 3 Â 10 6 model evaluations were performed to complete the uncertainty and sensitivity analysis study. The time to perform these simulations using the reduced model was about 45 min on a personal computer, which is about half the time to perform a single simulation of the full model on the computer cluster. Using 472 snapshots computed in the offline phase, we obtained a gain of about a factor of 5 Â 10 6 in the online computations with respect to the reference model. This demonstrates the advantage of ROM for such applications.

Conclusions
The developed ROM algorithm (aPOD) based on POD and the adaptive sparse grids method was applied to a coupled model of a test case for the Molten Salt Fast Reactor. We selected 27 input parameters to model their effect on the distribution of the flux and temperature, and the value of the multiplication factor. In a completely nonintrusive manner, aPOD was able to build a representative (1% accurate) ROM model with 4495 model evaluations. This number was effectively reduced by a factor of 10 with the setting l ¼ 0. This great reduction was successfully achieved due to the ability of the algorithm to automatically recognize that the 16 dimensions corresponding to b i and k i have no significant effect within the defined range. It was also observed that the initial power, thermal expansion coefficient, fission cross section of the fast 4 groups, lid velocity, and viscosity all have piecewise linear effect on the outputs. On the other hand, the fission cross section for the 2 lowest energy groups and the heat transfer coefficient have slight nonlinear effect. As a test of the model, 1000 Latin Hypercube Sampling points were tested and compared with respect to the reference model. The errors were found to be well within the defined tolerances for all outputs. The multipleoutputs approach was found to require more sampling points to satisfy the desired tolerances compared to a single separate run for each output. This can be explained by the fact that with the single-output ROM model, the POD modes are tailored to that output field and the algorithm only needs to sample points to satisfy the tolerance for that field. The multiple-outputs approach requires the composite POD modes to represent all output fields, which leads to more sampling points to satisfy the tolerances. However, because of the additional sampling in the construction of the reduced model, the error was found to be lower for the multiple approach compared to the single approach.
For an application of uncertainty and sensitivity analysis, we studied the effect of the 27 input parameters on the maximum temperature and the multiplication factor. The density histograms showed a normal distributions of these variables, which can be explained by the uniform distribution assumption of the selected parameters and the weak nonlinearity of the model with respect to the input parameters within the defined ranges. The maximum temperature was shown to be sensitive to the initial power and the heat transfer coefficient while the multiplication factor was mainly sensitive to the fission cross sections as expected. The uncertainty and sensitivity study was performed using a total of 3 million random points, which were completed in about half the time to run a single simulation of the reference model. The nonintrusive approach of the algorithm provides great potential for studies of complex coupled nuclear systems such as the molten salt reactor, particularity in applications of uncertainty quantification, sensitivity analysis, fuel management, design optimization, and control.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.