Analysis of the Molten Salt Fast Reactor using reduced-order models

In this paper, we present a reduced-order modeling approach to study the Molten Salt Fast Reactor (MSFR). Our approach is nonintrusive and based on the proper orthogonal decomposition method. We include adaptivity in selecting the sampling points both in time and parameter space. Steady-state and transient analysis were both performed using the developed models. In the steady-state analysis, we capture the effect of 30 model parameters on the spatial distributions of fission power and temperature, and on the multiplication factor. The dimensionality of the fission power was reduced from the 104288 nominal dimensions in the physical space to 10 dimensions in the reduced space, whereas the temperature was reduced from 220972 dimensions to 3. The reduced model was then used for uncertainty and sensitivity study of the maximum temperature in the reactor and the multiplication factor. In the transient analysis, the reduced model captured the effect of perturbations in the flow rate of salt in the intermediate circuit on the fission power density and temperature. The reduced models were successfully tested on a set of points that were not part of the snapshots used during the construction stage.


Introduction
Molten salt reactors have gained interest due to their potential safety, reliability, and sustainability (Generation IV International Forum, 2002).Different designs of this concept have been proposed in the literature Dolan (2017).In this work, we consider the Molten Salt Fast Reactor (MSFR) (Allibert et al., 2016).A key design feature of this unmoderated reactor is the use of a liquid salt fuel, which also plays the role of the coolant.This design introduces a unique modeling challenge because of the tightly coupled neutronics and thermal-hydraulics phenomena (e.g., transport of delayed neutron precursors, distributed thermal energy deposition directly in the coolant, a strong negative temperature feedback coefficient).To address these challenges, high-fidelity coupled models are used to provide an insight into the behavior of the reactor (e.g., Aufiero et al., 2014, Fiorina et al., 2014, Laureau et al., 2017, Cervi et al., 2019, Tiberga et al., 2020b).For safety assessment applications, an accurate and explicit quantification of the propagation of uncertainties through these complex models is required (International Atomic Energy Agency, 2019).Quantifying uncertainties and analyzing sensitivities in reactor physics can be accomplished using adjoint methods (Gilli et al., 2013;Perkó et al., 2013).However, adjoint-based methods require the availability of an adjoint solver, which might not always be feasible for coupled problems.Another approach is using forward-based methods, which requires repeated evaluation of the high-fidelity model for different parameter order to employ the projection.In cases where the high-fidelity model is a legacy coupled solver or closed-source, intrusive approaches are not applicable.Moreover, the stability of the high-fidelity model is not preserved in the reduced model with the use of projection-based POD methods (Amsallem and Farhat, 2012).
To overcome these issues, nonintrusive approaches have been proposed in the literature, which are techniques that can be wrapped around the high-fidelity solver and avoid the need to access the original set of equations.Nonintrusive spectral projection methods such as polynomial chaos expansion (PCE) is one approach for uncertainty quantification applications (Perkó et al., 2014a,b), which has been applied to a molten salt fueled system (Santanoceto et al., 2021b) and to perform a preliminary uncertainty and sensitivity analysis of the MSFR steady-state (Santanoceto et al., 2021a).However, such methods do not reduce the dimensionality of the high-fidelity model.Dynamic Mode Decomposition (DMD) is a nonintrusive model order reduction method that has been employed to model the MSFR (Ronco et al., 2020).DMD constructs a linear operator from a sequence of snapshots of the system.Then, a reduced model preserving the dynamics of the original system is formulated, which can be used for transient analysis and control applications (Tu et al., 2014).Propagating uncertainties in the initial conditions has also been studied using DMD (Takeishi et al., 2017;Abdo et al., 2019).However, for problems with parametrized models where uncertainty and sensitivity analysis for multiple model input parameters are to be studied, applying DMD is cumbersome.
An alternative ROM approach is through a nonintrusive POD implementation, where a surrogate model is trained to compute the POD coefficients instead of projecting the model equation onto the reduced basis.In nuclear reactor applications, a Range Finding Algorithm (RFA) has been used in Huang et al. (2017) to build the reducedbasis subspace(referred to as active subspace) combined with a simple polynomial surrogate for the POD coefficients.Hybrid approaches aiming to simultaneously reduce the dimensionality of both the state and the input variables using RFA have also been explored (Abdel-Khalik et al., 2013;Bang et al., 2015).We have proposed a nonintrusive adaptive POD algorithm in Alsayyari et al. (2019).A key difference in our developed algorithm is the adaptive sampling of the highfidelity model as opposed to random or a priori uniform sampling employed by the RFA.Our algorithm is based on adaptive sparse grids technique, which is suited for problems with higher dimensional input spaces.An uncertainty and sensitivity application of the algorithm on a simplified two-dimensional molten salt fueled system was presented in Alsayyari et al. (2020), where we have also addressed handling systems with multiple-outputs.We have presented an extension of the algorithm in Alsayyari et al. (2021) to deal with time-dependent parametrized problems.In this work, we investigate the MSFR reactor using the developed algorithm, which demonstrates the capability of the algorithm on a large-scale, three-dimensional MSFR model.We consider two applications.The first is a steady-state reduced model for uncertainty and sensitivity analysis of 30 model parameters.The second is a transient reduced model, which can be used for transient analysis and control applications.
The remainder of the paper is organized as follows: Section 2 presents an introduction to the problem formulation and the POD method.A summary of the adaptive sampling algorithm is presented in Section 3.Then, a description of the MSFR model is given in Section 4. The results for the steady-state analysis along with the uncertainty and sensitivity study are given in Section 5. We present the results of the transient study in Section 6.Finally, conclusions are discussed in Section 7.

Proper orthogonal decomposition
Consider a general high-fidelity model that is described by a function  ∶ R  → R  , which maps input parameter  ∈ R  to a state (field) vector  ∈ R  , where  is the size of the input space (representing, for example, material properties, and boundary conditions), and  is the size of the state space (representing, for example, space, angle, and energy).We consider  as a black-box model, which can either describe single physics or coupled multi-physics code, that is,  =  (). (1) The high-fidelity model  can be evaluated at any given   and produce an output   (i.e.,   =  (  )).We are interested in finding an approximation for   without the use of the computationally intensive  .For this task, we build a reduced model that produces an approximation φ , which takes the form of where   is the th basis vector in the reduced basis space,   () is the corresponding coefficient, which is a function of the parameter , and  is the size of the reduced basis space.The POD approach constructs a reduced basis space from snapshots of the state vectors  such that the error is minimized in a least square sense, that is, it solves the minimization problem The basis vectors are defined to be orthonormal (i.e., ⟨  ,  ℎ ⟩ =  ℎ , where ⟨., .⟩indicates the scalar product).The reduced basis space can be found using the method of snapshots.Let  be a matrix collecting snapshots of the solution  at some selected values of , that is, where   is the solution of the high-fidelity model ( ) at parameter value   , and  is the number of selected snapshots.Note that for transient problems where  computes a solution  that is time-dependent, we consider the time to be a parameter, rather than an independent variable -that is,  includes time as one of the parameters.Such a formulation allows for time adaptive sampling as described in Alsayyari et al. (2021).
The basis vectors are then found to be the left singular vectors of the singular value decomposition (SVD) applied to the matrix , that is, let  be decomposed using SVD as  =   ⊤ , then   ∀ 1 ≤  ≤  are the first  column vectors of  .The size of the reduced space is determined by , which can be chosen to truncate  such that the sum of the squared singular values () corresponding to the neglected singular vectors is below a predefined threshold  tr .
Once the basis vectors are known, we can use the orthogonality of the space to compute the coefficients at the sampled point (  ) as (6)

Adaptive sampling
Using the values of the coefficient   () at the sampled points, one can train a surrogate model to compute the solution of the system at any point.Different surrogate models can be employed for such a task.To deal with higher dimensional input spaces, we choose the Smolyak sparse grid interpolation technique.The hierarchical structure of the interpolant allows for the desired adaptive strategy, which can further reduce the burden of dealing with high dimensional spaces.This section is devoted to summarize the developed adaptive algorithm, which is detailed in Alsayyari et al. (2019).
Our implementation of the sparse grid technique uses localized adaptive sampling.Without loss of generality, the -dimensional space of the input parameter  is mapped to a unitary hypercube [0, 1]  .The sparse grid technique first generates unidimensional nodes along each dimension.Then, points in the parameter space are formed by a specific combination of the generated unidimensional nodes.Different choices for the unidimensional node generation rule are possible.However, we choose the equidistant rule to increase the separation between points.The unidimensional nodes are arranged in a tree structure as shown in Fig. 1.Each level is assigned an index  and contains several nodes.Nodes are added at each level at half the distances between the nodes from the previous levels.Each node is connected to two children at the next level ( + 1) and one father from the previous level ( − 1).The root node is considered a father for itself.There is, however, an exception at level  = 2, where each node has one child because these nodes mark the boundaries of the unit hypercube.
Each node in the tree can be uniquely identified with the level index  and an index .For each node from the tree    along dimension , we can build a basis function      () as where the dependence on the dimension  is dropped for notational convenience.The level parameter   is defined as follows: A point  in parameter space is formed by combining nodes from all dimensions (i.e.,  = ( )). Extending the tree structure to points in parameter space allows us to define forward points.The first forward points for a point  along the first dimension are (( ), where () is a function that returns the children of the node .In general, the forward points along the th dimension are ( Because each node has at most two children, each point has at most 2 forward points.We can also define backward points in the same manner by applying a function  * (), which returns the father of the node.Therefore, each point has at most  backward points.Recursively generating the forward points creates the set of descendant points.On the other hand, the ancestor points for a point  are formed by recursively generating backward points until the root point (0.5, … , 0.5) is reached.
The adaptive algorithm is iterative.In every iteration, a subset of points is selected and marked important from a set of trial points.Let the iteration number be  and a set that collects the selected important points at each iteration  be   .The union of the important points sets (  ) from all iterations up to  forms the set of selected grid points   , that is, A point   is marked as important if it satisfies one of two conditions.The first is to have an approximation error   above a defined threshold  int with all ancestors of that points included in   .The second is to have a descendant point with an error above  int .Points that do not meet any of the conditions are added to a set of inactive points.The error is computed in the  2 norm as where   is the solution returned by the high-fidelity code ( ) at point   , φ  is the approximation produced by the reduced model at the point   and iteration , and  is an offset introduced for cases when the norm of the solution is near-zero.
At iteration , a surrogate model for the POD coefficient   () is built as where  , ()() is the Smolyak interpolation operator applied to () that depends on the iteration  and dimension of the input space .
For the initialization ( = 0), we enforce  0, ()() to be zero.The term  , ()() is defined as where    is the cardinality of the set   .The function   () is the -variate basis function for the point The surplus    is defined as the difference between the interpolated value and the true value of the coefficient at   , At each iteration, a trial set is generated from the forward points of the important set.The model is tested on each point within the trial set.To provide an estimate for the global error in the reduced model at iteration  (  ), an appropriate norm can be chosen over the error defined in Eq. ( 10).In this work, we have chosen the maximum norm such as the reduced model error is defined by where   is the cardinality of the trial set at iteration .However,   only provides an estimate of the error in the reduced model and is not rigorous.

Algorithm
In the initialization step ( = 0), we select the root point at the center of the hypercube (0.5, 0.5, … , 0.5) and evaluate the high-fidelity model  .That point is added to the important set  0 , and the snapshot is added to the matrix .Then, set  = 1 and do 1.Generate the forward points of the points in  −1 ; 2. Add points that have all their backward points in  −1 (defined as in Eq. ( 9)) to the trial set; 3. Evaluate the high-fidelity model  at the points in the trial set; 4. Compute the coefficients   () at the points in the trial set using  −1, ()(); 5. Compute the error at every point in the trial set using Eq.(10); 6.If the maximum found error (  ) was below a prescribed global tolerance , terminate; Otherwise, continue; 7. Find the important points   ; 8. Add the snapshots corresponding to the important points in   to the matrix ; 9. Perform SVD on the matrix  to extract the reduced basis; 10.Use Eqs. ( 11)-( 14) to construct  , ()(); 11.Set  ←  + 1 and go to Step 1.
For a detailed description of the algorithm, see Alsayyari et al. (2019).

MSFR model
This section presents the MSFR main design features and modeling approach.A more detailed description of the model can be found in Tiberga et al. (2020b) and the references within.
A schematic illustration of the fuel circuit is shown in Fig. 2. The reactor core is a toroidal cavity where the liquid fuel salt (a mixture of lithium, thorium, and fissile nuclides fluorides) can flow freely without any moderator or control rod.Sixteen identical sectors branch out from the central cavity.Each sector contains a pump, a heat exchanger, and a unit for the separation and treatment of the helium gas dispersed in the fuel salt to control reactivity and remove metallic fission products (Delpech et al., 2009).A blanket with fertile salt surrounds the cavity while reflectors are placed at the top and bottom of the core.The heat exchangers transfer thermal energy to the intermediate circuit filled with inert molten salt, which in turn, delivers the heat to the energy conversion system consisting of a conventional Joule-Brayton cycle.Table 1 summarizes the main MSFR design parameters.
The reactor was modeled using an in-house multi-physics tool, which couples an incompressible Reynolds-averaged Navier-Stokes solver (DGFlows) (Tiberga et al., 2020a;Hennink et al., 2021) with an   multi-group neutronics code (PHANTOM-S  ) (Kópházi and Lathouwers, 2012).The latter is equipped with transport equations for the delayed neutron precursors to model their movement.Both codes employ the discontinuous Galerkin finite element method for space discretization and implicit backward differentiation formulae (BDF) time schemes.Fig. 3 displays the structure of the multi-physics tool and the data exchanged between the codes.PHANTOM-S  receives the average temperature on each element (  ) and corrects the cross sections accordingly with respect to a library at the reference temperature  0 .The corrections take into account the effects of density feedback (which has a linear dependence on temperature) and Doppler feedback (which has a logarithmic dependence).The velocity and turbulent viscosity fields ( and   ) are also taken from DGFlows as another input to solve the delayed neutron precursors equation.Then, the fission power density (   ) is transferred to DGFlows as it constitutes the right-hand side of the energy equation.The steady-state solution is sought by iterating the solvers until convergence (four iterations are typically sufficient).On the other hand, transient simulations are performed with a loose-coupling strategy, first computing a time-step with DGFlows and then calling PHANTOM-S  to solve the neutronics problem.We refer the reader to Tiberga et al. (2019Tiberga et al. ( , 2020b) ) for a more comprehensive description of the coupled code.
Because of the symmetry in the reactor design, we modeled only one recirculation loop.The geometry considered is reported in Fig. 4.While the neutron flux is calculated in the reflectors and the blanket, the CFD code (DGFlows) neglects heat transfer in these regions.Fig. 5 shows Data are exchanged at each iteration between the two solvers to model the coupled physics phenomena characterizing the MSFR (Tiberga et al., 2020b).the meshes used in this study.Neutronics calculations were performed on an unstructured mesh consisting of 26072 tetrahedra (of which 21489 are in the fuel salt domain).This mesh is finer in the core region, where neutron importance is high, and coarser in the external sector.This master mesh was then refined once uniformly in the outer-core region to obtain the CFD mesh, which consists of 55243 elements.A second-order polynomial was used for the velocity discretization, while a first order polynomial was used for all other quantities.
The fuel salt composition is reported in Table 2, along with some physical properties that were fixed and selected for the uncertainty analysis study of Section 5.The neutron energy range was condensed into six groups, with boundaries shown in Table 3. Delayed neutron precursors were grouped into eight families.All neutronics data were evaluated at temperature  0 = 900 × 10 0 K with Serpent (Leppänen et al., 2015) using the JEFF-3.1.1 data library (Santamarina et al., 2009).
For the CFD calculations, symmetry boundary conditions were imposed at the wedge sides, and standard wall-functions with adiabatic conditions were assumed at all walls.For the neutronics calculations, reflective boundary conditions were assumed at the sides of the wedge, while vacuum conditions were imposed elsewhere.For the transport of the precursors within the fuel salt, homogeneous Neumann and no-inflow conditions were imposed at all walls.
Lacking detailed design specifications for the primary heat exchanger, salt cooling was modeled via a volumetric heat sink term equal to ℎ int,0 (  int,0 −  ) , where ℎ int,0 and  int,0 are the nominal volumetric heat transfer coefficient and average temperature of the intermediate salt whose values are reported in Table 1.The fuel pump was modeled with a momentum source term, and buoyancy was taken into account through the Boussinesq approximation (the reference density and thermal expansion coefficient are reported in Table 2).We considered the flow field to be fixed at the nominal state because natural convection has a negligible contribution to the total nominal flow rate, and the pump specifications were fixed throughout the analyses in

Table 2
Properties of the fuel salt mixture (Benes et al., 2013).Melting point [K] 854.15 × 10 0 this work.Moreover, since we are interested in reactor steady-state and operational transients, decay heat was not taken into account.For the transient calculations described in Section 6, the time-dependent equations were discretized with the second-order BDF scheme with a fixed time-step size of 0.1 × 10 0 s.

Steady-state analysis
In this section, we study the effect of selected parameters of interest on three outputs: fission power density ( fiss ), temperature distribution ( ), and the multiplication factor,  ef f .The fission power density  fiss is computed on the neutronics mesh with 104288 unknowns (nominal dimensionality) while the temperature  is computed on the CFD mesh with 220972 unknowns.This analysis is performed on the steady-state model of the MSFR with a fixed reactor power of 3000 MW.However, since we are modeling 1/16th of the reactor, the observed nominal power is 187.5 MW.The selected parameters to be investigated were specific heat capacity of the fuel salt (  , influencing the temperature distribution), heat conductivity of the fuel salt (), fission cross sections for the six energy groups (  , ), capture cross sections for the six energy groups ( , ), delayed neutron fractions for eight families of precursors (  ), and corresponding decay constants (  ).The nominal values of the selected parameters are reported in Table 4.These nominal values are theoretical as the actual salt properties and the associated uncertainties are still under investigation.However, for the propose of this work, we assume an uncertainty following a Gaussian distribution for all parameters with a mean () equal to the nominal value and a standard deviation () taken as 5% of the mean.

Construction of the reduced model
Our implementation of the adaptive approach was developed for bounded input domains.For this reason, we consider the Gaussian distribution to be truncated.Truncating the Gaussian distribution, in this case, is a valid approximation because the parameter uncertainties we are considering are epistemic.Therefore, we are not altering a random process but rather limiting the scope of analysis to the region with the highest probability of having the true value.Moreover, truncating the distribution prevents unphysical values of the parameter, such as negative .The truncation was selected to be at 3, retaining 99.7% of the probability range.This implies that the range of variation for all parameters is set to be ±15% of the nominal value.We first constructed the reduced model to be uniformly accurate within the defined range, then use the reduced model for the uncertainty and sensitivity analysis  employing the corresponding probability distribution of each parameter.A separate reduced model was built for each output.The global relative tolerance  was set to be 1% for  fiss , 0.1% for  , and 50 pcm for  ef f .The interpolation threshold  int was 1 × 10 −3  fiss , 1 × 10 −2 for  , and 5 × 10 −5 for  ef f .The POD truncation threshold  tr was 10 −12 for all outputs.
After construction, each reduced model is tested on 1000 independent points that were not part of the snapshots generated during the constructions.Latin Hypercube Sampling (LHS) was used to draw the random testing points from the input space.Table 5 summarizes the test results for each model.It can be seen that all models resulted in a maximum relative error that was below the set tolerance .These results certify that the reduced models are an accurate representation of the high-fidelity MSFR model with in the desired tolerances.Fig. 6 compares the reduced model for  fiss with the reference model at the point with maximum error.The difference is seen to be maximum in the central region of the reactor core, where the flux is maximum.Fig. 7 shows the comparison at the maximum error for  .The maximum difference for this case is observed at the bottom of the core, where the relative error locally is around 0.3%.A single high-fidelity computation requires about 4 CPU-hours (performed on a Linux cluster) whereas the reduced model produces the results in less than a second.
Table 5 reports also the number of POD modes after truncation, representing the number of dimensions in the reduced space.The dimensionality of  fiss was reduced from the 104288 nominal dimensions in the physical space to 10 dimensions in the reduced space whereas  has 3 dimensions in the reduced space compared to 220972 nominal dimensions in the high-fidelity model.The low number of POD modes for  reflects the fact that the temperature profile in this reactor is fairly constant with respect to the perturbations considered.This is due to the fact that the reactor power is fixed, and the perturbations are applied homogeneously.Note that  ef f is a single-valued response, and reducing the dimensionality is not applicable for this model.
The total number of evaluations requested by the algorithm during the construction stage is also reported in Table 5.This number corresponds to the total number of executing the high-fidelity model computations.The models for  fiss and  needed 61, and 63 evaluations respectively, while the  ef f model needed 1639 evaluations.The increased number of evaluations reflects the strong nonlinearity of  ef f with respect to the input parameters compared to  fiss and  .This finding is supported by the number of unique nodes per dimension for each output, which is reported in Table 6.This number is the result of projecting the final grid points in   onto each dimension.It is a measure of the nonlinearity in the output with respect to the dimension as captured by the reduced model.A value of 1 indicates that the reduced model considered that dimension to have a negligible effect on the output.A value of 3 means that the output is linear or piecewise linear with respect to that input parameter because the reduced model used 3 nodes (the root node and the first two children) to construct a linear or a piecewise linear interpolant and no further refinement was required along that dimension.A value of 5 or more indicates that the interpolant along that dimension constructed a nonlinear interpolant between the nodes.In general, the nonlinearity in the interpolant is proportional to this number.Thus, Table 6 also explains the seemingly high difference between the number of evaluations required for the  ef f compared to  fiss and  , simply caused by the combination of the nonlinear dependence of  ef f on many of the input parameters simultaneously and the curse of dimensionality, leading to many more potential combinations of parameters the algorithm needs to check.
Table 6 shows that the fission power density was found to be linear with respect to the specific heat capacity   , which is explained by the temperature feedback effect.The fission cross sections of groups 2 to 6 also have a linear effect on the fission power density, which is expected since the fission power is proportional to   , .The fission cross section of the first group, however, was observed to have a negligible effect within the tolerances on  fiss .This can be explained by the lower magnitude of the first-group flux compared to the rest of the groups, as shown in the average flux value per group in Table 3.This has an effect of a lower weight in the calculation of the fission power density.The capture cross section of the most thermal group is the only capture cross section that has a linear effect on  fiss .This is due to its larger nominal value resulting in a larger range of variations compared to the other groups.The rest of the parameters had a negligible effect on  fiss and were considered as constants.The temperature is shown to be nonlinear with respect to   having 5 unique nodes and unaffected by the rest of the parameters within the defined tolerances.One can deduce that the 5 important points selected by the algorithm to be in   and reported in Table 5 correspond to the 5 points obtained by varying   only.The multiplication factor is nonlinear with respect to the cross section of groups 3 to 6 both fission and capture while having a piecewise linear interpolant with respect to the two most energetic groups.The specific heat capacity is also seen to have a linear effect on  eff through the temperature feedback.The delayed neutron fractions of families 2, 4, and 5 were the only families with a significant effect on  eff .This, however, is due to the larger nominal value of these parameters, which results in a larger range of variations compared to the rest of the families.The decay constant is shown to be taken as a constant for all families with 1 node each, which indicates that within the set tolerance of 50 pcm and a range of variation of ±15%,   had no effect on  eff .

Propagating uncertainties
We used the reduced models to propagate uncertainties in the parameters to the responses of interest.We consider two model responses; The first is the maximum temperature of the system, which has a value in the nominal state of 1084.8 × 10 0 K, and the second is the

Table 6
Number of unique nodes per dimension for each output.A value of 1 indicates that the output is constant with respect to the parameter.A value of 3 signals that the output is piecewise linear in the parameter.A higher value indicates that output is nonlinear in the parameter.
multiplication factor with a nominal value of 1.00999.To extract the probability density function (PDF) of the response, we run the reduced model at randomly sampled points drawn from the distribution of the input parameters.A histogram of the response is an approximation of the PDF.Fig. 8 shows the normalized density histograms for the maximum temperature and multiplication factor by using 100,000 random points.
The maximum temperature has a mean of 1085 × 10 0 K with a standard deviation of 6.9 × 10 0 K.The distribution is close to a normal distribution with a slight skew to lower temperatures.Fig. 9 shows the normality plot of the data, which is a measure of the degree of deviation of the data from the normal distribution.The normal probability plot shows that this skew is mostly observed in lower probabilities of temperatures between 1100 K to 1108 K.The tails deviate from the normal distribution due to the truncation of the normal distribution in the input parameters.The deviation from the normal distribution confirms the nonlinearity of the temperature with respect to   , as observed by the algorithm in Table 6.The multiplication factor has a mean of 1.01009, with a standard deviation of 0.01899.The distribution of  eff follows a perfect normal distribution, which is confirmed by the normal probability plot in Fig. 9.
In order to study the sensitivity of the response to the input parameters, the first order Sobol indices are computed from the reduced models.First order Sobol indices provide a measure for the partial contribution of a single parameter to the response of interest.However, from the number of unique nodes per dimension given in Table 6, the reduced model for temperature was shown to be sensitive only to  one parameter (  ).Therefore, we only compute Sobol indices for  eff .We used a quasi-random sampling Sobol sequence (Sobol, 2001).The sequence was generated on the unit hypercube, then mapped to the distribution of each parameter using inverse sampling of the cumulative distribution function (CDF).Two sets of sampling points each of size 10 5 were used to compute the Sobol indices as given by Saltelli et al. (2010).The Sobol indices in Fig. 10 show  eff to be sensitive only to the cross sections.The indices ranked the fission and capture cross sections of group 5 to be the most significant, followed by groups 3 and 4. The fast groups (groups 1 and 2) and the most thermal groups have   relatively lower impact on  eff .This can be explained by the fact that the reactor operates with an epithermal spectrum.This is confirmed from the average flux value per group in Table 3.The indices also show that, in general, fission cross sections have higher importance on  eff compared to the capture cross sections.This is expected because fission cross sections have a direct impact on the multiplication factor while the capture cross section impacts  eff through the loss term which also includes the system leakages.The sum of the first order Sobol indices was found to be close to 1, which suggests that the effect of higher-order interactions between parameters can be neglected.

Transient analysis
In this section, we consider the transient model of the MSFR where the reactor power is no longer fixed.The reactor is kept at steady-state by dividing the fission operator by nominal value of  ef f .Therefore, the initial conditions for the transients are the nominal values.We are interested in a reduced model for control and simulation purposes, capturing the dynamics of the fission power density and temperature with respect to perturbations in the salt flow rate of the intermediate circuit.This case is simulating an operational scenario where the reactor power is controlled through adjustments in the flow rate of the salt in the intermediate circuit.The intermediate circuit extracting heat from the heat exchangers was not modeled explicitly.However, for the transient analysis, we simulate the effect of controlling the reactor )) , (18) where  int () and ℎ int () are respectively the intermediate salt temperature and heat transfer coefficient as functions of time ,  int,0 is the nominal value of the intermediate salt temperature with a value of 908.15 × 10 0 K (as reported in Table 1), and ℎ int,0 is the nominal value of the heat transfer coefficient with a value of 19.95 × 10 0 MW∕m 3 K (Table 1).
Because we considered time as an input parameter in our formulation, this model has two input parameters, the change in the intermediate circuit flow rate () and time ().The outputs are the fission power density, which is computed on the neutronics mesh with 104288 unknowns, and the temperature, which is computed on the CFD mesh with 220972 unknowns.We considered the flow rate to range between −40% to +15% of the nominal values.The time range was taken  ∈ [0, 180] s -that is, we are interested in constructing a reduced model for the first 180 × 10 0 s after a perturbation in the intermediate flow rate .
A transient of increasing the flow rate by 15% is shown in Fig. 11, along with snapshots of the final solutions at  = 180 × 10 0 s.Increasing the flow rate of the intermediate circuit causes the temperature of salt flowing back into the cavity from the heat exchanger to drop, which is observed as a decrease of the minimum and average temperatures.The lower temperature in the reactor core introduces positive reactivity because of the strong negative temperature feedback of this reactor (Heuer et al., 2014).For this reason, the trends show an increase in the reactor power and the maximum temperature registered at the top  of the reactor cavity.The average temperature then gradually adjust the downward trend, and a new steady-state is reached.
To build the reduced model, we set the tolerances  = 1% and  int = 0.1% for the fission power density, and  = 0.1% and  int = 0.01% for the temperature.The POD truncation tolerance  tr was set to be 1 × 10 −12 for both models.A summary of the results is given in Table 7.The algorithm required 9 simulations (corresponding to 9 values of ) for the fission power density and selected a total of 78 snapshots to build the reduced model.For the temperature, the algorithm required 5 simulations and selected 62 snapshots.In the transient case, we define a snapshot to be the model output taken at a certain time and fixed parameter configuration whereas a simulation   to an a priori snapshot selection approach, which does not consider the actual dynamics of the system.
To test the constructed reduced models, we ran 8 high-fidelity test simulations with values of  at half the distances between the points selected during the construction phase.In this manner, we maximize the distance (along  dimension) between the test simulations and the simulations used for the construction of the reduced models.We then select 1000 random snapshots at different times from the testing simulations, which are used to test the reduced model.A histogram of the relative  2 norm error is shown in Fig. 13 for both models.The error in the temperature model was found to be well below the set tolerance, with the maximum error being 0.06%.For the fission power density model, most of the points resulted in an error below the set tolerance.However, 6 out of the 1000 points were above the tolerance, with the maximum being 1.8%.The point that resulted in the maximum error is compared with the high-fidelity solution in Fig. 14 for the fission power density and Fig. 15 for the temperature.In this case, the fission power density comparison shows the maximum difference to be at the center of the cavity, where the flux is maximum.The local relative error of at that point is 1.8%.The maximum temperature difference, on the other hand, is observed at the inlet of the reactor core where the temperature is minimum, with a local relative error of 0.1%.A full simulation to  = 180 × 10 0 s is completed with the reduced model in under 5 s, while the high-fidelity model needed 150 CPU hours on a high performance computing unit.

Conclusions
We have applied an adaptive POD approach to a large-scale, threedimensional model of the MSFR.The developed algorithm was able to construct reduced models for both steady-state and transient analysis of the reactor.The steady-state analysis considered higher dimensional input space with 30 parameters.Three reduced models were constructed to capture the effects of those parameters on the fission power density and temperature distributions, and on the multiplication factor.Each model was tested on 1000 random points that were not part of the snapshot selections.The model for fission power density required 61 F. Alsayyari et al. high-fidelity evaluations and resulted in a maximum relative  2 error of 0.24%.The model for the temperature required 63 evaluations and resulted in a maximum error of 0.02%.The multiplication factor model needed 1639 evaluations and the maximum error was 37 pcm.Evaluating the reduced models was completed in less than a second while the high-fidelity model needed 4 CPU-hours.
The reduced models were then used to propagate uncertainties in the input parameters to the maximum temperature of the reactor and the multiplication factor.The constructed PDF showed the multiplication factor to follow a normal distribution with a mean 1.01009 and standard deviation of 0.01899.The maximum temperature had a mean of 1085 × 10 0 K, a standard deviation of 6.9 × 10 0 K, and showed a distribution close to normal with a slight skew to lower temperatures.The sensitivity study concluded that the maximum temperature of the MSFR is only sensitive to the specific heat capacity within the defined tolerances.Heat conductivity and neutronics data had no significant effect on the temperature.Therefore, since the salt properties of the MSFR are still under investigation, we recommend prioritizing studies to reduce uncertainties in the specific heat capacity over heat conductivity.The multiplication factor was shown to be only sensitive to the cross sections.A standard deviation of 5% in the decay constant and delayed neutron fraction is sufficient to characterize the multiplication factor within 50 pcm error.However, a standard deviation of 5% in the cross sections resulted in a distribution of the multiplication factor with a standard deviation of 1899 pcm, which is significant.Therefore, uncertainties in cross section data should be below 5% standard deviation for this reactor.
The transient analysis studied the effect of controlling the flow of salt in the intermediate circuit on the fission power density and temperature distributions.Our approach considers time to be a parameter in input space in order to allow for an adaptive selection of the snapshots.The model for the fission power density needed 9 simulations and selected 78 snapshots.The test on 1000 independent points showed the maximum relative  2 error to be 1.8%.The model for the temperature required 5 simulations and selected 62 snapshots.The maximum error was found to be 0.06% for this model.The selected points showed that the algorithm sampled the high-fidelity model in regions of the beginning of the transient more than the steady-state.This allowed the algorithm to be efficient by simulating some of the points only to half the transient.A full simulation to the end of the transient required about 150 CPU-hours.Therefore, reducing the transient time for a point is a massive saving in computational resources.The constructed reduced models were able to produce a full simulation in less than 5 s.As follow-up work, this reduced model can be used to design a controller for the reactor.Other future applications could also include transients with significant primary circuit flow field changes (such as loss of pump power or pipe blockage events), where the temperature, neutronics and flow fields simultaneously have to be reduced.Since our algorithm is fully nonintrusive and can handle multiple fields, addressing such transients is straightforward and we expect our model to show comparable performance.
In all models, the number of points included in the final grid set   was a small fraction of the total snapshots generated by the algorithm.This fraction was about 50% for the transient models while for the steady-state models, it was found to be as low as 10%.This is an indication that the algorithm is still oversampling.The number of POD modes was found to be even a smaller fraction of the snapshots, which is a signal that most of the points were generated to train the surrogate model rather than discover new dynamics of the system.Therefore, the algorithm could be improved further with more advanced surrogate models to reduce the number of sampling points.The use of higherorder basis functions is a potential area of study for this propose.Moreover, our approach to the transient problem assumes the initial conditions to be parametrized.The approach is not yet applicable for problems with a time-dependent input signal, which cannot be parametrized.These areas of research are the subjects of future work.

Fig. 1 .
Fig.1.Illustration of the first 4 levels of the tree structure, where 0.5 is the root of the tree and nodes are added at half the distances between the previous nodes.Each node has 2 children except the nodes at level 2 where each has one child only.

Fig. 2 .
Fig. 2. A schematic illustration of the main design features of the MSFR fuel circuit.

Fig. 3 .
Fig. 3. Computational scheme of the multi-physics tool constituting the high-fidelity MSFR model.DGFlows is the CFD code, while PHANTOM-S  is the neutronics code.

Fig. 4 .
Fig. 4. Geometry of the MSFR recirculation loop used for simulations, showing the main regions considered in the model (Tiberga et al., 2020b).

Fig. 5 .
Fig.5.Mesh adopted for the MSFR model.The neutronics mesh (left), corresponding to the master-mesh, consists of 26072 tetrahedra (21489 in the fuel salt domain), while the CFD mesh (right) has 55243 elements.The latter is derived by refining the former once uniformly in the outer-core region.

F
.Alsayyari et al.

Fig. 6 .
Fig. 6.Fission power density distribution at the point of maximum error showing the reference model (top left), the ROM model (top right), and the distribution of the absolute difference (bottom).The relative  2 error was 0.24%.

Fig. 7 .
Fig. 7. Temperature distribution at the point of maximum error showing the reference model (top left), the ROM model (top right), and the distribution of the absolute difference (bottom).The relative  2 error was 0.02%.

Fig. 8 .
Fig. 8. Normalized density histograms for the maximum temperature (left) and the multiplication factor  ef f (right).The data was generated by sampling the reduced model with 100,000 random points drawn from the distribution of the input parameters.

F
.Alsayyari et al.

Fig. 9 .
Fig. 9. Normal probability plots for the maximum temperature (left) and the multiplication factor  ef f (right).

Fig. 10 .
Fig. 10.First order Sobol indices showing the first order sensitivities of  ef f to each input parameter.The sum of the first order Sobol indices is close to 1.
power through the salt flow rate in the intermediate circuit by adjusting the average intermediate salt temperature and heat transfer coefficient.Since we are interested in operational conditions, we employ a simple linear empirical model relating changes in the flow rate of the salt in the intermediate circuit to changes in the average intermediate salt temperature and heat transfer coefficients as int = −0.375, (16)ℎ int = 74825,(17)whereint is the change in intermediate salt temperature in Kelvin, ℎ int is the change in the heat transfer coefficients in units of MW∕m 3 K, and  is the percentage change in the flow rate of the salt in the intermediate circuit.In order to approximate the dynamics of controlling the salt flow rate, these perturbations are introduced in the model exponentially with a time constant of 10 × 10 0 s.Therefore, the models for the average intermediate salt temperature and heat transfer coefficient are  int () =  int,0 +  int (

Fig. 11 .
Fig. 11.A selected transient resulting from  = 15% showing the reactor power (top left), fission power density distribution at  = 180 × 10 0 s (top right), temperature trends for the maximum temperature  max , minimum temperature  min , and average temperature  avg (bottom left), and temperature distribution at  = 180 × 10 0 s (bottom right).

Fig. 12 .
Fig. 12.The generated points for the transient model of the fission power density (left) and temperature (right).The important points included in the final   are marked with a red circle.

Fig. 13 .
Fig. 13.Histogram of the relative  2 error in the fission power density model (left) and temperature model (right).A close up on the region above 1% is shown for the histogram of the fission power density model.

Fig. 14 .
Fig. 14.Fission power density distribution at the point of maximum error showing the reference model (top left), the ROM model (top right), and the distribution of the absolute difference (bottom).This point correspond to  = −36.5 and at time  = 2.4 × 10 0 s.The relative  2 error was 0.24%.

Fig. 15 .
Fig. 15.Fission power density distribution at the point of maximum error showing the reference model (top left), the ROM model (top right), and the distribution of the absolute difference (bottom).This point correspond to  = −36.5 and at time  = 1.5 × 10 0 s.The relative  2 error was 0.06%.

Table 3
Upper energy bound for each group and average flux of each group for the nominal case.

Table 5
Results for the steady-state reduced models for each output showing the maximum relative  2 error after a test on 1000 independent points, the number of POD modes after truncation, the number of evaluations to construct each model, and the final number of points in   .

Table 7
Summary of results for the transient reduced models corresponding to the fission power density  fiss and temperature  showing the number of required simulations to construct each model, number of snapshots selected, the final number of points in   , the number of POD modes after truncation, and the maximum relative  2 error after a test on 1000 random points.