Antonio Galia
email: antonio.galia@cea.fr

Richard Sanchez

Igor Zmijarevic
email: igor.zmijarevic@cea.fr

A METHOD OF DYNAMIC HOMOGENIZATION: APPLICATION TO 2D CORE CALCULATION

Keywords: Core Calculations, Homogenization Theory, Equivalence, Discontinuity Factors, Dynamic Homogenization, Domain Decomposition, PWR, APOLLO3 R

Three-dimensional deterministic core calculations are typically based on the classical two-step approach, where the homogenized cross sections of an assembly type are precalculated and then interpolated to the actual state in the core. On the other hand, 3D direct transport calculations and the 2D/1D Fusion method have recently been applied showing excellent agreement with reference Monte-Carlo, but still remaining computationally expensive for multiphysics applications and core depletion calculations. In the present work, we propose a method of Dynamic Homogenization (DH) as an alternative technique that lies between the classical and the direct approaches in terms of accuracy and performance, with no need of off-line calculations and avoiding expensive 3D transport calculations. The method consists in an iterative process between core and assembly calculations, whose main feature is to homogenize the cross sections taking into account the environment of each assembly. To do that, we define a reference homogenization problem that is solved with imposed core eigenvalue and fixed incoming boundary source and that is normalized to the core macro-currents. We have also analyzed the equivalence theory and flux discontinuity ratios techniques in the framework of non-conservative boundary conditions because particular care has to be taken in this case in order to preserve the total net assembly leakage. The current work presents analysis and tests of the Dynamic Homogenization technique for two-dimensional configurations. The NEA "PWR MOX/UO2 Core Transient Benchmark" was chosen to perform the calculations and compare the different methods.

INTRODUCTION

In nuclear reactor physics homogenization techniques play a fundamental role in core modeling because they allow to simplify such complicated systems, for which direct transport (DT) calculations are time demanding and require considerable computational resources. Moreover, the highly detailed information that DT can provide is generally not used for routine industrial applications. The objective of homogenization theory is then: i) to reduce the problem size, either substituting detailed heterogeneous geometries with homogeneous ones, or averaging the energy dependencies; ii) to construct a low-order operator that reproduces average transport quantities (typically the reaction rates) of the original system.

Since we do not know the reference fine transport solution in the core, the most delicate part of point i) is to define a Reference Homogenization Problem (RHP), that typically corresponds to an assembly calculation, whose solution is used as weighting function for cross-sections homogenization. In order to obtain accurate results, the solution of the RHP must then be representative of the actual core conditions.

In the classical two-step approach, the RHP is defined per assembly type with infinite-lattice conditions. The multiplication constant can then be very different from that of a steady-state reactor, which is the only information on the environment that is known at the homogenization stage. The condition k eff = 1 is therefore forced in the assembly transport calculation by introducing a critical leakage model so as to obtain a critical spectrum and the critical reaction rates to be preserved. This model is based on the assumption that the transport solution within an assembly is the product of a periodic solution and a macroscopic distribution: the so called infinite medium fundamental mode assumption. It has been shown that this approximation is correct in the asymptotic situation where the core is composed of a large number of identical assemblies [START_REF] Larsen | Diffusion theory as an asymptotic limit of transport theory for nearly critical systems with small mean free paths[END_REF]. It follows that the assemblies in the core that are the most sensitive to the classical model are those at the peripheral area next to the reflector and those next to other types of assembly such as the UOx-MOx or Rodded-Unrodded interfaces [START_REF] Palmtag | Two-group spectral corrections for MOX calculations[END_REF]. Hence, there is the interest of developing new methodologies that are more reliable for core analysis.

An issue with the homogenization process is the loss of the detailed information on the cross sections. Thus, the coarse operator resulting from the homogenization can at most reproduce reference macroscopic reaction rates in macroregions and macrogroups when provided with a concurrent macroscopic description of references sources and boundary conditions. Moreover, with the exception of very limited cases, flux homogenized cross sections fail to reproduce reference reaction rates and it is then necessary to recur to special techniques, such as equivalence theory (EQV) [START_REF] Kavenoky | The SPH homogenization method[END_REF][START_REF] Hébert | A Consistent Technique for the Pin-by-Pin Homogenization of a Pressurized Water Reactor Assembly[END_REF] or the use of flux discontinuity ratios (FDR) [START_REF] Smith | Assembly homogenization techniques for light water reactor analysis[END_REF][START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF] to reproduce the reference reaction rates in the RHP. Use of these techniques results in a marked improvement: for example, typical pin power errors are reduced from around 10-12% to roughly 5-7% in a LWR at nominal conditions for 2D configurations [START_REF] Mala | Nodal and Pin-by-pin Calculations Comparison with Codes SIMULATE-5 and DYN3D[END_REF].

The question that arises at this point is whether the main source of these errors is due to the inadequacy of the adjusted low-order operator to describe the physics of the full core problem, or that of the reference homogenization problem to provide a good approximation for the core flux, or even both. Concerning the former, in 2011 Grundmann and Mittag [START_REF] Grundmann | Super-homogenisation factors in pinwise calculations by the reactor dynamics code DYN3D[END_REF] computed the NEA PWR UO2/MOX Benchmark [START_REF] Kozlowski | PWR MOX/UO2 Core Transient Benchmark (Final Report)[END_REF] with a 8-group pin-by-pin simplified P 3 operator, with and without equivalence, showing better results with respect to the other benchmark participants who adopted the nodal two-group diffusion, halving the standard deviation of the assembly power relative error, with a maximum of pin power error around 7-8%. Concerning the latter point, many attempts have been made in the past to better account for the assembly environment at the homogenization stage. These can be divided into three classes: i) direct modeling, ii) cross-section correction techniques, and iii) iterative core-assembly calculations.

Direct modeling consists of running colorset calculations, combining each type of assembly in the core in order to homogenize with a more realistic environment. Typically, this approach is adopted for reflector homogenization [START_REF] Mala | EPR Fuel Cycle Depletion with pin-by-pin code Tortin and nodal code Simulate5[END_REF] to improve its response. However, if one has to consider many combinations of assembly types in several configurations of refueling strategy or thermal feedback, the method can be expensive and the homogenized nuclear data can increase considerably.

The second technique approximates the environment effect by correcting the homogenized nuclear data. Rahnema et al [START_REF] Rahnema | Leakage corrected spatial (assembly) homogenization technique[END_REF][START_REF] Rahnema | High-order cross-section homogenization method[END_REF] used a linear perturbation method to the assembly boundary conditions in order to calculate and tabulate the homogenized cross sections for a set of current-to-flux ratios. The latter were computed during the core calculation while the coefficients for the perturbation were determined from off-line assembly calculations, so that in the iterative process of crosssection correction, no lattice calculation needed to be performed. With the same philosophy of offline cross-section generation, Clarno et al. [START_REF] Clarno | Capturing the Effects of Unlike Neighbors in Single-Assembly Calculations[END_REF] used energy, angle, and space dependent albedo boundary conditions to "capture the effects of unlike neighbors in single-assembly calculations". Instead of using the core quantities, as first estimation these albedos were determined solving a 1D transport problem, that was representative of the 2D two-assemblies problem. However, the cross sections used by the 1D solver were homogenized in space using the 2D assembly flux at infinite lattice conditions. A similar approach adopted by Palmtag for UOx-MOx interfaces [START_REF] Palmtag | Two-group spectral corrections for MOX calculations[END_REF] consists of using the two-group homogenized cross sections, the spectral index and several fitting coefficients to evaluate two types of spatially dependent corrections: a leakage correction that accounts for fast spectrum changes due to the leakage effects, and a spectrum correction that accounts for the spectrum changes of the sources. The set of correction coefficients was determined once and for all from a reference case calculation and then used for different problems.

The third class of methods is based on a very different philosophy, because it does not require any off-line calculation. A lattice solver is used as an "on-the-fly" cross-section generator for the core solver. This approach is definitely more expensive but it can be easily parallelized.

The original idea of the iterative core-assembly calculations was, however, to eliminate the fundamental mode assumption, where the boundary conditions for the lattice calculations are calculated by the low-order operator and updated at each iteration. The method was firstly called Dynamic Homogenization (DH) was proposed by Mondot and Sanchez in 2003 [START_REF] Mondot | An iterative homogenization technique that preserves assembly core exchanges[END_REF] and tested for simple cases in 1D problems. In their work, the fine information on the environment came directly from the neighboring assemblies and from the core currents and eigenvalue. Each assembly has then different boundary conditions that account for its position in the core and for the isotopic content of the neighboring assemblies: both sources of gradients in a reactor.

In 2008 Takeda et al. [START_REF] Takeda | Leakage Dependent SPH Factor for PWR Whole Core Transport Calculation[END_REF] applied a similar iterative process but, to redefine the RHP at each iteration, they adopted albedo boundary conditions determined with the core currents instead of using the outgoing angular flux of the neighbor. It ensue that each lattice calculation is independent of the others because it only depends on core quantities, which entails that the spectrum of the neighboring assemblies cannot be properly taken into account. In 2014 Colameco et al. [START_REF] Colameco | Iterative Transport-Diffusion Methodology for LWR Core Analysis[END_REF] extended the DH method to a 2D configuration obtaining good agreement against reference transport calculation for a 2x2 cluster of UOx and MOx assemblies with reflection boundary conditions. The method was also adopted for Pebble Bed Reactors by Grimod and Sanchez in 2015 [START_REF] Grimod | A dynamic homogenization model for pebble bed reactors[END_REF] where the pebbles where depleted with their own fine-group fluxes.

To eliminate the error introduced by the interpolation on the parametrized cross section tables, Varin and Marleau proposed in 2006 [START_REF] Varin | CANDU reactor core simulations using fully coupled DRAGON and DONJON calculations[END_REF] the history based local parameter approach which consisted of dynamically homogenizing the assemblies with the usual infinite lattice model but using the actual physical parameters in the core required in multiphysics and depletion calculations. We consider interpolation errors to be much less important that the error introduced by the use of conservative boundary conditions. The former can be minimized in the classical two-step formalism by increasing the number of tabulated points; moreover, with the use of modern data compression techniques [START_REF] Tomatis | A multivariate representation of compressed pin-by-pin cross sections[END_REF] which faster reconstruction, which achieve between 80-90% reduction of the library size, such an approach would not have an important impact in the library size. In spite of using a dynamical homogenization, the history based technique is not in the spirit of our present approach, which consists of using more realistic boundary conditions for the assembly in order to correct the error introduced by the infinite lattice model. Clearly, our approach homogenizes the assembly using the actual core parameters.

The method of domain decomposition not only offers the numerical advantage of an easy implementation of moderate to massive parallelization, but also provides a natural paradigm for homogenization techniques for full core calculations using a low-order operator such as diffusion. By progressively increasing the degree of sophistication in the representation of the boundary condition: one starts at the lowest level with the classical infinite lattice model used in the two-step method, and move through intermediary homogenizations, such as the albedo method, to finally reach the method of dynamic homogenization presented in this work, where each assembly is provided with incoming fluxes that are obtained from the fluxes exiting neighboring assemblies or from the adjacent reflector. Moreover, the latter is also decomposed into "assembly" like domains which are homogenized in a similar fashion, providing thus for a physically correct reflector homogenization based on the transport detailed spectrum for the angular flux entering the reflector from peripheral assemblies Further up in the refinement of the boundary condition for the local assembly homogenization problem one finds the 2D-1D fusion method and the exact 3D core transport calculation, where the homogenization is used to construct a low-order nonlinear acceleration method for the core transport iterations, such construction is necessarily dynamic [START_REF] Joo | Dynamic implementation of the equivalence theory in the heterogeneous whole core transport calculaio[END_REF] because it has to be renewed each time the transport flux changes. However, although both homogenization techniques search to construct a global core low-order operator, there is an essential difference between them. While the acceleration technique forces the coarse operator to satisfy coarse transport balance, usually by modifying diffusion with a drift term, the dynamic homogenization, which is based in the original formulation in [START_REF] Mondot | An iterative homogenization technique that preserves assembly core exchanges[END_REF], modifies the transport assembly solution by normalizing the entering angular fluxes from neighboring assemblies so as to preserve the partial macro group currents predicted by the low-order core calculation, thereby preserving the low-order core eigenvalue and exchanges between assemblies. In few words, the acceleration methods converge to an approximated transport solution, in the case of the 2D-1D fusion, or to the exact 3D transport, while the homogenization techniques converge to a low-order core calculation. We note that, with the exception of the direct transport calculation, all other methods compute the low order core operator using assembly homogenization calculations for 2D assemblies. For 3D core calculations this might be a rather strong homogenization approximation.

The article is organized as follows: in Section 2 we review the most popular methods used for assembly homogenization, with emphasis on their application to the DH methodology; Section 3 discusses the main features of the classical two-step approach, showing the approximations introduced by the calculation scheme and the difficulty of reflector homogenization; a detailed description of the method of Dynamic Homogenization, its calculation scheme and the approximation used for the assembly boundary conditions are given in Section 4; numerical tests are performed on a 2D PWR core with different types of assemblies and burnups inspired by the NEA PWR MOX/UO2 Core Transient Benchmark [START_REF] Kozlowski | PWR MOX/UO2 Core Transient Benchmark (Final Report)[END_REF]. The results of these calculations are presented in Section 5 and used there for a comparison of the precision and the computational cost of DH against direct transport (considered here as reference) and the classical two-step approach; Finally, conclusions and further discussion of the advantages and disadvantages of DH and on the choice of a two-group diffusion operator as core solver are given in Section 6.

ASSEMBLY HOMOGENIZATION

A homogenization process requires, above all, the definition of a reference homogenization problem that provides the macroscopic quantities to be used in the construction of the low-order operator. As mentioned in the introduction, the objective of a homogenization technique is to construct a coarse operator that reproduces the average reaction rates and, therefore, the particle balance for the RHP, under the assumption that these quantities are close to those that would extant in the core. Regarding this last condition, a natural framework for the construction of the RHP is supplied by the exact solution for the flux in the RHP considering the RHP in its actual situation in the reactor core [START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF]. This flux is the solution of a source transport problem with fixed eigenvalue, where the eigenvalue is the core eigenvalue and the source is the angular flux entering the assembly. For the calculation of reactor cores in normal operation conditions a good approximation for the former is k ef f = 1, however the exact incoming angular flux can only be provided by a full exact solution of the entire core. Hence, a determining factor in the construction of the RHP for assembly homogenization will be to construct a realistic entering angular flux.

The RHP is typically solved with a fine-transport operator that can treat all the heterogeneities of the geometry, such as the Method of Characteristics (MOC) or Monte Carlo [START_REF] Ilas | A Monte Carlo based nodal diffusion model for criticality analysis of spent fuel storage lattices[END_REF]. The most common types of homogenization are full assembly or pin-by-pin. The flux-weighted homogenized cross sections associated to each coarse region are defined as:

Σ h,G x,i,R = g∈G r∈R Σ g x,i,r Φ g r V r Φ G R V R = τ G x,i,R Φ G R V R , (1)
Φ G R = g∈G r∈R Φ g r V r V R , (2)
where Σ h,G x,i,R is the homogenized macroscopic cross section for reaction x and isotope i in the coarse region R and coarse energy group G, while Σ g x,i,r is the macroscopic cross section in region r and group g of the fine discretization for the RHP calculation. The τ G R and the Φ G R defined in Eq. ( 2) are respectively the reference reaction rate and average scalar flux for R and G. A rigorous treatment would require the angular flux as weighting function but it would result in angulardependent homogenized cross sections, so generally it is not applied. Anisotropic transfer cross sections are typically homogenized with the scalar flux instead of using the angular moments.

A remark is necessary at this point. We note that the homogenized cross sections are obtained by dividing the reference reaction rate by the volume of the macro-region times the averaged flux in the macro-region. For the flux-weighted cross sections defined in Eq. (1) this flux is the reference flux but, as we discuss soon, it is often necessary to use a different flux, such as the flux computed with the low order operator. What is important to note is that, regardless of the flux used to obtain the homogenized cross sections for all reactions and isotopes, all these cross sections can be obtained from the set of total cross sections Σ h,G R . Indeed, for any reaction and isotope we have

Σ h,G x,i,R = τ G x,i,R τ G R Σ h,G R , (3)
where τ G R is the total reference reaction rate and the ratio of reference reaction rates is provided by the reference transport solution.

It follows from relations (3) that the only unknown coarse cross sections to be determined are the total cross sections Σ C,G R for all macro region and macro groups. These cross sections are conditioned by the associated reaction rate conservation equations

Σ C,G R V R Φ C,G R ( - → Σ C ) = τ G R , ∀R, G, (4)
where the fluxes Φ C,G R are obtained from the solution of the coarse RHP and are therefore functions of all the total coarse cross sections

- → Σ C = {Σ C,G R , ∀R, G}.
Moreover, coarse boundary conditions will be also defined so as to reproduce the macroscopic boundary values of the fine reference problem. Finally. if the RHP is an eigenvalue problem we will also request the preservation of the eigenvalue. However, the latter can be superfluous, albeit numerically advantageous, because, when the coarse operator is constructed so as to preserve net global leakage, then global neutron balance ensures that the eigenvalue is automatically preserved.

The starting point for our discussion is the global neutron balance equations in the entire geometrical domain D of the RHP. For a given macro group G these equations for the fine and the coarse RHP problems read:

S∈∂D A S J +,G S - S∈∂D A S J -,G S + R∈D Σ h,G R V R Φ G R = R∈D V R Q G R , (5)
S∈∂D A S J +,C,G S - S∈∂D A S J -,C,G S + R∈D Σ C,G R V R Φ C,G R = R∈D V R Q C,G R , (6)
Q = H 0 Φ + 1 λ F Φ, (7)
Q C = H C 0 Φ C + 1 λ C F C Φ C . (8)
In these equations, the boundary ∂D of D has been partitioned into a set of macro surfaces noted S, C stands for coarse, Q and Q C are the isotropic sources comprising scattering (H 0 ) and fission (F ) as in Eqs. ( 7) and ( 8) and J ± S are the outgoing (+) and incoming (-) average partial currents crossing surface S of area A S :

J ± S = 1 A S s∈S A s J ± s , (9)
where s denotes a boundary surface of the fine transport mesh. We also note that the balance equations have been obtained by direct integration over the entire geometric domain of the original equations (transport, diffusion, etc.) and express a relation between the cross sections and the solution (fluxes, currents) of the equations constructed with those cross sections. Therefore, these equations are satisfied for any regular set of cross sections for which the equations have a solution.

This means that the conservation relations are valid regardless of the constraints (3) and/or (4) or any other constraint that one might impose between the cross sections, as long as the equations admit a solution. In the following, we shall use total leakage and local leakage to refer to the leakage of the whole domain D or to the leakage associated to a coarse region R, respectively.

We discuss now the construction of the coarse operator whose cross sections satisfy constraints (3) and ( 4), for which the only independent cross sections are the total ones. Consider first the case when the coarse operator has the same total net leakage as the fine operator. Then, macro-region reaction rates conservation ensures that the eigenvalue is also preserved. Furthermore, it is also possible to show that the constraints in (4) are linearly dependent. We follow here an argument given in Ref. [START_REF] Sanchez | Diffusion Piecewise Homogenization via Flux Discontinuity Ratios[END_REF], whereupon the nonlinear problem defined by Eq. ( 4) is replaced by an equivalent homogenization problem where the detailed sources Q C,G (r) are those obtained from the solution of the original problem. Hence, the modified problem split into a set of one-group source problems each of which satisfies balance Eq. [START_REF] Mala | Nodal and Pin-by-pin Calculations Comparison with Codes SIMULATE-5 and DYN3D[END_REF]. It follows that the sum of the macro-region reactions rates is conserved for each macro-group,

R∈D Σ C,G R V R Φ C,G R = R∈D τ G R , ∀G, (10)
as a comparison of Eqs. ( 5) and ( 6) shows, and this conservation is satisfied regardless of the values adopted for the total homogenized cross sections. * The implication is that the homogenization problem has one degeneracy per macro group.

Conservation of total net leakage is always true for the traditional infinite lattice model, for which the leakage contributions vanish for both the reference and the coarse calculations. However, if the RHP has non conservative boundary conditions the preservation of the total net leakage cannot be always incorporated in the definition of the coarse RHP problem and, ultimately, it all depends on the type of boundary condition that the coarse operator can support. Transport-like operators require incoming fluxes and therefore the coarse boundary condition can preserve only the incoming reference currents in Eq. [START_REF] Smith | Assembly homogenization techniques for light water reactor analysis[END_REF]. For these operators it is the sum of the total leakage and of the total reaction rate per macro group that is unconditionally preserved and, therefore, constraints Eq. ( 4) are not degenerated. Note that the implication is that the solution of the homogenization problem will preserve the total exiting current but not the detailed exiting current over each surface. By imposing that the eigenvalue is preserved as well as the reaction rates, we ensure the existence of a single coarse operator. On the other hand for elliptic operators such as diffusion it is possible to use the total reference currents as boundary conditions and achieve then net leakage conservation.

Here, the eigenvalue is automatically preserved but we have again a degeneracy per macro group and the homogenization problem accepts an infinite † number of solutions.

If ones uses the customary flux-weighted homogenized cross sections in Eq. ( 1) for the construction of the coarse operator, then it is very unlikely that Eqs. (4) would be satisfied. The exception is full assembly homogenization or piecewise homogenization of an assembly consisting of a set of identical pins, in which case global balance (or local balance for the assembly with identical pins) ensures the conservation of reaction rates. However, for pin-by-pin or, more generally, piecewise homogenization, the individual macro-region reaction rates are most likely not preserved. One may conclude that, in general, the coarse operator constructed by flux-weighting homogenization fails to reproduce macroscopic reaction rates.

To overcome this problem one must construct a faithful coarse operator, meaning an operator that preserves reference reaction rates for each macro-region and each macro-group as in Eq. (4). Two techniques have been used to attain this goal: equivalence theory (EQV) and homogenization via Flux Discontinuity Factors (FDF). The former consists of directly computing an optimal set of cross sections, via iterative solution of nonlinear Eqs. 4. The latter, which applies only to elliptic coarse operators, consists of using flux-weighted cross sections, while introducing new degrees of freedom in the homogenization process, leading to a drastic modification of the coarse operator.

Equivalence theory

The idea behind equivalence theory, originally proposed by Kavenoky [START_REF] Kavenoky | The SPH homogenization method[END_REF] and later refined by Hébert [START_REF] Hébert | A Consistent Technique for the Pin-by-Pin Homogenization of a Pressurized Water Reactor Assembly[END_REF] and Sanchez et al. [START_REF] Sanchez | APOLLO2: a user oriented, portable, modular code for multigroup transport assembly calculations[END_REF], is to obtain a faithful homogenization by numerically searching for a set of homogenized cross sections, called equivalent cross sections Σ eq , which lead to a faithful coarse operator. Because of the dependence of the coarse fluxes on the coarse cross sections, the problem is nonlinear and it has been solved either i) by iteratively minimizing the functional

F( - → Σ eq ) = R,G 1 - τ C,G R ( - → Σ eq ) τ G R 2 , (11)
where -→ Σ eq is the set of equivalent cross sections for all macro regions and macro groups. The solution can be obtained for instance with a Newton algorithm [START_REF] Ortensi | A Newton solution for the Superhomogenization method: The PJFNK-SPH[END_REF], or ii) by fixed-point iterations with Eq. (4), which is then used to update the equivalent cross sections. In both methods, the iterations are initialized with the flux-weighted cross sections and the updated values of the equivalent cross sections are used to compute the new flux Φ C,G R ( -→ Σ eq ), solution of the associated coarse equation. The difference between the two methods consists in the way the cross sections are updated: while in the functional method the set of equivalent cross sections is updated by a global minimization of the functional, in the fixed-point method each equivalent cross section is obtained directly from its associated reaction rate according to Eq. (4).

As earlier explained, when conservative boundary conditions are applied, and in diffusion also for non conservative conditions, the homogenization problem is degenerated [START_REF] Ortensi | A Newton solution for the Superhomogenization method: The PJFNK-SPH[END_REF][START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF]. In order to choose one single solution, two different types of normalization are used in practice. The first consists of imposing the conservation of the domain's average reference scalar flux per macro-group, while in the second, known as Selengut's normalization [START_REF] Hébert | Development of a Third-Generation Superhomogénéisation Method for the Homogenization of a Pressurized Water Reactor Assembly[END_REF], it is the total transport incoming current per macro group that is preserved. In the particular case of full assembly homogenization with normalization of the average scalar flux, one has Σ eq = Σ h . The advantage of equivalence theory is that it can be used with whichever low-order operator is adopted for the coarse calculation, which can also be a transport operator. However, it has been shown that the fixed point iterations may not converge, in particular for cases where transport effects are dominant [START_REF] Ortensi | A Newton solution for the Superhomogenization method: The PJFNK-SPH[END_REF].

Flux discontinuity factors

The homogenization via flux discontinuity factors was analyzed by Koebke [START_REF] Koebke | A New Approach to Homogenization and Group Condensation[END_REF][START_REF] Koebke | On the Reconstruction of Local Homogeneous Neutron Flux and Current Distributions of Light Water Reactors from Nodal Schemes[END_REF], Smith [START_REF] Smith | Assembly homogenization techniques for light water reactor analysis[END_REF] and Sanchez [START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF]. This technique was initially applied for full assembly homogenization but we discuss it here in the general context of piecewise homogenization. The basic approach is to force the preservation of the reference net currents at each interface S of a macro-region so that, if fluxweighted homogenized cross sections are used in the coarse operator, the eigenvalue, the reactions rates and the local leakage are simultaneously reproduced for each macro-region and macro-group. However, as we mentioned earlier, in these conditions is very improbable that the resulting solution will be continuous at macro-region interfaces. This difficulty was circumvented thanks to Koebke's idea of introducing a relaxation in the continuity condition of the flux at each mesh interface [START_REF] Koebke | A New Approach to Homogenization and Group Condensation[END_REF] to allow for a discontinuity of the flux at the interface S between two adjacent macro-regions A and B:

f S A Φ C S A = f S B Φ C S B ⇒ Φ C S A = r AB Φ C S B , (12)
where r AB is called Flux Discontinuity Ratio (FDR). Note that r AB = 1 when no discontinuity is introduced, as well as at the boundaries with reflection conditions. Also, in Eq. ( 12)

f S A = Φ S Φ C S A (13)
is the Flux Discontinuity Factor ‡ (FDF) defined as the ratio of the reference to the coarse averaged interface fluxes on surface S.

Koebke [START_REF] Koebke | A New Approach to Homogenization and Group Condensation[END_REF] was the first to use local balance for a macro region to show that if one wanted to simultaneously use flux-weighted cross sections and preserve reaction-rates, then the reference local net currents should also be preserved. He also concluded that this will also preserve the averaged macro region fluxes and, furthermore, recognizes that most probably the resulting scalar flux will not be continuous at the interfaces between macro regions. The solution thus was to implement the preservation of the reference transport net current at each interface between two macro regions and accept flux discontinuity at the macro region surfaces. In his seminal work, Koebke decided to impose the somewhat artificial condition that the discontinuity factors should be equal on the two opposite sides of each homogeneous rectangular macro region. In order to achieve this he had to introduce macro region and direction dependent diffusion coefficients. For a two-dimensional homogenization mesh, the homogenization was thus achieved by adding four equivalence parameters per macro region: a direction-dependent diffusion coefficient and a bilateral discontinuity factor for each one of the directions. Also, the diffusion coefficients have to be iteratively determined.

In order to avoid this iterative process, Smith [START_REF] Smith | Assembly homogenization techniques for light water reactor analysis[END_REF] decided to use the typical volume averaged diffusion coefficient and defined his Generalized Equivalence Theory (GET) by relaxing the artificial symmetry constrain required by Koebke's homogenization and allowing for the two discontinuity factors associated to a given direction to be independent of each other. Hence, in Smith's GET a discontinuity factor is associated to each interface that a macro region may have with a neighbor macro region, which is simply calculated via its definition in Eq. [START_REF] Clarno | Capturing the Effects of Unlike Neighbors in Single-Assembly Calculations[END_REF]. For a rectangular homogenization mesh this amounts to four FDFs, two for each direction. ‡ Named heterogeneity factor by Koebke.

The question that arises next is how to determine these factors when the reference core solution is not known. Smith made a distinction between Reference Discontinuity Factors (RDF) and Assembly Discontinuity Factors (ADF), calculated using the core solution and the assembly solution of the RHP, respectively. Following Koebke's steps, Smith showed for a 1D problem that if the RDFs are used with the homogenized cross sections weighted with the reference core flux, then the coarse operator reproduces exactly the reference average quantities. This suggests the possibility to construct an operator, as CMFD, which can be used as a nonlinear acceleration for the transport solution.

However, the objective of a homogenization paradigm is to predict the reference heterogeneous solution without actually solving the entire core problem. As Koebke has previously done by comparing the performance of his FDFs calculated with different environments of a given assembly type, Smith compared the impact of the RDFs calculated for each assembly type in different core environments, and showed that the results obtained using the ADFs computed with an infinite lattice model gave a good estimate of the average of all the actual RDF corresponding to that type of assembly in the core [START_REF] Smith | Assembly homogenization techniques for light water reactor analysis[END_REF]. These ADFs cannot be exact because they do not account for the position or for the environment of the assembly. To determine the ADFs, he observed that for full assembly homogenization with periodic boundary conditions the coarse flux is constant. It follows that the coarse surface averaged flux in Eq. ( 13) is equal to the coarse volume averaged flux, which in turn is also equal to the reference volume averaged flux. Hence, Φ C S = Φ C = Φ so that the resulting ADFs are independent of the coarse solution for the RHP:

f GET = Φ S Φ RHP , (14)
where | RHP indicates evaluation for the RHP.

More recently, Sanchez [START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF] proposed a Black-Box (BB) homogenization model based on the use of flux-weighted cross sections wherein, instead of the reference averaged flux, it is the reference partial currents that are preserved at each interface. One of the advantages of this approach was to deal with diffusion as well as transport coarse operators. While in transport there were two Current Discontinuity Factors (CDFs) for each of the two opposite interface partial currents, for the diffusion operator one could use a single equivalent FDF, so that

J ±,C S = f Φ C S 4 ± J S 2 . ( 15
)
to achieve the conservation of the partial currents at the interface: J ±,C S = J ± S . From these conservations conditions one derives the value of the equivalent FDF:

f S A = 2 J + S + J - S Φ C S A . ( 16
)
to be compared with Koebke and Smith's definition in Eq. [START_REF] Clarno | Capturing the Effects of Unlike Neighbors in Single-Assembly Calculations[END_REF]. Moreover, in the case of infinite lattice boundary conditions one has J S = 0 so that J + S = J - S and for full assembly homogenization one obtains the simple formula

f BB = 4 J ± S Φ RHP , (17)
which is also independent of the discretization used with the coarse operator. The GET and BB homogenization models in Eqs. ( 14) and ( 17) give equal FDFs only in the extreme cases when the initial assembly is homogeneous or when the reference operator is also diffusion.

The result of the equivalent FDF is to introduce a jump from the partial currents at the interface of the homogenized macro region:

∆J ±,C S = ±(J ±,C S -lim →0 + J ±,C S ) = ± (f -1)Φ C S 4 RHP , (18)
where S denotes the surface parallel at the interface at a distance towards the interior of the macro region. As was previously observed [START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF], this jump is equivalent to inserting in the coarse model a singular surface anisotropic source emitting ∆J ±,C S neutrons in the outward direction and as many in the inward direction. We note that the GET FDF leads also to the interface values in Eq. ( 15) and therefore exhibits also a partial current jump as in Eq. ( 18), the only difference in the jump values is in the value of the flux discontinuity factor f . The general case of piecewise homogenization, where the coarse operator has both internal and boundary interfaces, was recently investigated by Sanchez et al. [START_REF] Sanchez | Diffusion Piecewise Homogenization via Flux Discontinuity Ratios[END_REF] who computed the flux discontinuity ratios in Eq. ( 12) induced by the simultaneous preservation of the net averaged currents at all the internal interfaces of the RHP, while using the reference surface averaged currents at the external interfaces as boundary conditions for the coarse RHP. These FDRs can be calculated from the minimization of the functional

F( - → r ) = S,G J C S ( - → r ) J S -1 2 , (19)
where -→ r is the set of all FDRs, one per internal interface and macro group. Moreover, they considered the case when the diffusion discretization mesh is equal to the homogenization mesh and the numerical discretization of the diffusion operator has only one degree of freedom per interface, and computed an analytical solution for finite differences and nodal transverse discretizations. For these numerical schemes the degree of freedom can be identified with the interior limit of the averaged interface flux and these surface averaged fluxes can then be calculated from the solutions of an independent set of local diffusion problems, one per macro region, where each macro region is solved with the boundary conditions provided by the reference interface averaged currents. Finally, the FDRs are calculated from the surface averaged fluxes of the two macro regions that share each interface.

When the diffusion discretization mesh is finer than the homogenization one, the solution cannot be computed analytically and has to be determined from the iterative minimization of the nonlinear functional in Eq. [START_REF] Tomatis | A multivariate representation of compressed pin-by-pin cross sections[END_REF]. Numerical calculations showed that, while with finite differences the FDRs slowly converged with the number of sub meshes, this was not the case for the transverse nodal method, given thus a numerical advantage to the latter for which analytically determined FDRs for one sub mesh are a very good approximation for any other number of sub meshes.

In the context of piecewise homogenization, the diffusion solution depends only on the flux discontinuity ratios and not on the individual flux discontinuity factors at both sides of each internal interface. Therefore, whether one uses GET FDFs as in Eq. ( 13) of BB FDFs as in Eq. ( 16) does not make any difference in the final values of the FDRs. This follows from the fact that the numerators of these two FDFs evaluations, 2(J + S + J - S ) and Φ S , are reference transport values which depend only on the interface and not on the adjacent macro regions and, therefore, simplify out when evaluating the FDR. It follows that the FDRs depend only on the method and the mesh adopted for the discretization of the coarse problem. However, for assemblies homogenized in separate RHPs one has to provide also the flux discontinuity factors at the boundary surfaces which will be used in the final whole core diffusion calculation, and it is here that the use of GET or BB makes a difference. The conclusion to this discussion is that while the FDRs are mathematically defined, the FDFs are based in a physical assumption regarding the behavior of the coarse flux at the interface between two adjacent assemblies: the GET model enforces the continuity of the reference flux, while the BB model preserves the exchanges between the assemblies as described by the reference partial currents.

When compared with equivalence theory, the homogenization technique based on FDFs has not only the advantage of improving the solution of core calculations, but it also offers the simplicity of using the physical transport flux for cross section weighting, plus the potential for a direct evaluation of the homogenization parameters that avoids problematic iterative solutions of nonlinear problems that might not always converge.

CMFD acceleration

With the same incentive of preserving the reference net currents at each mesh interface, which allows to use flux weighted cross sections while preserving reaction rates, another widely used technique is the Coarse Mesh Finite Difference method (CMFD) [START_REF] Smith | Determination of homogenized diffusion theory parameter for coarse-mesh nodal analysis[END_REF][START_REF] Smith | Nodal method storage reduction by nonlinear iteration[END_REF][START_REF] Aragonés | A Linear Discontinuous Finite Difference Formulation for Synthetic Coarse-Mesh Few-Group Diffusion Calculations[END_REF][START_REF] Smith | Full Core, 2-D, LWR Core Calculations with CASMO-4E[END_REF][START_REF] Kelley | CMFD acceleration of spatial domain-decomposed neutron transport problems[END_REF][START_REF] Jung | Practical numerical reactor employing direct whole core neutron transport and subchannel thermal/hydraulic solvers[END_REF], which today it is often used as nonlinear acceleration for transport due to its fast solution and efficiency. As opposed to the FDF technique, the CMFD flux is continuous at the interfaces while an artificial drift term is added to the diffusion operator in order to preserve local balance. This is mediated by correcting the definition of the net current:

J AB = -D AB (Φ B -Φ A ) -DAB (Φ B + Φ A ). (20)
In this equation the first term is the expression for the traditional Fick's law obtained with a meshcentered finite differences scheme,

D AB = 2D A D B D A ∆x B + D B ∆x A , (21)
where D Z and ∆x Z for Z ∈ A, B are the diffusion coefficient and the size of a mesh element, respectively, and A and B stand for two adjacent meshes with B to the right of A. The diffusion coefficient D in a region is typically calculated by homogenizing in space and energy 1 3Σtr by fluxvolume weighting according to Eq. ( 1). The coefficient DAB of the corrective term is evaluated instead by substituting the reference transport currents and scalar fluxes of the previous iteration in Eq. [START_REF] Joo | Dynamic implementation of the equivalence theory in the heterogeneous whole core transport calculaio[END_REF]. Like in the FDF technique, for a given diffusion coefficient in a mesh, one has to define, respectively for 2D or 3D geometries, four or six equivalence parameters (one per interface). This does not apply, however, to the partial current-based Coarse Mesh Finite Difference (p-CMFD) method [START_REF] Yuk | Whole-Core Transport Solutions with 2-D/1-D Fusion Kernel via p-CMFD Acceleration and p-CMFD Embedding of Nonoverlapping Local/Global Iterations[END_REF], where two corrective coefficients per interface are introduced in order to independently preserve the two partial currents, and therefore the net current. When the CMFD is used as an acceleration, the transport fluxes on the fine spatial and energy meshes are calculated at the end of each power iteration via rebalancing:

ψ g,l+1 r,h = ψ g,l+1/2 r,h Φ C,G,l+1 R g∈G r∈R Φ g,l+1/2 r , ∀r ∈ R, ∀R ∈ D, ∀h, ∀g ∈ G, (22)
where r,h,g are respectively the fine region, the angular moment and the energy group indexes of the angular flux ψ, l is the iteration index and Φ is the scalar flux. The sources for the next iterations are calculated with the accelerated flux and at convergence the flux ratio of Eq. ( 22) must be equal to one. When a method of domain decomposition is used to solve the transport problem, the same rebalancing equation is applied to the outgoing angular flux of a subdomain before transmitting the information to the its neighbors:

ψ +,g,l+1 s,d = ψ +,g,l+1/2 s,d Φ C,G,l+1 R g∈G r∈R Φ g,l+1/2 r , ∀s ∈ ∂R, ∀∂R ∈ ∂D, ∀d, ∀g ∈ G, (23)
where + stands for outgoing, d is the direction index and s is the fine-transport surface element at the boundaries of the coarse region ∂R.

In our analysis we will compare the classical and the DH approaches using the EQV and the FDF techniques, which are the most popular in industrial applications, while for the reference direct calculation CMFD is used for the acceleration of the outer iterations.

CLASSICAL TWO-STEP APPROACH

The classical two-step calculation scheme has been widely used for core modeling and design thanks to its efficacy of producing fast results at different core configurations. The calculation scheme is shown in Fig. 1. The first of the two steps consists of "off-line" 2D RHP calculations for each assembly, performed, as earlier mentioned, on fine spatial and energy meshes with high order angular dependency. A library is produced where the homogenized microscopic cross sections are stored for different physical states and for the most significant isotopes that take part in the fuel depletion process, plus a residual that takes into account the remaining isotopes. Each state is associated to a set of physical parameters, such as burnup, moderator density, fuel temperature, etc. In the second step, a 3D diffusion core solver is used for the actual 3D core depletion calculation. At this stage, the code searches and interpolates the library data to the actual state of the reactor, and performs the 3D calculation on coarse spatial and energy meshes.

Typically, the equivalent cross sections or the flux discontinuity ratios are determined and stored at the end of each lattice calculation and, in the second step, the parameterized library containing all the data necessary for the construction of the coarse operator is fed into the coarse core calculation, so that there is no need for the "equivalence box" in Fig. 1. However, this requires that at the homogenization stage one has to know in advance which low-order operator will be used in the second step and, for the general case of piecewise homogenization with EQV or FDR with submeshing, the necessary coarse solver routines must be included in the same package than the lattice code. Another possibility to organize the two-step procedure is the one that is fully illustrated in the left of Fig. 1 and in which the "equivalence box" is active. Here, the first step is used to compute a first parameterized library containing all transport results necessary for the homogenization stage, including volume and surface averaged fluxes, surface averaged partial currents and flux-weighted cross sections. In the second step, this library is used to perform the homogenization process and fabricate the final parameterized library that is then fed into the coarse core calculation. This operation is represented by the "equivalence box" which has direct access to the appropriate coarse solver routines available in the coarse core code. This procedure is more flexible in that it leaves the user free to adopt different core solvers, with the only limitation of fixing the number of coarse groups in the first step.

Let us enumerate the approximations introduced by the two-step approach. At the RHP calculation stage, the surrounding environment of the assembly is not known, so conservative boundary conditions are typically imposed, as if the motif was repeated to infinity (the so called infinite lattice calculation). This is the first approximation of the two-step approach, since 1) it assumes that the net exchanges between assemblies are negligible. It ensues that the solution is symmetric in space and that the eigenvalue of the problem can be very different from that of a steady-state core. Therefore, the weighting spectrum is typically adjusted by introducing a homogeneous critical-leakage model, so as to obtain an assembly eigenvalue equal to one [START_REF] Sanchez | APOLLO2 Year 2010[END_REF]. This model was developed on the assumption that 2) the transport solution is the product of a periodic fine distribution and a macroscopic distribution, the so called infinite medium fundamental mode assumption.

The sequence of self-shielding and critical lattice calculations is then repeated for each physical state to produce homogenized data. It is a usual practice 3) to perform the depletion calculation for only one assembly history with physical parameters that are typically the nominal ones, and successively using the isotopic concentrations of each burnup step for all the other branch calculations. That is why in Fig. 1 the burnup parameter has a separate loop, associated to the depletion step, while an internal loop is used to do all the computations varying the values of the other parameters with fixed burnup value. This approach may fail when an assembly is actually depleted at physical conditions that are far from the nominal ones, like in the case of a rodded assembly [START_REF] Tomatis | Quantification of history effects in PWR modelling[END_REF]. This is the case also for BWRs where the vapor quality significantly changes axially, which requires to compute depletion histories for different values of that parameter.

An important issue that arises in the classical calculation scheme is how to define an adequate RHP for the reflector homogenization. Indeed, the reflector requires a different homogenization approach because no multiplication occurs and vacuum conditions are imposed on at least one boundary. Regarding this problem, there is a large number of procedures that have been adopted in the literature, such as 1D models, 2D calculations with one fuel assembly next to a reflector "assembly", or one row of 8 or 9 fuel assemblies next to each other including the reflector at the end, or even bigger motifs of 5 × 5 assemblies containing a reflector layer surrounding a mini-core. However, 4) the use of these reflector cross sections introduces further approximations because the RHPs used to compute them may contain assemblies with different fuels or irradiated at different conditions than those in the actual configuration in the core.

In addition to σ h,G x,i,R , one may want to store in the library the isotopic concentrations at each burnup step as well as the transport scalar flux integrated in energy and space, in order to be able to recalculate the reference reaction rates when performing the core calculation. Moreover, in order to calculate the flux discontinuity factors, one also needs to store surface quantities. In case of full assembly homogenization, the pin power form factors may also be stored, in order to reconstruct the power within a coarse mesh by 5) applying power reconstruction techniques that preserve the coarse nodal pin power results. Clearly, the size of the library containing all this information from the fine-transport calculation can quickly become huge. That is why in practice only a small set of values for each physical parameter are only considered, and 6) interpolation techniques are applied because the actual state of the assembly in core can hardly be represented by the values of the parameters used to construct the library.

DYNAMIC HOMOGENIZATION APPROACH

The method of Dynamic Homogenization is an iterative process between assembly and core calculations which does not require any "off-line" computations. The purpose of the method is to eliminate all the approximations and limitations of the classical two-step scheme discussed in Section 3, while avoiding expensive 3D transport calculations; it offers thus an interesting and promising alternative technique for the two-step core calculations.

The calculation scheme of DH is shown in Fig. 2. In the present work only static calculations for a given set of macroscopic cross-sections are considered and self-shielding and depletion calculations in the external loops of Fig. 2 are not studied. The main feature of the method is the use of a reference homogenization problem that accounts for the environment of each assembly in the core, including also reflector and external structure "assemblies". As we shall see, in 2D the method can be seen as a mirror image of a nonlinear acceleration scheme for a direct transport calculation based on the Domain Decomposition Method (DDM).

We start by reviewing the use of the DDM for the solution of the whole core transport equation. This discussion is also based on the diagram in Fig. 2. Let D be the global geometrical domain, typically the reactor core including surrounding materials and reflector, with boundary Γ and N D the number of subdomains, which comprise fuel and reflector "assemblies". Thanks to the DDM, the global core transport problem becomes a set of coupled local assembly transport problems which is iteratively solved in parallel. At each iteration the solution of the local assembly problems with incoming boundary conditions is computed. In order to reproduce the reference solution obtained without DDM, two new conditions have to be imposed: 1) force the same eigenvalue in all the subdomains and 2) enforce the continuity of the angular flux at the interfaces between subdomains. The first condition is not necessary but ensures a robust convergence. The second condition is achieved by iteratively replacing the fluxes entering each assembly (ψ in,i ) with those exiting the neighboring assemblies (ψ out,j ):

     L i ψ i = H i ψ i + 1 λ P i ψ i , r ∈ D i , ψ in,i = ψ out,j , r ∈ Γ ij = Γ i ∩ Γ j , ψ i = β i ψ i + ψ in,i , r ∈ Γ i ∩ Γ, (24)
where D = ∪ i=1,N D D i is the partition of the core into N D subdomains, D i and Γ i are, respectively, the subdomain i and its boundary and Γ ij = Γ i ∩ Γ j is the common interface between neighboring assemblies i and j. Also, in this equation L = Ω • ∇ + Σ, operators H and P stand for scattering and fission production, β accounts for an albedo boundary condition and r is a generic position vector. The outermost loop of the iterative process, where interface conditions and eigenvalue are updated, consists of global iterations, while the fission source is locally updated at the subdomain level in the outer-iteration loop; in this loop a multigroup assembly transport problem with fixed eigenvalue and incoming angular fluxes is solved in parallel for each subdomain. In practice, the maximum number of outer iterations is fixed to one because numerical experimentation has shown that in the DDM it is not worth to converge the multigroup problem, in particular if the boundary source is far from its converged value.

A multigroup full core coarse operator can be used between two global iterations for the acceleration of the eigenvalue and of the fission and boundary sources. For nonlinear acceleration it may require a modification of the coarse operator to ensure that the acceleration process converges, as it is the case for the popular CMFD acceleration where an artificial drift term is added to Fick's law. The multigroup acceleration is mediated by using the solution of the global coarse eigenvalue problem in the rebalancing Eqs. ( 22) and [START_REF] Sanchez | APOLLO2: a user oriented, portable, modular code for multigroup transport assembly calculations[END_REF]. To summarize, each subdomain calculation is solved with imposed multiplication constant λ C , and fixed incoming boundary sources ψ in,i , where the former, being an integral parameter of the whole domain, is determined at the core level:

λ C,k+1 = λ C,k w, F C Φ k+1 w, F C Φ k , (25)
where C stands for coarse, F for fission operator, Φ for scalar flux, w for a weight function and k is the power iteration index for the coarse core operator.

We turn now to a discussion of the changes to be introduced to convert the full core DDM calculation just outlined into a dynamic homogenization scheme. As mentioned in Section 1, the method of domain decomposition provides a natural paradigm to define a RHP, leading to homogenization techniques which can be applied to the local problem at each global iteration to construct a coarse operator that reproduces exactly the average transport quantities on coarse spatial and energy meshes. This homogenization is necessarily a dynamic process because the homogenized cross sections and the equivalence parameters used to construct the coarse operator depend on the transport solution and, therefore, as long as the transport has not converged, the coarse operator keeps changing.

With the DH method we invert the usual roles of the operators, where the coarse operator serves as an acceleration for the full core transport solution, and consider instead the low-order operator as the one that gives the solution of the core problem and the transport operator as the generator of homogenization parameters. This approach can open new possibilities for reducing the computational time of the global problem at the cost of introducing some approximations. For instance, the coarse operator does not necessarily have to reproduce the exact average transport quantities, and the local RHP that generates the homogenization parameters can be simplified with respect to the true situation in the core. A clear example of the latter is to adopt a 3D low-order operator and a 2D transport problem for the RHP.

At this point an issue comes up because each local transport problem is solved with a boundary source coming from its neighbors and with the imposed eigenvalue from the low-order operator, which now does not necessarily reproduce the eigenvalue of the full transport calculation. Moreover, because of the subdomain coupling, it is not possible anymore to simultaneously assure that 1) the core eigenvalue is the same everywhere in D and 2) the fine angular flux at the interfaces between subdomains is continuous. Therefore, to construct an iterative process where both the assembly homogenization and the global coarse operator converge, it is necessary to introduce a relaxation in the continuity condition of the angular flux by using a normalization factor that preserves the coarse incoming partial current:

ψ -,g i (r, Ω) | l+1 Γ ij = ψ +,g j (r, Ω) | l Γ ij × J -,C,G i (r) | l+1 core g∈G | n • Ω | ψ +,g j (r, Ω)dΩ | l Γ ij ∀g ∈ G, (26)
where l denotes the global iteration index, n the unit vector normal to the interface and g and G the fine (transport) and coarse energy groups. Note that the ψ + from the neighbor gives the fine distribution in space, energy and angle that allows to define the transport boundary conditions of each local RHP, normalized so as to preserve the coarse incoming partial current J -defined per macro-group and macro-surface. In this formula, we have used the same global iteration index for the incoming angular flux and the coarse incoming current to express that for each local RHP the transport and the coarse calculations are normalized by the same boundary source.

It is important to point out that the level of details employed for the boundary source in Eq. ( 26) is equivalent to that of a 2D direct transport calculation solved by DDM. This means that we are providing to each subdomain a very detailed information (ψ + ) but with an accuracy corrupted by the current normalization (J -). However, by introducing the discontinuity of the entering angular flux at the interfaces, we do not only ensure the convergence of the iteration process by preserving the eigenvalue of the low-order operator, but we also allow for the use of a simplified, rough distribution in space and angle of the incoming angular flux in Eq. [START_REF] Koebke | A New Approach to Homogenization and Group Condensation[END_REF]. The shape of the energy distribution as given by the angular flux exiting the adjacent subdomain is, on the other hand, mandatory.

At convergence, the current ratio in Eq. ( 26) is typically not equal to one but there exists two particular cases where it is. The first and less interesting one is when the coarse and the transport operators are exactly the same, so no homogenization arises. The second case is when a coarse operator is constructed so as to preserve the transport coarse partial currents at each coarse surface of the homogenized problem including the surfaces on the interfaces between assemblies. This entails that at the convergence of the iterative process the angular transport flux becomes continuous at the interfaces and, therefore, one obtains the exact full core transport solution. Clearly, in these circumstances the coarse operator acts as a nonlinear acceleration for the transport solution. In all other cases, the solution of the global problem depends on the low-order operator and the process is a truly dynamic homogenization method.

To ensure convergence of the DH method we check the following criteria between two successive global iterations l and l + 1:

• Eigenvalue λ C : |λ l C -λ l+1 C | < ε λ ; • Fission Integral F r per transport region r: |1 -F l r F l+1 r | < ε f ; • Incoming flux ψ -,g
s,d per surface element s, direction d and group g:

1 - ψ -,g,l s,d ψ -,g,l+1 s,d
2 < ε ψ , where ε λ , ε f and ε ψ are the tolerances accepted for convergence.

Finally, because the DH solution depends on the low-order operator and the transport information is locally available and can be used at negligible cost, it is expedient to introduce the advanced homogenization techniques which we presented in Section 2. We revisit next these techniques to analyze them in the context of the DH-like RHPs with non zero surface leakage.

Equivalence Theory in DH

With the presence of non-zero surface leakage in a reference homogenization problem, the preservation of the reaction rates per macro-group and macro-region is insufficient to ensure that the homogenized RHP will have the same eigenvalue as the reference transport calculation, since the eigenvalue also depends on the total leakage, which is not necessarily reproduced by the low-order operator. In order to preserve the reference balance, one must add an additional condition to the equivalence problem. The latter can be either the preservation of the reference net current at the boundaries per macro-surface and macro-group, which is only possible in diffusion theory, or the preservation of the reference eigenvalue with same incoming boundary source per macro-surface and macro-group, which is possible for any operator, let it be transport or diffusion like.

In the DH framework, the equivalence problem is naturally defined because the RHP is solved with imposed core eigenvalue and fixed incoming boundary source that is normalized so as to preserve the core coarse partial current. Hence, in the balance Eqs. ( 5) and ( 6) we have then J -= J -,C = J - core for each S ∈ ∂D and G, and λ = λ C = k core eff . The boundary source dictates the norm of both the transport and the homogenized problems, so if the method converges there exists a unique coarse operator that reproduces the reference balance for a given core state. It ensues that the homogenized problem will also reproduce the total outgoing current per macrogroup S∈∂D A S J + S , but not necessarily for each surface element S. Because of this, the fine angular flux is still discontinuous at the interfaces with neighbors, but the method offers the advantage of having the same local coarse power density for both transport and coarse problems, which can be used to perform a depletion calculation on the fine transport mesh.

Flux Discontinuity Ratios in DH

We remind that this method can only be used with coarse operators that are diffusion like. The basic idea here is to compute FDRs for all internal coarse surfaces of a subdomain. As for the equivalence homogenization, the RHP for the transport problem has an incoming boundary condition, but the RHP coarse problem is provided with reference net transport currents. Both problems have the same fixed eigenvalue. In these conditions, for the general case when the coarse operator has a submeshing the FDRs are computed from the iterative minimization of the functional in Eq. [START_REF] Tomatis | A multivariate representation of compressed pin-by-pin cross sections[END_REF], whereas in the case with no submeshing they can be obtained by solving a set of independent local problems, one per coarse region, with net currents boundary conditions. Regardless of the case, the computation requires the values of the net reference currents on all coarse surfaces and the final coarse core operator has discontinuity conditions at all internal coarse surfaces.

Recall that the boundary condition for the coarse operator which consists of taking the coarse net current equal to its transport counterpart allows to preserve the reference balance and therefore the reference reaction rates. However, as for any other local homogenization model, at convergence the transport currents are not necessarily those that would extant for the full core transport calculation. In this case, it would follow that transport and coarse full core solutions have different assembly power densities. The reason behind this is that the normalization in Eq. ( 26) makes the transport solution to be discontinuous at the subdomain interfaces.

We examine now two different implementations of the FDR technique. These implementations differ on the coarse conditions used at the interface between subdomains.

In the first approach we use continuity conditions for the flux of the low-order operator at the assembly interfaces and this leads to a typical DH method. Indeed, because of the normalization in Eq. ( 26), at convergence the final transport solution is discontinuous at the interfaces and so are then the transport flux and net current. This implies that the current used as boundary conditions for the local subdomains are not those of a full core transport calculation. Note that this also applies if we compute FDFs at the interface using either the GET or BB prescriptions in Eqs. [START_REF] Clarno | Capturing the Effects of Unlike Neighbors in Single-Assembly Calculations[END_REF] and [START_REF] Colameco | Iterative Transport-Diffusion Methodology for LWR Core Analysis[END_REF] and use them to construct FDRs to be used as assembly interface conditions for the full coarse operator.

In the second approach instead, we enforce the continuity of the transport flux at the assembly interfaces and use them to construct a fully coherent coarse core operator with FDR conditions at the interfaces between subdomains. Since the interface FDRs are computed using as coarse boundary conditions continuous transport currents, the iterations converge to the exact transport solution and the converged coarse operator reproduces exactly the coarse transport reaction rates. During iterations, and for some coarse operators also at convergence, the coarse incoming interface currents are different from the transport ones, so the incoming flux normalization in Eq. ( 26) acts as a nonlinear acceleration for the transport solution if the BB definition in Eq. ( 16) is used, since the latter preserves the average transport partial currents. When convergence is achieved, the coarse operator reproduces simultaneously the eigenvalue, the reactions rates and the partial currents of the fine transport operator on coarse energy and spatial meshes.

Let us explain how one enforces continuity of the angular flux at the interfaces while accelerating the incoming angular flux via Eq. [START_REF] Koebke | A New Approach to Homogenization and Group Condensation[END_REF]. Consider the interface between two subdomains A and B. At the l + 1 global iteration the transport solution for subdomain A is computed using the normalized incoming angular flux from the previous iteration, the ψ -,l+1 A given in Eq. ( 26), and produces an exiting angular flux ψ +,l+1 A . These are the angular fluxes which are used to compute the current boundary conditions for the coarse operator in the first approach discussed in the previous paragraph. Clearly, the transport flux is not continuous at the interface between A and B because ψ ±,l+1

A = ψ ∓,l+1

B

. In order to re-establish flux continuity, in the present implementation we proceed to sweep the incoming fluxes before constructing the coarse operator so that now ψ -,l+1

A = ψ +,l+1 B and ψ -,l+1 B = ψ +,l+1
A . Because in the parallel DDM implementation the sweep operation is done by an exchange of data between subdomains, the main difference between the two approaches based on the use of FDRs is that in the former the sweep is done after the construction of the coarse operator, while in the second it is prior to it. The second difference resides in that in the second approach one computes the FDRs at the interfaces.

20

The problem that we analyzed was inspired by the NEA PWR MOX/UO2 Core Transient Benchmark [START_REF] Kozlowski | PWR MOX/UO2 Core Transient Benchmark (Final Report)[END_REF], adopting some simplifications in order to have smaller data library to process for preliminary calculations. We considered only two types of assemblies instead of four, the UOx with 4.5% enrichment and the MOx with an average content of Pu-fissile of 4.3%. Both assemblies have been irradiated at critical infinite-lattice conditions to obtain the necessary macroscopic cross sections at three average burnup levels: fresh fuel at 0 GWd/t, once-burned at 20 GWd/t and twice-burned at 35 GWd/t. The core layout is shown in Fig. 3. The UOx assemblies contain guide tubes, fuel pins with and without IFBA (Integral Fuel Burnable Absorber). The MOx assemblies contain guide tubes with and without WABA rods (Wet Annular Burnable Absorber), and a three-zone MOx fuel with 5.0%, 3.0% and 2.5% fissile Pu (in respectively the interior, periphery and corner positions). Nominal power conditions have been considered to generate burnup dependent fuel and absorber compositions and the corresponding self-shielded cross sections assuming a unique average power density for each assembly. Three approaches are compared: Direct Transport (DT) considered as reference, two-step (DB2), and DH. To make a fair comparison, the same fine-transport operator was used for the three cases, and the same diffusion operator as coarse operator for the two-step and DH approaches. The twodimensional model problem in this core configuration is not critical. In order to avoid additional criticality searches in each of these calculations, via boron concentration or control rods, the reference multiplication factor is not equal to one. This does not compromise the validity or comparison of the presented methods as the DB2 leakage model, where employed, accounts for reference k eff value.

The same domain decomposition is also applied to the three approaches so that each subdomain corresponds to an assembly in the core including the reflector assemblies. The benchmark problem is symmetric with respect to the diagonal, so we computed one eighth of the reactor geometry with a total of 40 subdomains to be treated. We run all the calculations in a MPI parallel environment using 40 processors, so that each processor task concerns the computation and the homogenization of only one subdomain.

Our code organization was conceived so that both transport and coarse operators rely on the same spatial decomposition, in order keep all data only locally available in memory and, therefore, not accessible by the other processors. It ensues that the only information that is exchanged between processors involves the following quantities at the interfaces between subdomains: the angular flux for transport, the scalar flux and the homogenization parameters for CMFD and the partial currents for parabolic nodal diffusion (NEM2). In addition, the quartic nodal diffusion (NEM4) requires both partial currents and transverse leakage of the regions adjacent to the interfaces. In Appendix A we describe in details the Nodal Expansion Method that was adopted in our calculations and we also present in Appendix B how the method is used as nonlinear acceleration for transport.

One may notice that the few-group diffusion calculation could be performed efficiently without DDM. This is kept to preserve the consistency in the code organization, especially important for the future 3D calculations.

Direct Transport Approach

The standalone solver IDT of the code APOLLO3 R was used as transport solver, which is a multigroup discrete ordinates short characteristics solver [START_REF] Sanchez | APOLLO2 Year 2010[END_REF][START_REF] Zmijarevic | Multidimensional Discrete Ordinates Nodal and Characteristics Methods for the APOLLO2 Code[END_REF][START_REF] Masiello | New Numerical Solution with the Method of Short Characteristics for 2-D Heterogeneous Cartesian Cells in the APOLLO2 Code: Numerical Analysis and Tests[END_REF][START_REF] Ban | Code-to-code comparisons on spatial solution capabilities and performances between nTRACER and the standalone IDT solver of APOLLO3 R[END_REF]. The source iterations are carried out using the standard transport sweep, CMFD is used as nonlinear acceleration for the global iterations and the Boundary Projection Acceleration (BPA) for the inners. In our core application, the transport solver used a macroscopic self-shielded cross-section library in 172 energy groups with transport corrected P 0 scattering. A S 16 , level symmetric quadrature formula was adopted for the angular discretization and the linear characteristics expansion for the spatial representation of the angular flux. We divided each side of a pin-cell into 3 sub-surfaces for the projection on the polynomial basis of the surface flux. Note that the fission source is updated at each outer iteration using the accelerated flux moments.

A preliminary analysis was performed to find the optimal choice of coarse energy groups to be used for a pin-by-pin CMFD acceleration. We compared the overall calculation time using 172, 26, 6 and 2 group flux rebalancing according to Eqs. ( 22) and [START_REF] Sanchez | APOLLO2: a user oriented, portable, modular code for multigroup transport assembly calculations[END_REF]. Results showed that transport accelerated by two-group pin-by-pin CMFD had the lowest computational time, thus we considered this case to compare the performance with the DH approach. We refer to this reference case as DT-CMFD. This choice is in agreement with the analysis effectuated by Cho et al. in [START_REF] Cho | A Comparison of Coarse Mesh Rebalance (CMR) and Coarse Mesh Finite Difference (CMFD) Acceleration Methods for the Neutron Transport Calculations[END_REF] and Lenain et al. in [START_REF] Lenain | Coarse-Grained Parallelism For Full-Core Transport Calculations[END_REF]. The former showed that CMFD is effective when the optical thickness of the coarse mesh is smaller than one mean-free-path and it becomes unstable for the optically larger cells. The latter also showed that the collapsing in energy does not significantly affect the convergence properties of the iterative scheme, allowing actually to perform faster coarse calculations.

In a further analysis, we compared two DT calculations with same transport operator, but accelerated using Eqs. ( 22) and ( 23) by two different coarse operators: in the first case by two-group pin-by-pin CMFD and in the second case by two-group pin-by-pin Parabolic Nodal Diffusion with Flux Discontinuity Ratios (DT-NEM2+FDR). We observed that the number of outer iterations was the same in both calculations [START_REF] Rahnema | High-order cross-section homogenization method[END_REF] and that the run-times and the total number of inner iterations in transport were very close. Despite the fact that the computational cost of the two corrected diffusion operators is negligible with respect to transport, (up to 20 seconds over 58 minutes of whole run-time), CMFD is computationally less expensive then NEM2+FDR since the number of unknowns is lower. So we did not observed any particular advantage of NEM2+FDR over CMFD, but further studies will be done in future works.

In this work then, we will consider transport accelerated by two-group pin-by-pin CMFD as the best option for the reference calculation in our problem. Moreover, we also considered a coarse mesh homogenization for this low-order operator in order to fairly compare the convergence properties against other cases described later. The geometrical details for the coarse mesh homogenization will be presented with more details in the next paragraph and we refer to this case as DT-CMFD3x3, since the coarse mesh was constructed so that the assembly is homogenized on a threeby-three grid. Using this homogenization option, we also analyzed a calculation where transport is accelerated by a NEM4 operator with FDR according to Eqs. ( 22) and ( 23) (DT-NEM4+FDR).

The reference power per pin and per assembly is presented in Fig. 4 after normalization such that the average power per fuel pin or per fuel assembly is equal to one. The hot point can be found in a once-burned UOx assembly at the position B2 of the core (see Fig. 3). The distribution of relative power errors per pin or per assembly and the root mean square with respect to the reference are calculated using Eq. ( 27):

e i = P i -P i ref P i ref , RMS = N i i e 2 i N i , ( 27
)
where i is either the fuel pin, either assembly index, with N i being the total number.

DH Approach

For the DH calculations we used the same transport operator described in Section 5.1 based on the IDT solver and nodal diffusion as coarse operator. The transport solution is locally accelerated by BPA for the inners and by CMFD for the outers, where the latter uses the multigroup structure of transport. The number of local (per subdomain) outer iterations was fixed to one. Moreover, we decided to homogenize in energy from 172 to 2 groups in order to compare the results with the classical two-step scheme that typically applies two-group diffusion for the core calculation.

The diffusion coefficient in the fuel assemblies is equal to the leakage coefficient resulting from the homogeneous leakage model that corresponds to the imposed reactor eigenvalue, while in the reflector it was equal to 1 3Σtr .

The critical leakage model in the two-step calculation scheme is based on the fundamental mode assumption, and the experience has showed that the leakage coefficient used as diffusion coefficient produces more accurate results than 1 3Σtr [START_REF] Ragusa | Steady State Analysis of Multigroup Diffusion and SPN Methods in an MSLBlike Situation using the CRONOS Code System[END_REF]. The role of the leakage model in the classical approach is to adjust the heterogeneous spectrum so as to obtain the critical fundamental mode used as weighting function for cross-section homogenization. In order to obtain critical conditions, it is necessary to iterate between the equivalent homogeneous medium and the heterogeneous problem, because the former is constructed using an initial flux guess. On the other hand, the heterogeneous problem in DH has already realistic core conditions, so there is no need for flux adjustment. However, it is possible to create an equivalent homogeneous medium using the heterogeneous solution, and subsequently compute the homogeneous leakage rate that gives a multiplication constant equal to the given k core eff , where the latter is equal for both RHP and the core problem. The purpose of using a homogeneous leakage model in the DH calculation scheme is uniquely to compute the leakage coefficient and use it as diffusion coefficient. In a DH framework, the leakage D g B 2 ϕ g , with ϕ g the fine energy distribution of the equivalent homogeneous medium, is not necessarily equal for each energy group to the actual surface leakage computed in the heterogeneous problem. However, the energy-integrated leakage must correspond to the heterogeneous situation in order to preserve the balance and, therefore, the core eigenvalue. This may explain why the leakage coefficient contains more physical insight than 1 3Σtr and produces results of better quality. Our results corroborate this statement as well as the conclusion of [START_REF] Ragusa | Steady State Analysis of Multigroup Diffusion and SPN Methods in an MSLBlike Situation using the CRONOS Code System[END_REF].

We explored two types of spatial homogenization: the pin-by-pin and the coarse mesh.

With pin-by-pin homogenization we applied the NEM2 diffusion operator and we checked that the spatial mesh was fine enough and there was no need for higher order expansion. We preferred this operator rather then finite differences because we observed that the fixed point iterations used for the search of the equivalent cross sections showed much slower convergence rate than that of NEM2. With finite difference, it was possible to reach the 10 -4 pin-wise precision in reaction rates only by sub-meshing the pin-cells. This suggests that the NEM2 has a spatial representation of the within-cell flux much closer to the transport solution than the finite differences. Moreover, as said earlier, the computational cost of diffusion is negligible with respect to transport, so the choice of one of the two diffusion operators does not compromise the overall run-time.

We recall here that this type of spatial homogenization provides a detailed power distribution per pin, such that the impact of different homogenization options will be analyzed by comparing the pin power distributions against the reference. Here will be presented the diffusion solutions with the following options: a homogenization using only the flux-volume weighting (DH-NEM2), the one with equivalence (DH-NEM2+EQV), using flux discontinuity ratios at the interior surfaces (DH-NEM2+FDR), and also with the Black-Box discontinuity factors at the external surfaces (DH-NEM2+FDR+BB). The latter is the case when DH becomes a nonlinear acceleration and it differs from the cases discussed in Section 5.1 in two aspects: i) it does not apply the rebalancing Eqs. ( 22) and ( 23) but rather Eq. ( 26) only; ii) the fission source is accelerated locally using a 172-group pin-by-pin CMFD operator.

The coarse mesh homogenization is performed on a three-by-three spatial grid per assembly, where the assembly subdivision is done such that 17 pin-cells along the assembly axis is divided in groups of 5, 7 and 5 cells. In this type of homogenization we applied the NEM4 diffusion operator. For the comparison purposes, we are allowed to apply a straightforward reconstruction of the pin power distribution with the aid of the local transport solution, only if the coarse operator preserves the transport reaction rates in each coarse spatial mesh. This allows us to to analyze finely the homogenization options with equivalence (DH-NEM4+EQV) and with flux discontinuity factors applied to all the coarse surfaces (DH-NEM4+FDR+BB).

The stopping criteria of the iterative process introduced in Section 4 were set to: ε λ = 10 -5 , ε f = 10 -4 and ε ψ = 10 -4 . The same criteria were used for DT calculations.

Two-step Approach

The two-step calculation scheme is considered only in its pin-by-pin homogenization variant, since for the comparison, a coarse mesh homogenization would require a power reconstruction technique, whose study and comparison was not the objective of this work. We adopted the NEM2 diffusion operator for the reasons explained in Section 5.2. To homogenize the cross sections we applied the same transport operator as for the other approaches and, in order to have a consistent comparison, we forced the homogeneous critical leakage model to reproduce the reference core eigenvalue that is k eff = 0.97449. As in the DH approach, the diffusion coefficient in the fuel assemblies was equal to the leakage coefficient determined by the homogeneous leakage model, while in the reflector it was equal to 1 3Σtr .

One of the important points in the two-step scheme, as discussed in Section 3, is the reflector homogenization because it requires a separate calculation. In our work we investigated two different models for the reflector, and none of them uses cross sections weighted by the reference transport flux, because it would mean knowing a priori the solution, which is not always the case. The first model performs a transport calculation of a motif composed of one UOx assembly burned at 20 GWd/t next to a reflector assembly. Vacuum boundary conditions were imposed at the outermost side of the reflector. We produced the two-group cross sections for only one "output" medium mixing iron and water in the reflector. The homogenized data of this model were used for the whole reflector surrounding the core.

The second model considers four different motifs, representative of the periphery of the core as showed in Fig. 5, and this time we produced for each motif two "output" media, one for the iron and one for the water. The first two motifs are clusters of 2x1, where the fuel assembly is once-burned UOx and twice-burned MOx respectively for the first and second motif, and the reflector assembly has vacuum at the outer side. The other two motifs are clusters of 2x2 assemblies. In particular, motif 3 has one UOx assembly burned at 20 GWd/t and three reflector assemblies with vacuum boundary conditions on the outer sides. This motif is used to homogenize the part of the reflector on top-right. Motif 4, instead, has two once-burned UOx assemblies on top-left and bottom-right, one fresh MOx assembly on bottom-left, one reflector assembly on top-right, and conservative boundary conditions on all sides. Finally, for each part of the reflector surrounding the core we matched the cross sections homogenized using the most representative motif among the four. At this point, we wanted to analyze the impact of both reflector models on the power distribution.

To do that, for the two cases, where only the reflector cross sections differ, we used the fuel assembly cross sections obtained with the classical two-step approach applying the equivalence theory. We compared the power distribution of the two cases with the direct calculation, and observed that the one-medium reflector model exhibits higher pin power errors at the interface with the reflector, with a maximum of around 18% and a global root mean square of of 3.02%.

Using the four-motifs model instead, the maximum pin power error was around 6% and the global root mean square 1.91%. Fig. 6 and Fig. 7 show respectively the pin power error and the assembly power error for the two reflector models of the two-step scheme.

This difference of errors at the reflector interface was predictable, since one model utilizes more details than the other, and we also think that this position of the core may have a lower importance since the power is low. However, it is important to discuss that comparing the two cases that have same fuel cross sections and different reflector ones, the reflector modeling results in a different diffusion power distribution even in the center of the core. As the periphery assemblies are greater in number, the change in local power next to the reflector induce a significant change on power density in the center of the core if the total power of the reactor is maintained constant. In the UOx-20 assembly at position B2, which is the hot point of the core, the assembly power varies from 2.35% to 0.98%, while the power of the UOx-00 assembly at position B1 from 4.66% to 3.21%.

Clearly, the core eigenvalue problem redistributes the power so that a higher underestimation next to the reflector produces a higher overestimation in the center of the core. It follows that if the solution that we obtain is different for each reflector model that we adopt, the prediction of the hot points in the center of the core becomes more uncertain, which is an undesirable issue. For all the following comparisons, when referring to two-step, the reflector cross sections are those produced with the four-motif model. Like we did for DH, in order to show the impact of each homogenization option on the final solution we considered the following cases: leakage model only (DB2), leakage model with equivalence (DB2+EQV) and leakage model with Flux Discontinuity Factors (DB2+FDF). Note that for the latter we used the acronym FDF and not FDR because the discontinuities are introduced in all surfaces including the assembly interfaces. Moreover, for the reflector assemblies no FDFs were applied.

Figure 6: Effects of reflector modeling in the two-step scheme. Pin power errors.

DH vs DB2

Fig. 8 shows the relative pin power errors respectively for the two-step and DH approaches, applying flux-volume homogenization only, the equivalence and the flux discontinuity factors techniques. In Table 1 the maximum, minimum and RMS of the relative pin power errors are presented for each assembly in the core. The errors of the DH cases with Black-Box discontinuity factors at the interfaces between assemblies are not presented since the reference solution is exactly reproduced. First of all, for both of the approaches we can clearly see that the choice of a homogenization option that adjusts the low-order operator, such as EQV or FDF, does globally improve the solution of the two-group diffusion operator. The RMS of the relative pin power errors in the two-step approach is 2.2%, 1.91% and 1.7% respectively for the flux-volume homogenization, the equivalence, and the flux discontinuity factors options. In the DH approach with pin-by-pin homogenization, instead, the RMS is 1.9%, 0.94%, 0.65% and 0% respectively for the options flux-volume only, EQV, FDR and FDR+BB. The EQV and the FDF techniques result to be more effective in DH than in twostep thanks to the environment information. There is no guarantee, in fact, that an equivalence technique used in the two-step approach always improves the solution everywhere in the core, because the distribution of the reaction rates determined by solving the RHP can be very different from the actual one in the core, so forcing the low-order operator to reproduce it can be not always advantageous. An example of this can be seen for the fresh MOx assemblies at peripheral positions G2, G4 and F6, that show a peak of error in the central pins of the assembly for the case DB2+EQV, while for the case DB2 they are not as prominent. In other words, the choice of RHP may be inadequate, that is, at the homogenization stage the assembly solution was symmetric while in the core close to the reflector they experience strong gradients. More rigorous procedure would be to generate the homogenized cross sections for the peripheral assemblies using a fuel-reflector motif as RHP, like it was done for the reflector zones, but in general, while preparing the cross sections for the core calculation the position of a specific type of assembly is in general not known in advance. However, for local phenomena like the strong absorption in the fresh UOx assemblies containing IFBA, we can clearly observe that EQV and FDF properly adjust the diffusion solution in both two-step and DH approaches.

In the two-step scheme, some other error peaks are found because of the absence of the environment information in the model, such as the UOx-MOx interfaces, but also for the UOx-UOx interfaces with different burnup. In particular, the fresh UOx assemblies containing IFBA in positions B1, D1, E2 and D3, which belong to the interior part of the core, show a sensitivity to the environment influence with a maximum relative pin power error between 3% and 4%. These errors, however, result to be lower in the case of FDF option, probably because of the discontinuities at the assembly interfaces. On the other hand, in DH the errors introduced by the infinite lattice approximation at the interfaces between different types of assembly were eliminated and the maximum error remains lower than 2%, that typically arises in the corner pins of the assemblies in the inner part of the core.

The power distribution in the peripheral part of the core is predicted by two-step with a comparable quality using the three homogenization options, with peaks of errors between 4.5% and 6.5%, while in DH the behavior is different. In our DH calculations the reflector cross sections were homogenized into a 17x17 grid geometry, with the same mesh size as that of a fuel cell. We observed that in DH-NEM2+EQV the equivalence did not converge in the reflector subdomains, either due to the inefficiency of fixed point iterations in a problem with the strong flux gradients or to the weak values of the flux at the outer boundary. Therefore, the flux-volume weighted homogenized cross sections were used in the reflector coarse geometry instead of the equivalent cross sections. The case DH-NEM2+EQV shows an improved solution close to the reflector with respect to two-step, but this region still remains the most difficult to solve, since it strongly depends on the reflector properties that might not be well modeled by two-group diffusion, especially if no equivalence is applied. However, an enormous advantage of DH is the possibility to calculate directly the FDR also in this region of the reactor without applying any iterative process that may not converge, thanks to the fact that the geometry used for homogenization is locally available in each subdomain composing the reflector, while it is not possible in the two-step scheme. Finally, because of the correction of the diffusion operator also in reflector assemblies, the case DH-NEM2+FDR exhibits the best reflector response and, therefore, the best solution. This definitely shows how important the reflector homogenization is in reactor modeling.

In Fig. 9 the coarse mesh homogenization option is shown for the DH approach with equivalence, where the pin power has been reconstructed using the local transport flux of each assembly. Recall that the errors of DH-NEM4+FDR+BB are not shown since the reference solution is exactly reproduced. The coarse mesh homogenization option with equivalence shows a similar behavior of DH-NEM2+EQV in the interior part of the core, with slightly higher errors at the corners of the interfaces between assemblies, due to the coarser representation of the currents. On the other hand, the reflector response is not as well predicted as in the pin-by-pin. In this case the fixed point iterations for the search of the equivalent cross section converged also in the reflector subdomains but the EQV technique does not show any significant improvement in the quality of the solution at the peripheral area. This can be explained by two reasons: the first one is related to the type of operator, as the Fick's law fails to describe the physics of the problem in this region; the second one can be related to the choice of the mesh for cross-section homogenization, which is too coarse to represent the rapid spectrum change in the reflector.

DH-NEM4+EQV

Figure 9: Dynamic homogenization with equivalence and coarse mesh homogenization.

Relative errors [%] of the reconstructed pin power.

Convergence Rate and Run-times

Table 2 presents the eigenvalue error, the RMS of pin power errors, the number of global iterations and calculation run-time for each case that we investigated. All calculations were performed with 40 parallel process, where each process computes only one assembly.

As expected and explained in Section 5.1, CMFD constructed on coarse spatial mesh is less effective than pin-by-pin, resulting in a total of 18 global iterations instead of 12, and therefore a calculation run-time of 81 minutes instead of 58. This is due to two reasons: i) the average fission source is accelerated on a coarser mesh, so the transport operator needs more iterations for the pointwise convergence on the fine mesh; ii) the linear flux approximation of the finite difference scheme is insufficient for this mesh size to describe the flux gradients everywhere in the core, which ensues that the average flux is not well predicted and the acceleration becomes less effective. Therefore, we wanted to analyze the case DT-NEM4+FDR for the following reasons: the algorithm of the iterative scheme is identical for both calculations, the transport flux moments are accelerated using the same equations and same spatial and energy meshes, and the only difference is that the coarse flux is represented by a higher order of polynomial expansion which is sufficient for the size of the chosen coarse mesh. This acceleration shows a convergence rate of the transport fission source, which situates between the two CMFD operators, resulting in 13 global iterations with a run-time of 68 minutes. Figs. 10a and10b show the error decay as a function of the global iteration number, respectively for the transport fission source and the boundary angular fluxes.

In Fig. 10c we have plotted the number of core iterations needed by the coarse operator to converge as a function of the global iteration number, since at convergence it equals one. Note that: the global coarse problem is solved using the same DDM like in transport; as for the core iteration we refer to a power iteration solved by DDM, where at the end of each we update the coarse fission source, the eigenvalue and the coarse boundary conditions of each subdomain; our two-group diffusion operator are not accelerated.

At the beginning of the iterative process, we can see that the pin-by-pin CMFD for the first five global iterations reaches the maximum number of core iterations that was fixed to 500. A similar behavior is observed with the CMFD3x3 operator for the first two global iterations, while the NEM4+FDR operator exhibits a better convergence rate with less than 100 core iterations already after the first global iteration. On the other hand, during the last global iterations, where the solution is close to the convergence, we observe that the number of core iterations decreases slowly in the case of two coarse mesh operators with respect to the pin-by-pin operator, which entails that they are less effective in accelerating the pointwise convergence of the transport fission source.

The fact that pin-by-pin CMFD and coarse mesh NEM4+FDR have good properties in different phases of the iterative process may be of interest to develop a nonlinear acceleration with an adaptive mesh scheme, in order to take advantage of both operators and improve the convergence rate of the overall process. We suppose that this strategy may be convenient especially for 3D calculations, since the pin-by-pin CMFD is no more computationally negligible. In this work we only show the feasibility of using different operators as nonlinear acceleration and further analysis on the properties of discontinuous operators will be investigated in future.

We now analyze the cases DH-NEM2+FDR+BB and DH-NEM4+FDR+BB. First of all, the former converges in 11 global iterations while the latter in 13, resulting respectively in 60 minutes and 75 minutes of run-time. This can be explained, as said earlier, by the use of a coarser mesh to accelerate the transport convergence but, unlike the DT cases, here the coarse scalar flux is not used at all to accelerate, neither the fission source nor the outgoing angular flux, but rather the coarse partial currents to accelerate the incoming boundary source. In Fig. 10b we observe that the boundary source of DH-NEM2+FDR+BB converges faster than all the other cases, in part because of the use of the Eq. ( 26) instead of Eq. ( 23), since the former has more physical insight, and in part because of the multigroup local acceleration. When comparing DH-NEM2+FDR+BB with DT-CMFD and DH-NEM4+FDR+BB with DT-NEM4+FDR, we can clearly see in Fig. 10a that the convergence rate of the fission source does not significantly change. However, the fission source is accelerated differently, either using the local information (DH) or the global information from the diffusion solution (DT). This is explained by the fact that fission is an internal source while the eigenvalue and the boundary flux depend on the external environment, and in both cases (DT and DH) the latter information is accelerated using coarse quantities that somehow contain the environment information. Nevertheless, we can observe a difference of the two approaches for non-fissile subdomains. In Fig. 10d we have plotted the total number of transport inner iterations, where total stand for the sum of all subdomains, fissile and not. For the first five global iterations the transport operator in all the cases reaches the maximum number of inner iterations which was fixed to three in fissile subdomains and 50 in non-fissile subdomains, since in the latter there is no outer loop on fission source. We can observe that DH-NEM2+FDR+BB requires more inners than DT-CMFD and DH-NEM4+FDR+BB more than DT-NEM4+FDR. This is explained by the lack of using Eq. ( 22) in the DH cases, since it gives the advantage of accelerating the convergence of the multigroup problem also in non-fissile subdomains. We checked, in fact, that the increase of the number of inner iterations occurs mainly in the reflector subdomains. It follows that the overall run-time of DH-NEM2+FDR+BB and DH-NEM4+FDR+BB is slightly higher than the respective DT cases with same coarse mesh homogenization option.

All of Dynamic Homogenization cases show on average a run-time that is comparable to Direct Transport accelerated by pin-by-pin CMFD. The fastest cases are DH-NEM2 and DH-NEM2+EQV mainly due to the lesser number of global iterations (10 instead of 12) with a run-time of 50 minutes for both, since EQV iterations have a negligible cost. We have observed that in these cases the diffusion solution after 5 iterations does not change significantly, which means that cross sections are close to the converged values, but the process needs 5 more iterations to converge the transport quantities. We think that the slightly lower number of global iterations is due to the fact that when using a low-order operator that preserves net currents, even if the homogenized cross sections do not change very much, its matrix coefficients also depend on the homogenized transport currents and surface fluxes per coarse region, adding then additional dependencies from transport. This is why the convergence rate of the cases DH-NEM2+FDR, DH-NEM2+FDR+BB and DT-CMFD cannot be very different.

Comparing the cases DH-NEM2+EQV and DH-NEM4+EQV we observe that they have same eigenvalue, even if the power distribution of the former and the reconstructed pin power of the latter are different with the RMS of relative errors being respectively 0.94% and 1.47%. As we observed for the direct calculations, the coarse mesh homogenization option exhibits a higher number of inner iterations in transport as showed in Fig. 10d. Therefore, even if the coarse diffusion calculation (NEM4+EQV) is computationally less expensive than the pin-by-pin diffusion (NEM2+EQV) calculation (respectively 3 seconds and 10 seconds cumulated at convergence), the overall run-time of the latter is lower because almost all of the computing time is spent for transport.

Finally, the two-step calculations, also shown in Table 2, are the most rapid and the least accurate as expected. The run-time of DB2 cases includes only the diffusion calculation. For the interested reader, the initialization of DH and DT cases consists of infinite lattice calculations that produce the first weighting flux for cross-section homogenization and, therefore, the zero-th global iteration can be considered as a two-step calculation. The initialization has a computational cost of roughly 10 minutes and it is, evidently, included in the total run-times. In our two-step cases, the lattice calculation takes 12 minutes instead of 10 because of the critical buckling search and the critical fundamental mode that is not performed in DH and DT. However, it is hard to include this additional cost of generating the homogenized cross sections in the overall two-step computational cost for several reasons: i) one eighth of assembly is typically computed to take advantage of the geometrical symmetries, while in our cases for simplicity of implementation we computed the whole assembly also for two-step; ii) the separate calculation for reflector homogenization is not included in the 12 minutes discussed earlier, and in our two-step case the code read an external library for only the reflector cross sections. The highest computational cost for generating the reflector library among the 4 motifs was 25 minutes associated to motif 3 (see Fig. 5), that was solved using a 2x2 domain decomposition, 4 parallel process, same transport operator accelerated by two-group pin-by-pin CMFD.

In summary, the results show that the DH approach does not exhibit any considerable advantage for 2D configurations in terms of computational cost with respect to a two-step approach, which is obvious, nor to a direct transport with a non-linear acceleration. We also investigated the possibility of simplifying the information transmitted at the boundaries so as to relax the constraints of the iterative process, with the purpose of reducing the number of global iterations and, therefore, the computational cost. We introduced the approximations to the boundary conditions by reducing the level of details of the fine distribution in angle and space of the incoming angular flux in Eq. [START_REF] Koebke | A New Approach to Homogenization and Group Condensation[END_REF].

To do that, we considered the case DH-NEM2+EQV, since it showed the lowest number of globals, and imposed an isotropic boundary flux averaged in space per each coarse surface. With the solver options that we used, this results in one value per pin-cell side and per group at the assembly interfaces that is transmitted to the neighbors instead of 432 values per cell side and per group. The latter accounts for 72 directions, 3 sub-surfaces of flux projection over a linear polynomial basis (2 surface moments). This calculation exhibited a RMS of the relative pin power errors of 0.97% instead of 0.94% for the case without approximations at the boundaries, which means that in our configuration the latter do not drastically impact the accuracy of the solution. The total number of inner iterations was slightly reduced, while the number of global iterations did not change at all (10 iterations). We also observed that the diffusion solutions of the cases DH-NEM2 and DH-NEM2+EQV did not change significantly after five global iterations, which entails that the homogenized cross sections are close to the converging values, but the transport operator needs more iterations to converge the local multigroup problem. We suppose then that the convergence rate of the iterative process is more sensitive to the convergence of the local fission source rather then the convergence of the incoming source. spect to a direct calculation accelerated by two-group pin-by-pin CMFD. This means also that the calculation run-times are very close. The biggest but not appreciable difference is 10 global iterations instead of 12, which results in 50 minutes instead of 58, associated to the DH cases with flux-volume homogenization and equivalence options. Therefore, using the same computational resources of a direct transport calculation, the DH method does not show considerable advantages for 2D configurations in terms of performances.

Case

However, the objective of this work was to investigate, first of all, that the solution of the DH approach has an excellent agreement in the radial plane with the reference calculation. Moreover, results demonstrate that the choice of a two-group diffusion operator with appropriate homogenization parameters is largely valid for core modeling and design even for high heterogeneous configurations. We do think that improvements on the performance are promising for the 3D core problem via DH where no 3D transport computation is needed.

For equal homogenization options, DH always shows a better quality of the solution with respect to the classical two-step approach, mainly due to the elimination of the typical error peaks at the interface between assemblies of different type and to the improvement of the reflector response, where two-step easily fails. The RMS of the DH diffusion relative pin power errors is, in fact, 1.9%, 0.94% and 0.65% respectively with flux-volume homogenization option, with equivalence and with flux discontinuity ratios in the interior surfaces of the assembly. The RMS of the diffusion pin power errors for the two-step approach is instead 2.2%, 1.9% and 1.7% respectively for the same aforementioned homogenization options. Note that the pin power errors are calculated with respect to a direct transport calculation considered as reference.

We also explored a coarse mesh homogenization, where the assembly is homogenized in a 3-by-3 grid. Especially for this type of homogenization, the transport flux that is locally available can be used to reconstruct the power within the pin. A straightforward reconstruction requires the average transport power and the diffusion power to be equal for each coarse mesh. This is possible only if the equivalence or the flux discontinuity ratios in all surfaces are applied. The latter is a a nonlinear acceleration, while in the former the RMS of relative errors of the reconstructed power is 1.47%. The higher errors of this case with respect to the pin-by-pin homogenization arise mainly at the reflector interfaces, which entails these this homogenization options (coarse mesh with equivalence) may not be sufficient to model the reflector response and the inadequacy comes from the nature of the coarse operator.

When the discontinuities are introduced also at the assembly interfaces so as to preserve the transport net currents in each macro-surface, the diffusion operator reproduces exactly the transport reaction rates, which entails that it can be used to accelerate the convergence of the fission and boundary transport sources. A comparison between CMFD and parabolic nodal diffusion operator with two-group pin-by-pin homogenization showed comparable transport convergence properties in our configuration. However, the two-group CMFD constructed on coarse spatial mesh is less effective than a quartic nodal diffusion operator with same homogenization options, with a total number of global iterations of 18 and 13 respectively. This is due to the linear flux approximation of the finite difference scheme that might be insufficient to describe the flux gradients everywhere in the core for the mesh size that we considered (an assembly homogenized on a 3-by-3 grid), which ensues that the average flux is not as well predicted as the higher order flux expansion of the nodal operator.

An interesting technical aspect of our work is that the same code organization based on domain decomposition allows to perform calculations with all the three approaches: two-step, dynamic homogenization and direct transport. This means that the DH or the direct iteration strategies can become a solver option of the same code with no significant modifications.

Homogenization with net boundary leakage poses new problems as to the use of homogenization methods and particular care has to be taken in order to preserve the total net assembly leakage.. Therefore, we have analyzed and adapted the equivalence theory and the flux discontinuity ratios technique in the framework of non-conservative boundary conditions . Another aspect of this problem is the cross section self shielding calculation with, to be consistent with the DH homogenization technique, will have to be done with non conservative boundary conditions. In this work, we did not address the question of self-shielding, which should be performed at each depletion step. One has to be aware that many traditional self-shielding methods rely on a heterogeneoushomogeneous equivalence that relates a real problem to the reaction rates of an infinite homogeneous medium. Treating the assembly as an open domain with incident surface sources brings up the need for a self-shielding technique that allows for a net current at the boundaries, a condition that eliminates the direct use of most traditional methods for multigroup self shielding. However, a good candidate is the conventional sub group formulation [START_REF] Jung | Practical numerical reactor employing direct whole core neutron transport and subchannel thermal/hydraulic solvers[END_REF], where an escape probability is defined which helps to represent the net leakage as the difference of an escape term, proportional to the local flux, minus a source term which is proportional to the infinite medium flux. This technique is used in the 2D-1D fusion calculation scheme, which perform self shielding calculations over the whole radial plane eliminating the fundamental-mode assumption. This is very similar to what we do in DH except that our "transport" solution is discontinuous at assembly interfaces, while preserving the shape of the neutron spectrum. Such a work, which will be a generalization of the 2D-1D sellf shielding model to a heterogeneous region, will be done in future work, allowing thus the evaluation of the error introduced by our present use of traditional schemes for multigroup self shielding.

In future works, we will extend the method to 3D configurations and will test the core depletion calculation in order to quantify the bias introduced by the two-step approach to the assembly isotopic content.

with:

Φ x (x) = 1 ∆y∆z ∆y ∆z dydzΦ(x, y, z), L x (x) = L xy (x) + L xz (x) = 1 ∆y∆z ∆y ∆z dydz( d dy J y (x, y, z) + d dz J z (x, y, z)), q x (x) = HΦ x (x) + 1 λ F Φ x (x),
where Σ r is the removal cross section, H and F are the operators respectively for the scattering and the fission sources and λ the eigenvalue. Note that the group index has been omitted. The one-dimensional flux Φ x (x) is then approximated using an expansion of Legendre polynomials truncated at the N -th order and defined on a node k of size ∆x:

Φ x (x) = Φ + N n=1 φ x,n P n (µ), (29)
where

x ∈ [- ∆x 2 , ∆x 2 ], µ = 2x ∆x , µ ∈ [-1, 1].
From the Lengendre polynomials properties we have:

P n (±1) = (±) n ,
so the expression of the flux at the interfaces x ± of a node k, where + andstand respectively for right and left, is:

Φ x ± = Φ + N n=1 (±) n φ x,n , (30)
while knowing that:

J x (x) = -D∂ x φ(x), ∂ x P (µ) = 2 ∆x P (µ), P n (±1) = (±) n+1 n(n + 1) 2 ,
we can derive an expression for the current at the node interfaces:

J x ± = - 2D ∆x N n=1 (±) n+1 n(n + 1)φ x,n , (31)
A.1 Parabolic Expansion (NEM2)

We now want to eliminate the dependency of the flux moments in favor of the response matrix formulation for the partial currents exiting (J out ) and entering (J in ) a node. Substituting the following equations

Φ x ± = 2(J out x ± + J in x ± ), (32)
into Eq. ( 30), we obtain an expression for the flux moments as a function of the partial currents:

φ x,1 = 1 2 (Φ x + -Φ x -) = (J out x + + J in x + ) -(J out x -+ J in x -), φ x,2 = 1 2 (Φ x + + Φ x -) -Φ = (J out x + + J in x + ) + (J out x -+ J in x -) -Φ.
Replacing these expressions in Eq. ( 31), it follows:

J x + = J out x + -J in x + = - D ∆x [8(J out x + + J in x + ) + 4(J out x -+ J in x -) -6Φ], (33a)
J x -= -J out x -+ J in x -= D ∆x [4(J out x + + J in x + ) + 8(J out x -+ J in x -) -6Φ]. (33b)
Now, by averaging Eq. ( 28) over a node k, the discretized balance equation becomes:

1 ∆x [J x + -J x -] + Σ r Φ = Q x , (34)
with

Q x = 1 ∆x ∆x dx[q x (x) -L x (x)] = q -L x , and
L x = L xy + L xz = 1 ∆y [J y + -J y -] + 1 ∆z [J z + -J z -].
Finally, from the balance Eq. ( 34) we explicit the average scalar flux in a node:

Φ = 1 Σ r [Q x - 1 ∆x (J out x + -J in x + + J out x --J in x -)], (35)
and substitute it into the definitions of the interface net currents (Eq. ( 33)) to obtain:

A x J out x = B x J in x + Q x , ∀k ∈ [1, N x ], (36)
where

J out x = J out x - J out x + , J in x = J in x - J in x + , Q x = γQ x γQ x , A x = 1 + 8α + 6β 4α + 6β 4α + 6β 1 + 8α + 6β , B x = 1 -8α + 6β -4α + 6β -4α + 6β 1 -8α + 6β , α = D ∆x , β = D Σ r (∆x) 2 , γ = 6D Σ r ∆x .
Eq. ( 36) is a system of 2N x equations, with N x the number of nodes for an axis x, that is solved for the outgoing currents by inverting a tridiagonal matrix. Note that for the continuity condition of the partial current, the incoming current is equal to the current exiting the neighboring nodes, therefore:

J in x = J out x nk and J out x nk = J out x k-1 + J out x k+1 - , (37)
where nk stands for neighbors of node k. The 3D solution is achieved by sweeping along the three axes and iterating on the transverse leakage sources L x , L y , L z , that are updated at the end of each sweep (inner iterations loop). When L x , L y , L z converge, the new average flux is determined with Eq. ( 35) using the last scattering and fission sources that are fixed for the whole inner iterations loop. The diffusion NEM2 operator can admit boundary conditions either for the net currents, or for the incoming partial currents, or also for the boundary flux. For specular boundary conditions, clearly, zero net current has to be imposed, while for vacuum boundary conditions typically the incoming current is set equal to zero. When the NEM2 operator is spatially decomposed, the only information required by each subdomain to perform the sweep is the partial currents exiting the neighbors according to Eq. ( 36), somehow similarly to transport that involves only the boundary angular fluxes and no information of the neighboring geometry is needed.

A.2 Quartic Expansion (NEM4)

The interest of a higher order flux expansion arises when the geometry contains large coarse meshes, because it avoids the increase in the number of computational meshes used to solve the discretized problem. In this case the two conditions of Eq. ( 32) are not sufficient to eliminate the dependency of the four flux moments in favor of the partial currents. In order to find the two additional constraints we can take advantage of the Legendre polynomials properties. We substitute Eqs. ( 29) and ( 31) into Eq. ( 28) and we write the balance equation for a node k:

- 4D (∆x) 2 N n=1 φ x,n P n (µ) + Σ r (Φ + N n=1 φ x,n P n (µ)) = q x + N n=1 q x,n P n (µ) -L x , (38)
where we have assumed that the scattering and fission sources have same expansion of the flux, and the diffusion coefficient is constant within the node. We recall here that the orthogonality of the Legendre polynomials yields:

< P m , P n >= 1 ∆x ∆x/2 -∆x/2 P m (µ)P n (µ)dx = 1 2 1 -1 P m (µ)P n (µ)dµ = 1 2m + 1 δ m,n ,
where δ m,n is the Kronecker function, and we observed that:

< P m , P n >= (2m + 3)δ m+2,n .

By projecting Eq. ( 38) on P m and using the latter properties we obtain:

-4D (∆x) 2 (2m + 3)φ x,m+2 + Σ r 1 2m + 1 φ x,m = 1 2m + 1 q x,m -< P m , L x >, which allows to have an expression for the higher moments as a function of the lower order moments of the flux and the source :

φ x,m+2 = α m Σ r φ m -Q x,m , (39)
where α m = 1 4(2m + 1)(2m + 3) (∆x) 2 D and Q x,m = α m (q x,m -L x,m ), L x,m = (2m + 1) < P m , L x >, m > 0.

We now substitute the expressions for the moments of order 3 and 4 given by Eq. ( 39) into Eqs. [START_REF] Aragonés | A Linear Discontinuous Finite Difference Formulation for Synthetic Coarse-Mesh Few-Group Diffusion Calculations[END_REF] and [START_REF] Smith | Full Core, 2-D, LWR Core Calculations with CASMO-4E[END_REF]:

Φ x ± = Φ ± (1 + α 1 Σ r φ x,1 ) + (1 + α 2 Σ r φ x,2 ) ∓ Q x,1 -Q x,2 , (40)
J x ± = - 2D ∆x [(1 + 6α 1 Σ r φ x,1 ) ± (3 + 10α 2 Σ r φ x,2 ) -6Q x,1 ∓ 10Q x,2 , (41)
and by combining Eqs. ( 32) and ( 40) we obtain the expressions for the first two flux moments as a function of the partial currents:

φ k x,1 = 1 1 + α 1 Σ r [(J out x + + J in x + ) -(J out x -+ J in x -) + Q x,1 ], φ x,2 = 1 1 + α 2 Σ r [(J out x + + J in x + ) + (J out x -+ J in x -) + Q x,2 -Φ].
We then substitute the latter expressions and Eq. ( 35) into Eq. ( 41) to obtain:

J x ± = ±(J out x ± -J in x ± ) = - D ∆x {c 1 [(J out x + + J in x + ) -(J out x -+ J in x -)] ± c 2 [(J out x + + J in x + ) + (J out x -+ J in x -)] ± c 2 Σ r ∆x [(J out x + -J in x + ) + (J out x --J in x -)] ∓ c 2 Σ r Q x -(12 -c 1 )Q x,1 ∓ (20 -c 2 )Q x,2 },
where

c 1 = 2 1 + 6α 1 Σ r 1 + α 1 Σ r , c 2 = 2 3 + 10α 2 Σ r 1 + α 2 Σ r , Q x = (q -L x ).
The response matrix formulation for one node k is then:

A x J out x = B x J in x + Q x , ∀k ∈ [1, N x ], (42)
where

J out x = J out x - J out x + , J in x = J in x - J in x + , Q x = Q x,- Q x,+ ,
where the source terms are:

Q x,± = D ∆x [ c 2 Σ r Q x ± (12 -c 1 )Q x,1 + (20 -c 2 )Q x,2 ],
and the matrix coefficients are: At this point we need to find the closure equation that allows to determine the moments of the transverse leakage source. We recall that the latter has two components corresponding to each transverse axis: L x (x) = L xy (x) + L xz (x), and both are functions of x. The projection on then polynomial basis of the latter can be determined if the shape of the functions is known. We introduce then a model that approximates the leakage shape assuming: i) a quadratic expansion over a node k of the following form: J k y + -J k y - ∆y k , in order to take advantage of the orthogonality of the Legendre polynomials in the projected equation; it ensues that in order to define the parabolic function, two additional conditions are needed and the latter are given by the assumption that ii) the average of this function over the adjacent nodes results in the average transverse leakage: xy , where bbd stands for beyond the boundaries and k bd for boundary node. For vacuum boundary conditions the leakage bbd is typically imposed equal to zero. When the NEM4 operator is spatially decomposed, each subdomain requires the exchange with its neighbors of the incoming partial currents, the average leakage at the boundary nodes for each transverse axis and the spatial mesh size of the latter according to Eqs. ( 42) and (43).

Appendix B Discontinuous NEM and acceleration (NEM+FDR)

When the Flux Discontinuity Factors technique is applied, both scalar flux and partial currents are no more continuous and at the interface of two nodes k and k -1 it yields:

f x k-1 + Φ x k-1 + = f x k -Φ x k -⇒ f x k-1 + (J out x k-1 + + J in x k-1 + ) = f x k -(J out x k - + J in x k - ), (44)
while the net current is continuous:

(J out x k-1 + -J in x k-1 + ) = -(J out x k - -J in x k - ). (45)
Using Eqs. ( 44) and (45) we can express the incoming partial current as a function of the outgoing currents of the neighboring nodes and that of the node itself:

J in x k - = 2f x k-1 + f x k-1 + + f x k - J out x k-1 + + f x k-1 + -f x k - f x k-1 + + f x k - J out x k - .
In practice, only the Flux Discontinuity Ratios

r x k ± = f x k±1 ∓ f x k
± are used to solve the system, and the last equation written for both right and left surfaces of the node becomes:

J in x k ± = 2r x k ± r x k ± + 1 J out x k±1 ∓ + r x k ± -1 r x k ± + 1 J out x k ± . (46)
Clearly, if the two discontinuity factors in an interface are the same, there is no discontinuity so the ratio equals one, the second term of Eq. ( 46) vanishes and we find again the continuity condition of the partial current. We rewrite Eqs. ( 36) and ( 42) so as to take into account the discontinuities:

(A x -B x • F1 x )J out x = B x • F2 x J out x nk + Q x , ∀k ∈ [1, N x ], (47)
where

F1 x =   rx --1 rx -+1 0 0 rx + -1 rx + +1   , F2 x =   2rx - rx -+1 0 0 2rx + rx + +1   , ∀k ∈ [1, N x ].
At this point we need to determine the FDRs to feed the operators F1 x and F2 x in Eq. (47). When submeshing is not applied to a zone of homogenization, the FDFs can be calculated directly and a priori according to Eq. ( 13) or Eq. ( 16) without employing any iterative process. In the latter equations, the numerators are transport quantities, which are available at the homogenization stage, while one needs an expression for the coarse surface flux Φ C S at the denominator as a function of the average transport flux in the node and the average transport net current at the surfaces, since both are the quantities to be preserved. Note that for a direct computation of the FDR, only the knowledge of coarse surface flux Φ C S is required.

For the NEM2 operator this expression can be determined by combining Eqs. ( 32) and ( 33):

Φ C=NEM2 S ± = Φ R ∓ ∆x R 6D R (2J S ± + J S ∓ ). (48)
In Eq. (48) we adopted the notation of Section 2 on the homogenization theory, where S ± are the right (+) and the left (-) surfaces of a coarse region R, ∆x R is its size along an axis x, Φ R and J S ± are respectively the average transport scalar flux in R and the transport net currents across S ± , and the diffusion coefficient D R can be arbitrary chosen. Note that for a given D R a set of six FDFs is completely defined (one for each surface in a 3D mesh) and when it is used in Eq. (47), the transport solution is exactly reproduced in average.

For the NEM4 operator an approach to determine Φ C S ± is to use Eq. ( 30) that requires the knowledge of average transport scalar flux and the values of the four flux moments. The former is one of the quantities to be preserved, while the latter can be computed using the following four conditions: preservation of the net currents on the right and left surfaces of a node along a direction; preservation of the average transverse leakage of the right and left adjacent nodes. Because of Eq. ( 39), the flux moments are coupled by the scattering and fission sources, which entails that a multigroup system of equations has to be solved:

A G x Φ G x = Q G x + J G x , ∀G, (49)
where

A G x =     2 6 12 20 2 -6 12 -20 Σ r 0 60d 0 0 Σ r 0 140d     G R , Φ G x =     φ x,1 φ x,2 φ x,3 φ x,4     G R , Q G x =     0 0 q x,1 q x,2     G R , J G x =     d∆xJ x + d∆xJ x - -L x,1 -L x,2     G R
, where

d = - D (∆x) 2 , q G x,n = G =G Σ G →G s0 φ G x,n + χ G λ G νΣ G f φ G x,n .
In Eq. ( 49) the indexes R and G mean that the vector and matrix elements are evaluated for the macro-group G and coarse region R. The first two equations of the system are obtained from Eq. ( 31), while the last two equations from Eq. ( 39). The moments of the leakage source L x,1 and L x,2 are calculated according to Eq. ( 43), where the average leakage of the adjacent nodes are determined with the transport net currents in the transverse direction, while to evalualate the source moments q G x,n the reference eignevalue has to be imposed, adding an additional condition to the equivalence problem. Note that the classical notation was used for the cross sections and the fission spectrum. The system of Eq. ( 49) is composed of 4N G equations, with N G the number of macro-groups, and it can be solved by direct inversion or by iterating on the source Q x . Once the flux moments are computed, the coarse surface flux for the NEM4 operator is detemined with:

Φ C=NEM4 S ± = Φ R + 4 n=1 (±) n φ R,x,n , (50)
where we have omitted the group index.

Although a diffusion NEM operator can admit boundary conditions either for the net currents, or for the partial currents, or for the boundary flux, we recall here that when a RHP has non-zero surface leakage it is not possible to determine the FDF beyond the boundaries so, in order to reproduce the transport balance, one has to impose the preservation of the transport boundary net currents.

When Eq. ( 47) is used as nonlinear acceleration, the matrix coefficients depend on the transport flux and net currents, thus they change at each iteration:

(A x -B x • F1 x ) l+1/2 J out,l+1 x = (B x • F2 x ) l+1/2 J out,l+1
x nk

+ Q l+1 x , ∀k ∈ [1, N x ], and
Φ l+1 = 1 Σ l+1/2 r [Q l+1 x - 1 ∆x (J l+1 x + -J l+1 x -)].

Figure 1 :

 1 Figure 1: Two-step calculation scheme.

Figure 2 :

 2 Figure 2: Dynamic Homogenization calculation scheme.

Figure 3 :

 3 Figure 3: Core layout of NEA Benchmark [9]. It contains two types of assembly (UOx, MOx) at three burnup (0, 20, and 35 GWd/t).

Figure 4 :

 4 Figure 4: Normalized reference power distribution.

Figure 5 :

 5 Figure 5: Reflector homogenization. Four-motif model.

Figure 7 :

 7 Figure 7: Effects of reflector modeling in the two-step scheme. Assembly power errors.

Figure 8 :

 8 Figure 8: Relative pin power errors.

  Core Iterations per Global Iteration (d) Total Number of Transport Inner Iterations

Figure 10 :

 10 Figure 10: Convergence Properties

A x = 1 +

 1 a + e -b + e -b + e 1 + a + e , B x = 1 -a + e b + e b + e 1 -a + e , with a = D ∆x (c 1 + c 2 ), b = D ∆x (c 1 -c 2 ), e = Dc 2 Σ r (∆x) 2 .

) 2 + 3α k±1 + 2 2 .

 22 L k xy,n , where, by adopting the change of variable x = ±∆x 2 µ , which implies the change µ → ±µ , the expansion coefficients are: k±1 ∆x k . It follows that the moments of the transverse leakage source can be determinated using: Despite the boundary conditions either for the net currents, or for the incoming partial currents, or for the boundary flux, the NEM4 operator requires an additional condition for the transverse leakage source at the boundary nodes, since the latter depend on the leakage of the adjacent nodes. For specular boundary conditions one has to impose L bbd xy = L k bd

Table 1 :

 1

	max e i					0.	1.24							0.	3.94
	min e i					6.24 4.3							4.38 4.8
	RMS					3.99 3.14							2.78 2.95
					0.	1.62	0.						0.	1.83	1.77
					2.65	4.55	4.13						2.79	4.94	4.15
					1.39	2.56	3.23						1.36	2.29	2.56
				1.89	2.94	0.	0.	1.01				1.42	2.79	0.	0.	4.31
				2.09	3.03	3.29	5.45	4.27				1.52	2.72	3.43	3.49	4.8
				1.12	2.04	1.95	3.05	2.96				0.62	1.9	1.7	1.83	2.66
			2.15	4.53	0.55	1.62	2.51	1.6			2.32	4.28	0.	1.77	3.18	3.32
			0.41	1.71	1.72	3.13	3.02	3.83			0.71	1.97	1.81	1.48	3.86	2.92
			1.25	2.44	0.69	1.25	1.9	1.78			0.99	2.14	0.83	0.64	2.16	1.38
		2.16	2.93	3.34	4.5	1.5	0.57	0.		2.33	2.72	3.04	4.35	1.44	0.7	2.98
		0.	0.	1.47	1.37	0.48	4.64	2.86		0.	0.	0.48	1.92	1.54	2.44	3.13
		1.63	2.28	1.98	2.51	0.62	2.07	2.12		1.29	1.78	1.35	2.23	0.72	1.06	1.89
	2.69 5.4	2.01	5.09	1.97	2.38	3.19	0.39	2.2	5.03 1.85	4.82	1.67	2.3	3.56	2.67
	0.	0.	0.	0.36	0.	2.31	2.8	1.97	0.	0.18	0.	1.27	0.22	0.58	3.48	2.32
	1.98	3.25	1.51	3.02	1.28	1.4	1.93	1.24	1.41	2.77	1.18	2.59	0.91	1.07	2.14	1.22
					DB 2							DH-NEM2
						0.89	2.46							1.28 3.24
						5.46	3.3							1.91	3.04
						2.62	2.1							1.	1.72
					1.75	0.95	0.						1.66	0.91	1.09
					2.05	3.15	4.02						0.	1.15	3.43
					0.72	0.78	2.92						1.01	0.35	1.58
				4.13	3.13	0.	0.38	2.2				2.23	1.79	0.92	0.	2.11
				0.07	1.14	5.23 5.82 3.29				0.	0.45	1.46	2.37 3.61
				2.11	1.84	2.9	3.	2.19				1.28	0.74	0.79	1.55 2.05
			0.99	2.94	0.8	0.62	0.13 6.4			0.88	1.39	1.03	0.25	0.7	1.48
			0.09	0.	1.65	2.91	3.48	2.93			0.17	0.	0.59	0.97	0.89	1.24
			0.57	2.41	0.75	1.26	0.86	2.17			0.13	0.69	0.32	0.41	0.3	0.71
		2.43	1.59	2.67	3.07	0.28	0.36	0.		1.6	0.19	2.03	1.59	1.11	0.05	1.42
		0.	0.	0.	0.59	3.63	4.94	3.31		0.04	0.2	0.	0.	0.57	1.19	2.08
		0.98	1.26	1.45	2.09	1.49	2.24	2.28		0.32	0.07	0.66	0.94	0.4	0.71	1.39
	2.83	3.86	1.41	3.26	1.63	1.27	1.09	0.22	1.75	1.54	0.75	1.57	1.53	0.69	1.04	1.46
	0.	0.	0.	0	0.18	1.97	2.34	1.81	0.	0.	0.17	0.	0.	0.3	0.52	1.13
	2.17 3.22 0.79	2.82	0.87	0.69	0.54	1.18	1.41	0.78	0.26	0.88	1.16	0.31	0.31	0.67
				DB 2 +EQV					DH-NEM2+EQV
						0.	1.28							0.	0.46
						4.71 3.67							1.99 2.64
						2.84	2.81							0.67	0.83
					1.2	0.28	0.						0.81	0.81	0.54
					1.07	2.19	4.1						2.52	1.44	2.18
					0.6	0.97	2.97						0.67	0.62	0.71
				2.46 1.82	0.	0.	0.52				0.7	0.79	0.4	0.	0.44
				0.68	0.65	2.08	4.54	4.41				1.97	1.06	2.3	2.18	2.57
				1.87	0.87	1.31	2.71 3.37				0.54	0.6	0.64	0.77	0.97
			0.87	2.08	0.74	0.7	0.55	0.24			0.59	0.8	0.64	0.35	0.62	0.48
			0.	0.	0.61	2.47	1.25	3.46			2.15	1.07	2.16	2.03	1.49	2.25
			0.57	1.41	0.35	0.62	0.69	1.67			0.52	0.61	0.57	0.5	0.63	0.69
		1.5	1.94	2.17	2.25	1.42	0.	0.		0.66	0.55	0.55	0.95	0.59	0.04	0.
		0.	0.	1.03	0.	0.25	3.75	3.43		1.87	1.95	1.9	1.08	2.3	2.11	2.55
		0.83	1.25	1.52	1.33	0.56	1.71	2.69		0.53	0.56	0.44	0.69	0.68	0.55 1.02
	2.28	2.2	1.14	2.33	1.94	1.23	0.82	0.	1.04 1.03	0.8	0.92	0.86	0.42	0.77	0.71
	0.	0.	0.	0.	0.	1.57	1.41	2.44	1.96	0.85	1.63	0.97	2.21	1.88	1.41	2.3
	2.15	1.67	0.72	1.71	1.5	0.72	0.54	1.63	0.82	0.72	0.57	0.67	0.68	0.46	0.64	0.66
				DB 2 +FDF						DH-NEM2+FDR

Respectively MAX, MIN and RMS of relative pin power error [%] within assembly. 31

Table 2 :

 2 Run-time comparison. All calculations use 40 parallel process.

		k-k ref	RMS N. global iterations Run-Time
	DT-CMFD	k ref =0.97449	ref	12	58 min
	DT-CMFD3x3	ref	ref	18	81 min
	DT-NEM4+FDR	ref	ref	13	68 min
	DH-NEM2+FDR+BB	ref	ref	11	60 min
	DH-NEM4+FDR+BB	ref	ref	13	75 min
	DH-NEM2+FDR	-21 pcm	0.65%	11	57 min
	DH-NEM2+EQV	-35 pcm	0.94%	10	50 min
	DH-NEM4+EQV	-35 pcm	1.47%	11	56 min
	DH-NEM2	-131 pcm	1.9%	10	50 min
	DB2+FDF	-20 pcm	1.7%	0	10 sec
	DB2+EQV	-38 pcm	1.9%	0	6 sec
	DB2	-104 pcm	2.2%	0	6 sec

(a) One-medium reflector model. (b) Four-motif reflector model

CONCLUSIONS

In this work we have investigated a method of Dynamic Homogenization (DH) that we propose as an alternative technique for full reactor calculations and analysis. Objective of the method is to solve the full core problem with a one-step approach via an iterative process, between fine transport assembly calculations and two-group diffusion core calculation, preserving assembly macro-exchanges. In contrast with an iterative scheme used in a nonlinear acceleration for the transport problem, in this approach the transport operator does not serve as core solver, but rather as generator of homogenization parameters, which entails that the quality of the solution depends on the coarse core operator. However, the advantage of knowing the actual conditions of the core and the environment of each assembly, allows the DH method to seek for a better reference homogenization problem, that provides a weighting flux for cross sections homogenization which is closer to the real situation than the one given by the critical infinite lattice conditions adopted in the classical two-step calculation scheme.

Besides avoiding the infinite lattice approximation and the critical leakage model, the DH approach has other advantages such as follow: (i) no need for a power reconstruction technique, since the transport solution is locally available in every assembly; (ii) no need for a multiparametrized cross-section library and, therefore, no need for interpolation of the homogenized data, since the assembly transport calculation can be performed at the actual core conditions. Moreover, (iii) the depletion calculation can be performed using transport fluxes at the pin-level and imposing the assembly power density from the actual power distribution. Another advantage of the method is (iv) the possibility to homogenize the reflector using a heterogeneous geometry and taking into account the real environment of the core. Furthermore, (v) the method avoids expensive 3D transport calculations. This makes definitely DH a good candidate as one-step calculation scheme for core analysis and design. On the other hand, because each assembly in the core has to be computed, the DH approach is feasible only in a parallel framework.

In this work we have applied the method to a core calculation for a two-dimensional configuration. Before incurring to the 3D, with the present analysis we firstly wanted to verify whether the main source of the typical two-step errors is due to the inadequacy of the adjusted low-order operator to describe the physics of the full core problem, or that of the reference homogenization problem to provide a good approximation for the core flux, or even both. A nuclear reactor is radially more heterogeneous than axially, therefore a diffusion operator can more easily fail. Moreover, for our 2D core calculations 40 process were sufficient, so they turned in a desktop machine, but a real 3D core calculation largely requires more computational resources.

Despite the influence on the core solution of the environment information accounted at the homogenization stage, we studied the applicability and the impact in DH of the two main techniques used to preserve the transport reaction rates: the equivalence theory and the flux discontinuity ratios. The evidence that comes out from our calculations is that these techniques result to be more effective in DH than in two-step, where the homogenization parameters are affected by the fundamental mode approximation. We also showed how any coarse operator constructed with the flux discontinuity ratios technique can be used as nonlinear acceleration for the full core transport problem.

In 2D configurations the DH method exhibits a convergence rate which is comparable with re-

Appendix A Nodal Diffusion (NEM)

The Nodal Expansion Method (NEM) that we used in our calculations has been implemented by Sanchez et al. in [START_REF] Sanchez | Diffusion Piecewise Homogenization via Flux Discontinuity Ratios[END_REF]. NEM is a class of methods that takes advantage of the Cartesian geometry by applying a transverse operator of the form 1 ∆y∆z ∆y ∆z dydz• to the 3D diffusion equation (called transverse integration procedure) in order to obtain a system of three one-dimensional equations coupled by the transverse leakage source term. Each one-dimensional equation have the following form defined for a general axis x:

Note that the coefficients of Q x also depend on the last transport flux but the source is computed using the coarse flux. The angular flux solved by DDM is accelerated with: