Application of multiphysics model order reduction to doppler/neutronic feedback

In this paper, a proper orthogonal decomposition based reduced-order model is presented for parametrized multiphysics computations. Our application physics is Doppler feedback in a simplified model of the molten salt fast reactor concept. The reduced model is created using the method of snapshots where the offline training set is obtained by exercising a full-order model created with the OpenFOAM based multiphysics solver, GeN-Foam. The steady state models solve the multi-group diffusion k-eigenvalue equations with moving precursors together with the energy equation. A fixed velocity field is assumed throughout the computations, hence the momentum and continuity equations are not solved. The discrete empirical interpolation method is used for the efficient coupling of the ROM solvers, while the input parameter space is surveyed using the improved distributed latin hypercube sampling algorithm.


Introduction
Molten salt reactor (MSR) designs were originally developed in the mid-1950s at Oak Ridge National Laboratory (ORNL, USA) [1][2][3]. In MSRs, the nuclear fuel is in liquid form, dissolved in a salt. Salt compositions vary, but are typically based on fluorides or chlorides and include one or more of the following compounds: LiF, NaF, BeF 2 , ZrF 4 , KF, NaCl, MgCl 2 . Diverse variations on that reactor concept were investigated in the 1960s and 1970s, including graphite-moderated thermal-spectrum reactors at ORNL [4,5] as well as fast-spectrum burner reactors at Argonne National Laboratory [6,7]. In the early 1970s, MSR research was in competition for U.S. federal funding with sodium-cooled fast reactor systems and the MSR research in the USA dwindled down to lowpriority, low-funding efforts over the following decades, while the thermal-spectrum light water-cooled pressurized and boiling water reactors and the sodium-cooled fast reactors became the reference baseline worldwide for thermal and fast spectrum systems, respectively. Nonetheless, electricity-production priorities and safety requirements in the nuclear sector have evolved over the last 50 years and MSRs are now one of the six concepts selected for further investigation in the frame of the Generation-4 International Forum [8]. MSRs have highly promising features in terms of sustainability and safety features. Indeed liquidfueled MSRs can be designed to have strong negative-only * e-mail: jean.ragusa@tamu.edu reactivity feedbacks; they operate at atmospheric pressure; they allow for an online removal of gaseous fission products; and they give the possibility to drain the fuel salt in passively cooled and critically-safe tanks in case of emergency. Most fast-spectrum MSRs currently under development are based on pumped-loop designs, where the fuel salt is pumped outside of the primary vessel and transfers heat to a secondary coolant in separate heat exchangers. Examples of fast-spectrum molten salt reactor designs include: the MOlten Salt Actinide Recycler and Transforming (MOSART) project [9], the Molten Salt Fast Reactor (MSFR) concept based on fluoride salt, developed in the EVOL (Evaluation and Viability of Liquid Fuel Fast Reactors) [10][11][12] and then SAMO-FAR (Safety Assessment of the Molten Salt Fast Reactor) [13] programs under the auspices of EURATOM, and the Molten Chloride Fast Reactor (MCFR), currently developed by Terrapower [14]. Loop-type fast-spectrum molten salt reactors present new modeling challenges: -the fuel is in liquid form, yielding a more complex level of multiphysics coupling than traditional light water reactors (e.g., velocity fields needed to assess space/time location of fuel and delayed neutron precursors); -turbulent fuel-salt flow, leading to a large impact of turbulence modeling (thermal flow mixing in the core, effects of nozzle inlets,... [15]); -presence of gas bubbles in the salt, leading to compressibility and reactivity effects; -high solidification temperature of the salt, requiring the modeling of solidification/melting zones (wall solidification, frozen drain functionality, see [15]); -high operating temperature in the salt, necessitating the study of thermal radiation heat transfer.
The modeling challenges in fast-spectrum MSRs largely prohibit the use of simulation packages developed and tailored specifically for Light Water Reactors (LWRs). For instance, the Virtual Environment for Reactor Applications (VERA), developed as part of the Consortium for Advanced Simulation of Light water reactors (CASL) is being leveraged for MSR modeling but is only at an early stage of development for MSR and mostly focused on thermal-spectrum reactors, where salt flows mostly uni-directionally in graphite channels [16,17]. Hence, highfidelity models, based on first-principle physics, become the only available tool to gain understanding of the significance of foreseeable phenomena unique to molten salt and circulating fuel systems. Driven by the need for high-fidelity computational fluid dynamics (CFD), many research teams developing fast-spectrum MSRs rely on the open-source OpenFOAM CFD platform, either directly [18,19] or embedded in reactor physics packages such as GeN-Foam (Generalized Nuclear Foam) [20,21]. However, the simulation of such complex systems requires the solution of coupled partial differential equations, which can be computationally expensive to obtain. Model-order reduction comprises a set of empirical and mathematical techniques that can be used for lowering the computational complexity of the Full-Order Models (FOM) by creating Reduced-Order Models (ROM) that sacrifice a modest amount of accuracy for large gains in compute time. An effective mathematical technique for model-order reduction is the Reduced Basis (RB) method, which is based on the assumption that the solution of a complex model lives in a relatively small subspace. When parametric studies are to be performed, one should expand that subspace to account for solution variability. One manner by which such a subspace is "discovered" is through multiple full-order solutions, with adequate sampling of the parameter (input) space. This is known as the method of snapshots. The information contained in these snapshots is then extracted by means of Proper Orthogonal Decomposition (POD) [22,23] (e.g., via correlation matrix or singular value decompositions). The FOM is then Galerkin-projected on the obtained basis functions to yield the ROM. In this work, a POD-based ROM is created for parametrized multiphysics computations on the Doppler feedback effect in the Molten Salt Fast Reactor (MSFR) [24,25]. In this technique, the FOM is projected onto an appropriate subspace obtained via a POD of solutions determined by exercising the FOM itself. In other words, snapshots are taken at different states of the FOM and a POD is used to build a suitable basis for projection-based reduction. In the nuclear engineering community, applications of POD-based ROMs can be noted in reactor kinetics problems [26,27], fixed source, steady state neutral particle transport applications [28]. An approach for POD-ROM of eigenvalue problems connected to reactor physics has already been developed in [29][30][31][32][33].
However, the utilization of POD-based ROM for multiphysics applications can be challenging because the operators in the full-order models often do not have an affine decomposition, meaning that the ROMs involve costly, full-order operations, possibly yielding negligible savings in computation time. A method for the reduced-order modeling of multiphysics problems in nuclear engineering has also been derived in [30]. The approach described in the present work serves as an alternative. Instead of assuming an overall linear temperature dependence for the cross sections, here we opt for a hyper-reduction technique, namely, the Discrete Empirical Interpolation Method (DEIM) [34]. This allows the handling of an arbitrary, non-linear temperature dependence of the cross sections (in our case, in the logarithm of temperature). Our approach is exemplified using a liquid nuclear fuel system, more specifically, for a simplified 2D model of the Molten Salt Fast Reactor. The methodology has been integrated into GeN-Foam. For different applications of DEIM in other fields of engineering, we refer the reader to [35][36][37].
The rest of the paper is organized as follows. In Section 2, projection-based model order reduction is reviewed, first for linear systems, then for nonlinear systems, such as the ones encountered in multiphysics applications. In Section 3, a full-order simulation model is provided for a fast-spectrum molten salt reactor and its associated reduced-order model is derived. Section 4 describes the chosen MSFR computational model and its parametric input space. Results comparing the FOM and ROM models are provided in Section 5. We conclude and propose future work in Section 6.

Model order reduction: background
In this section, we review some of the basics of projectionbased model-order reduction. For brevity of the exposition, the reduced basis will always be sought through a POD decomposition and the reduced system will always be obtained using Galerkin projection. We refer the reader to [38,39] for other reduced basis approaches and to [40] for Petrov-Galerkin projection-based ROM.

Model order reduction for linear systems
First, we consider a linear parameterized steady-state FOM. After discretization, the FOM can be written as where A(µ) ∈ R N ×N is the discretized linear operator (a possibly large system of dimension N ), x(µ) ∈ R N is the solution vector and the parametric dependence of the model is denoted by the d-dimensional input parameter vector µ = [µ 1 , . . . , µ d ] T . In order to discover the lower dimension manifold in which the parametric solutions can be adequately represented, the FOM is exercised for a certain number of realizations of the input parameter vector where only the r most dominant modes are retained. A computationally effective manner to do this consists in generating the correlation matrix Q of the snapshots as Then, the eigenvalue decomposition of the correlation matrix is obtained where W = [w 1 , . . . , w Ns ] is the matrix of eigenvectors and Λ is a diagonal matrix containing the Λ i eigenvalues. The orthonormal basis vectors in V can then be constructed using the snapshots by Once this offline stage has been completed, the reduced order model is obtained by projection of the full order system: with the reduced linear operator and the reduced state vector x r (µ) ∈ R n . This linear system of size r N , is solved for x r , which are the coefficients of the full-order solution in the basis V . Hence, the full-order solution is reconstructed as For additional information, we refer the reader to [41,42].

Model order reduction for nonlinear systems
Next, we consider a nonlinear full-order model (FOM) which, after discretization, can be written as where a nonlinear vector-valued function has been added, Henceforth, for the sake of easier readability, the dependence of the operators on µ is only assumed and not shown explicitly. If one carries out the procedures described in Section 2.1, the resulting nonlinear reduced-order model is Even though this equation is expressed in terms of the vector of reduced unknowns, x r , solving it requires evaluating the full-order nonlinear vector-valued function, of size N r. This can be computational very expensive, often resulting in limited CPU time savings when solving the reduced order model, compared to the original full-order model. To remedy this, the Discrete Empirical Interpolation Method (DEIM) [34] is used in order to approximate F in a low-dimensional space by sampling it at only m N components. Hence, during the offline stage, a second matrix of snapshots is collected for the nonlinear function values, Similarly to the method described in equations (2)-(4), the POD of S F is computed and m of the basis functions for the nonlinear vector-valued function are retained to subse- and the resulting nonlinear reduced-order model system is Several remarks are in order: be pre-computed once for all; -P T F extracts only a few (actually m, with m N ) components of F that need to be evaluated. This is a large savings afforded by the use of the DEIM. Because discretization of partial differential equations usually results in local connectivity between an unknown and its neighbors, evaluating P T F (V x(µ)) is therefore computationally inexpensive when m N .
For additional information on DEIM, we refer the reader to [34,43].

Governing laws
In this section, we select mathematical models (partial differential equations) that will be used as FOM and derive their reduced-order counterparts. These models are said to be parametric in the sense that certain parameters are uncertain in the governing equations and will be sampled from appropriate probability density functions. In this work, the following parameter are assumed uncertain: the total reactor power (P th ), the heat exchange coefficient in the heat sink, i.e., in the heat exchanger (α ext ), the ultimate heat sink temperature (T ext ), and the volumetric area of the heat sink (A V ). The vector containing all the uncertain parameters will be denoted by µ = (α ext , A V , T ext , P th ) T . The mathematical models are selected to be representative of a simplified molten salt reactor design. The molten salt reactor application is discussed in a later section.

Full-Order Model (FOM)
In order to represent the physics of a MSFR, we define a multiphysics model comprising of -neutronics: six-group k-eigenvalue diffusion theory, with delayed neutron precursor balance equation with a drift (advection) term, and -thermal-hydraulics: here, we only consider an energy balance equation, with a given (but spatial varying) flow field computed using nominal values of the uncertain parameters. For the types of parametric modifications employed here, the effect of flow perturbations should be small. For the reduced-order modeling for turbulent flows, we refer the reader to [44,45].
The cross sections in the neutronics models are temperature-dependent. The multi-group diffusion keigenvalue problem [46] can be described as a coupled partial differential equation system expressing the balance of neutrons in each energy bin as where G denotes the number of energy groups, I the number of delayed neutron precursor groups, Φ g is the neutron scalar flux, D g is the diffusion coefficient, Σ r,g is the macroscopic removal cross-section and νΣ f,g is the total fission neutron yield times the macroscopic fission cross-section in energy group g. Furthermore, Σ g →g is the macroscopic scattering cross-section from group g to g, χ p,g and χ d,g describe the fraction of prompt and delayed neutrons released in group g, while β is the effective delayed neutron yield. Moreover, λ z denotes the decay constant corresponding to the delayed neutron precursor concentration C z of group precursor group z.
The problem is supplemented with boundary conditions: reflective boundary conditions ( ∇Φ g · n = 0, g = 1, . . . , G) are applied on the symmetry planes and zero-incoming current boundaries are applied on the other boundaries (−D g ∇Φ g · n = 1 2 Φ g , g = 1, . . . , G). These equations are coupled with the steady-state balance equations for the delayed neutron precursors, commonly expressed as where u is a precomputed (stationary) velocity field, β z is the delayed neutron yield for precursor group z and α eff is an effective diffusion coefficient that takes into account the effects of turbulent mixing as well. For simplicity, zero gradient boundary conditions ( ∇C z · n = 0 for z = 1, . . . , I) are used for each of the precursor equations on every wall. We stress that u is not uniform, but taken from a CFD simulation using the nominal values of the uncertain parameters (the RANS continuity and momentum equations were solved using k− model for turbulence [21]). Thus, u and α eff are known fields during the generation of the reduced operators. We leave model-order reduction of turbulent hydrodynamics for subsequent work. The multigroup neutron fluxes are normalized using the group-wise power cross section Σ p,g to ensure a certain total reactor power, P th : The cross sections are assumed to be temperaturedependent. For fast spectrum reactors, as shown in [47], a logarithmic interpolation between pre-computed databases (for example at 900 K and 1200 K) of the group constants yields good results: where T denotes the temperature field. Finally, to be able to account for the temperature feedback, the following energy balance equation is solved: where ρ is the density, c p is the heat capacity, k T the effective thermal conductivity of the fluid, while α is the heat transfer coefficient and A V is the volumetric area of the heat sink. All of the parameters in the equation above are assumed to be constant. Furthermore, the heat transfer coefficient α is computed as the harmonic mean of two coefficients, one that characterizes the heat transfer between the salt and the structure of the heat exchanger (α salt ) and a second describing the heat flow between the heat exchanger and the external heat sink (α ext ). For simplicity, zero gradient boundary conditions ( ∇T · n = 0) are used for every surface in the model. The iteration scheme used for the solution of the coupled problem is discussed in details in [21]. α ext , A V and T ext are assumed uncertain; hence this model is parametric in those input parameters.

Reduced-Order Model (ROM)
The ROMs are constructed by physics-wise (equationwise) projection of the FOM equations onto suitable reduced bases obtained applying POD to the solution snapshots. During the offline phase, the solution fields are collected into corresponding snapshot matrices and a Proper Orthogonal Decomposition is carried out for each snapshot matrix separately. This segregated approach has been proven to be effective for k-eigenvalue multigroup problems [33]. The approximate solutions in the reduced P. German et al.,: EPJ Nuclear Sci. Technol. 5, 17 (2019) 5 spaces can be written as where ψ g,i is the ith basis vector of the subspace selected for the neutron flux in group g, C z,i is the ith basis vector for the precursor group z, τ i denotes the basis vector of temperature and L i is the ith basis vector of the logarithmic temperature. Moreover, f g , c d , t and l vectors contain the coordinates of the approximated flux in group g, approximated precursor concentration in group d, approximated temperature and logarithmic temperature within their corresponding reduced spaces. Using these, the ROMs for the multi-group diffusion equations can be described as where k r is the largest eigenvalue of the reduced system and the reduced operators are computed using the approximation mentioned in equation (15) together with a Galerkin projection onto the corresponding subspace. Thus, the elements of the reduced diffusion operator can be expressed as where · denotes the volumetric integral of the given scalar fields that can be carried out numerically. Even though it is not shown explicitly, this reduced operator takes into account the boundary conditions imposed on the scalar flux by incorporating a ψ g,i , 1 2 ψ g,j Γ term (in case of vacuum boundary) for every boundary face Γ of the computational domain. Again, it must be mentioned, that D g depends on the approximate logarithmic temperature (see Eq. (16)). By translating this linear dependence into the ROM, the expressions of the reduced operators become slightly more convoluted and can be written as where the values of D const g and D var g can be determined using equation (13) and the coefficients of the logarithmic temperature (l i ) are computed using the coefficients of the temperature (t i ) through the Discrete Empirical Interpolation Method, described in Section 2.2 in detail, as follows: It should be noted that, in this case, it is enough the carry out the DEIM in terms of the logarithm of the temperature due to the fact that the nonlinearity in the model is caused by the temperature-dependence of the cross sections, and the non-linear function involves the logarithm of the temperature. As noted in the previous Section, we recall that the reduced operator can be precomputed in a tensor form and every time the operator has to be reconstructed due to the changing temperature field, it only requires the summation of reduced matrices making the coupling of the reduced models extremely efficient. The same treatment is applied to the other reduced operators containing cross-sections, even though in the following definitions it is not shown explicitly. The additional reduced operators in equation (17) are computed as One can notice that these reduced matrices may be rectangular depending on the number of POD modes used for the different reduced bases. Similarly, the reduced form of the precursor equations can be expressed as where the entries of the reduced operators are computed as As a last step, the reduction of the energy equation is carried out as where the entries of the reduced matrices and sink vector are given by Finally, we note that the elements of input parameter vector µ = (α ext , A V , T ext , P th ) T are simply factors in the reduced order equations as well, thus the reconstruction of the reduced matrices is not necessary for repeating simulations with the ROM. For the solution of the coupled problem, the standard GeN-Foam iteration scheme is used which is discussed in paper [48] in detail. The same iterative scheme is used to solve both the FOM and the ROM models. To compare the solutions of the ROM with those of the FOM we use the following two error indicators: that describes the absolute difference between the largest eigenvalue of FOM and the ROM, and that gives the relative L 2 error in the field variables, where ζ can be Φ g (g = 1, . . . , G), C z (z = 1, . . . , I) or T .

MSFR computational model
A 2D axisymmetric model of the MSFR has been created using the available information in [18,48,49]. The dimensions of the geometry are provided in Figure 1 together with the pre-computed velocity field. The velocity field has been obtained by a standalone steady state solve of the incompressible porous Navier-Stokes equations with k-turbulence model [21]. The heat exchanger (HX) is modeled as a porous medium responsible for flow resistance and heat sink, while the pump (P) is modeled as a simple volumetric momentum source. For more information about the semi-empirical expressions used to compute the parameters of the flow resistance and heat sink, the reader may refer to [21]. The mesh used for the computations has been created using SALOME [50] and contains 11,064 cells. For the neutronics computations, six energy and eight precursor groups are used, meaning that the full-order model has 165,960 degrees of freedom. The corresponding group constants are generated using Serpent 2 Monte Carlo Transport code [51] for two salt temperatures, 900 K and 1200 K. The bounds of the used energy group structure are presented in Table 1.
Altogether 20 parameter vectors are sampled for snapshot generation using the Improved Distributed Latin Hypercube Sampling (IHS) method [53]. The uncertain parameters in these vectors are varied in a ±20% interval around their mean valuesμ = (10 5 W m −2 K, 100 m 2 m −3 , 900 K, 1440 MW th ). Using these snapshots, altogether 16 reduced spaces are created using POD. One for the neutron flux in each energy group, one for the precursor concentration in each group, one for the temperature and one for the logarithmic temperature. The decay in the eigenvalue of the correlation matrices (defined in Eq. (2)) are presented in Figure 2; the eigenvalues are normalized to their largest one. It is visible that for all field variables, the eigenvalues decay rapidly, suggesting that only a few modes are enough to approximate the full-order solution. It can also be observed that, as the half-life of the precursor group decreases (from group 1 to 8), the decay in the eigenvalues of the corresponding correlation matrices is slower.
The dimensions of the reduced spaces for the different solution fields are summarized in Table 2. These numbers are acquired using a strategy based on the energy-retention limit defined as where λ i are the eigenvalues of the correlation matrix built using the snapshots of a selected field, r is the number of POD modes retained for this field and (1 − lim ) is the energy-retention limit. This characterizes the error in the reconstruction of the snapshots originating from discarding the rest of the POD modes (r + 1, . . . , N s ). When lim = 0, all of the POD modes are used, hence the snapshots can be reconstructed exactly. In this work, lim = 10 −7 was used.  The relatively low dimensions can be explained by the fact that all the uncertain parameters appear in the thermal balance equation and changing them does not influence the distributions of flux and temperature fields considerably. This also means that the 165,960 degrees of freedom in the FOM can be reduced to 27 unknowns in the ROM.

Results
After building the reduced-order model, a new realization of the input parameter vector, µ * = (9.3 × 10 4 W m −2 K, 111.1 m 2 m −3 , 833.3 K, 1320 MW th ) is selected and simulations are performed with both the ROM and the FOM in order to test the accuracy and efficacy of the ROM. We stress that µ * was not included in the training set. Figure 3 shows the reconstructed scalar flux in energy group one from the ROM together with its absolute deviation from the full order solution. It is visible that the maximum error is more than three orders of magnitude lower than the maximum value of the full order solution. Figure 4 compares the precursor concentration in precursor group eight. Again, it is visible that the maximum absolute deviation is approximately three orders of magnitude lower than the maximum of the original solution. Figure 5 shows the reconstructed temperature profile of the ROM with its relative deviation from the FOM. Again, the maximum relative error is slightly above 1%. Furthermore, the deviation between the effective multiplication  factors coming from the ROM and FOM was 8.5 pcm, which is also satisfactory. It is worth noting that loosening the energy-retention limit to 10 −6 resulted a 77.6 pcm difference in the eigenvalues, while tightening it to 10 −8 gave 22.3 pcm due to the inclusion of a third POD mode for the logarithmic temperature that decreased the accuracy in the temperature coefficients. By further increasing the energy-retention limit, we observed no further change from the 22 pcm difference level.
Furthermore, the relative L 2 errors in the field variables, defined at equation (30), are summarized in Table 3. It can be observed that none of the fields have a relative L 2 error above 1%. The speed-up factor for the solution of the problem was 1550 for this specific example.

Conclusions
An efficient way of coupling parametrized neutronics and energy equations has been developed for POD-based ROMs. It utilizes the group-wise projection of the multigroup diffusion equations, including the precursor balance equations, and the enthalpy equation in moving fuel systems. In the current state of the development, the velocity equations are not solved; a constant velocity field is assumed. The non-linear temperature dependence of the group constants is handled using DEIM. To demonstrate the viability of the method, it has been implemented in OpenFOAM based multiphysics solver, GeN-Foam. Moreover, a ROM has been created for the MSFR and the corresponding reduced spaces have been generated by exercising a simplified 2D axisymmetric FOM. The limited number of parameters and the underlying physics resulted in reduced spaces with dimensions between 1-3. A test has been performed with four uncertain parameters in the energy equation and the results of the ROM have been compared to those of the FOM. A good agreement is observed both in terms of the solution vectors and the effective multiplication factors. The experienced speed up factor was above 1,500.