Paolo Cattaneo
email: paolo.cattaneo@cea.fr

Roland Lenain

Elsa Merle

Cyril Patricot

Didier Schneider

Numerical Optimization of a Multiphysics Calculation Scheme Based on Partial Convergence

Keywords: Multiphysics, Fixed-Point, Anderson, Partial-Convergence, APOLLO3 ®, FLICA4, MENDEL, PWR

This paper deals with the numerical optimization of a multiphysics calculation scheme. The purpose of this tool is the fine-scale neutronic and thermalhydraulic modelling of Pressurized Water Reactors (PWRs), under steady-state nominal conditions and fission products equilibrium concentrations. The neutronic model follows a two-steps approach with pin-cell homogenization. The considered coupling scheme includes the neutronic core model, the subchannel thermal-hydraulic and heat conduction one and the isotopic depletion one. From the numerical standpoint, this problem is a large system of non-linear equations.

For its resolution, two standard numerical methods are considered, the damped fixed-point method and the Anderson acceleration. In particular, this work focuses on the application of the partial-convergences within these methods. The partial-convergence technique deals with the research of a progressive internal convergence of the considered neutronic and thermal-hydraulic solvers. Within this approach, the degree of convergence is controlled either by limiting the number of single-solver iterations or by adjusting the required precisions on the respective key variables. To realise a large number of simulations in an affordable time, the chosen case study is a mini-core (5x5 PWR fuel assemblies plus

Background and motivation

This work concerns the numerical optimization of a steady-state coupling scheme for the simulations of Pressurized Water Reactors (PWRs) under normal operating conditions. In almost every scenario, except at zero power, this type of modelling is regarded as a coupled multiphysics problem. In fact, the heat generated by the neutron interactions with matter affects the density and temperature of the materials, producing Doppler and moderator effects that in turn condition neutron reaction probabilities which acts to modify the distribution of the power. Recently, many research groups worldwide have demonstrated an increasing interest for the subject [START_REF] Kochunas | VERA core simulator methodology for PWR cycle depletion[END_REF][START_REF] Ryu | Solution of the BEAVRS benchmark using the nTRACER direct whole core calculation code[END_REF][START_REF]Advanced multi-physics simulation for reactor safety in the framework of the NURESAFE project[END_REF]. With respect to the neutronic models, [START_REF] Kochunas | VERA core simulator methodology for PWR cycle depletion[END_REF][START_REF] Ryu | Solution of the BEAVRS benchmark using the nTRACER direct whole core calculation code[END_REF] deal with high fidelity simulations based on direct core calculations and massive parallelization, while the present work and [START_REF]Advanced multi-physics simulation for reactor safety in the framework of the NURESAFE project[END_REF] still rely on the homogenization process, hence require less computing power. In particular, in [START_REF]Advanced multi-physics simulation for reactor safety in the framework of the NURESAFE project[END_REF], the pin power distribution is obtained combining a coarse mesh homogenization with a form function (referring to [START_REF] Gomez-Torres | Implementation of a fast running full core pin power reconstruction method in DYN3D[END_REF] for a description of this method).

In the current paper, pin-cell homogenization is adopted, therefore, in terms of computing time and accuracy this choice could be considered in between the two other approaches. In this work, in the cited sources, subchannel thermalhydraulic modelling and heat conduction in every fuel rod slice are chosen.

From the point of view of numerics, the coupled problem constitutes a large nonlinear system of equations. It may be described by a residual function F(x) = 0 or equivalently by a fixed-point function of the type G(x) = x, with F, G : R n → R n , n ∈ N. In most of the considered applications, evaluating any of these functions is equivalently computationally expensive and the calculation of the derivatives is prohibitive. For this reason, the cost of forming the Jacobian matrix hinders the implementation of the exact Newton algorithm. There is therefore a rich literature in terms of numerical methods that do not require the computation of the derivatives, for example [START_REF] Brent | Algorithms for minimization without derivatives[END_REF]. The common approach to solve the coupled problem is to combine pre-existing specialized solvers that treat only a subset of the physics. In particular, if it is chosen to avoid intrusive modifications of the solvers, the technique takes the name of black-box coupling [START_REF] Keyes | Multiphysics simulations: Challenges and opportunities[END_REF]. The advantages in terms of simplicity and modularity are clear. Whilst, the main inherent disadvantage is that, in this way, it is possible to access only to a limited set of variables. Hence, some advanced algorithms conceived specifically for a given multiphysics coupling scheme may not be viable because of the need for unavailable internal quantities (an example is given in [START_REF] Zhang | The comparison between nonlinear and linear preconditioning JFNK method for transient neutronics/thermal-hydraulics coupling problem[END_REF]). This could constitute a limit in terms of stability and robustness for very strongly coupled problems, but this is generally not the case for PWRs under steady-state nominal conditions, as confirmed by the successful applications reported in [START_REF] Hamilton | An assessment of coupling algorithms for nuclear reactor core physics simulations[END_REF][START_REF] Toth | Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system[END_REF].

In this framework, the damped fixed-point iterations have been widely used in reactor physics simulations, for instance in [START_REF] Schmidt | Foundational development of an advanced nuclear reactor integrated safety code[END_REF][START_REF] Zhang | The improvement of coupling method in TINTE by fully implicit scheme[END_REF]. The relaxation factor used to accelerate and stabilize the fixed-point iterations is very simple and effective.

However, optimal values of this parameter are case dependent and difficult to determine, most often they are derived empirically. The strong dependency on the relaxation factor is one of the major problems of the fixed-point method.

Given the high cost needed for the evaluation of the multiphysics function, it is important to try to exploit to the greatest extent the information produced during iterations. For this reason, many works focus on more advanced methods, the most common are Anderson acceleration [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF] and Jacobian-Free Newton-Krylov (JFNK) [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF]; some references to their application in reactor modelling are [START_REF] Toth | Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system[END_REF][START_REF] Hamilton | An assessment of coupling algorithms for nuclear reactor core physics simulations[END_REF][START_REF] Zhang | The comparison between nonlinear and linear preconditioning JFNK method for transient neutronics/thermal-hydraulics coupling problem[END_REF]. Nevertheless, as observed in [START_REF] Schmidt | Foundational development of an advanced nuclear reactor integrated safety code[END_REF], complex schemes like JFNK do not always perform better than the fixed-point method, in part due to the fact that their efficiency is highly dependent on the preconditioning. An alternative, rather widespread in the industry, but not so commonly found in literature is the generalised fixed-point with partial-convergences. This method is based on the research of a progressive internal convergence of the solvers by limiting the maximum number of single-physics iterations performed at each multiphysics evaluation. Most frequently the partial-convergences strategy is applied only to the neutronic solver as the algorithms used for the steady-state eigenvalue problem are often particularly suitable for this type of truncation. An example of generalised fixed-point with partial-convergences is given by [START_REF] Clarno | High fidelity modeling of pellet-clad interaction using the CASL virtual 40 environment for reactor applications[END_REF]. In particular, the focus of this work is on the application of the partial-convergences to the damped fixed-point and Anderson algorithms. A mathematical analysis of the convergence process of the fixed-point and Anderson iterations in case of degradation of the internal convergence criteria of the specialised solvers is provided by [START_REF] Toth | Local improvement results for Anderson acceleration with inaccurate function evaluations[END_REF].

The first three sections are respectively dedicated to the problem formalization, the modelling details and the introduction to the damped fixed-point algorithm. To apply the numerical methods on a practical test, a case study is defined in the fifth section. In the sixth section, the generalised fixed-point algorithm is defined and its performance is analysed. In particular, a wide range of limitations on both the neutronic and the thermal-hydraulic iterations is explored. In the seventh section, the implementation of the Anderson algorithm is described and its performance is compared to that of the generalised fixed-point with partial-convergences. Finally, the integration of the partial-convergences to the Anderson algorithm is studied.

Problem formalization

The objective of this particular multiphysics calculation scheme is to find the steady-state solution of the coupled system of neutronics, thermal-hydraulics and fission products evolution. The problem is formulated as it appears in the system of Eq. [START_REF] Kochunas | VERA core simulator methodology for PWR cycle depletion[END_REF].

           (ρ w , T f ) = T H(q w , q f ) (φ, q w , q f ) = N (ρ w , T f , C F P ) C F P = F P (φ) (1)
Where TH is the thermal-hydraulic operator, which determines the water density (ρ w ) and the fuel temperature (T f ) corresponding to a given heat source.

N is the neutronic function, which solves the k -eigenvalue criticality problem and computes the multi-groups scalar flux (φ) and the correlated power generation in the fuel (q f ) and in the water (q w ). The neutron scalar flux and the power distributions are normalised to keep the total power as defined by the case study. The shape of the power distribution depends on the other physics via the macroscopic cross-sections, which are highly influenced by the density and temperature fields and by the fission products concentrations. Indeed, these quantities affect both the particle number densities and the neutron energy spectrum. The role of FP is to find the equilibrium concentrations of the fission products (C F P ) corresponding to (φ). The research of the critical boron concentration and of the critical control rod insertion might be included in the problem, but their application is not considered in this work.

Modelling details

This work is carried out using a combination of specialised solvers and exploiting the tools of the SALOME platform [START_REF] Ribes | SALOME platform component model for numerical simulation[END_REF], which facilitate the manipulation, the exchange and the storage of data.

The neutronic solvers used in this study belong to APOLLO3 ® [START_REF] Schneider | APOLLO3: CEA/DEN deterministic multi-purpose code for reactor physics analysis[END_REF]. Given the magnitude of the problem, the neutronic modelling follows the two-steps approach. During the first step, lattice calculations, the dimension of the problem is reduced by computing very precisely several 2D fuel assemblies under simplified boundary conditions and for a set of possible core parameters. In this way, a set of pin-cell homogenized cross-sections is produced in order to define an equivalent problem on coarser energy and spatial meshes. The chosen parameters are burnup, fuel Doppler temperature and moderator density. Via this parametrization it is possible to synthesize, under simplifying hypotheses, how the neutronic cross-sections are influenced by the other physics. This part of the scheme is carried out independently from the coupling scheme and it could be seen as a preparation phase. For this stage, the Method Of Characteristics solver TDT [START_REF] Santandrea | Accelerated polynomial axial expansions for full 3D neutron transport moc in the APOLLO3 ® code system as applied to the ASTRID fast breeder reactor[END_REF] is used. In order to account for the fact that the cross-sections are to be utilised by a solver which is different from the lattice one, the SuPer-Homogenization technique [START_REF] Kavenoky | The SPH homogenization method[END_REF] is adopted. The second step of the neutronic model, the core calculation, deals with the calculation of the entire domain.

After a comparative performance analysis [START_REF] Cattaneo | Development of a Multiphysics Best-Estimate Approach for LWR Reference Calculation[END_REF], it has been chosen to use the Simplified Spherical Harmonics [START_REF] Baudron | MINOS: a simplified PN solver for core calculation[END_REF] of order three with eight energy groups. The reflector's cross-sections are homogenized with the simplifed approach of 1D traverses [START_REF] Sandrin | An analysis of reflector homogenization techniques for full core diffusion calculations[END_REF].

The research of equilibrium concentration of the fission products is performed by the MENDEL-solver library [START_REF] Lahaye | First verification and validation steps of MENDEL release 1.0 cycle code system[END_REF]. The evolution of the eight isotopes that impact the most on the reactivity is considered: xenon-135 and its precursor iodine-135, samarium-149, its precursors neodymium-147, promethium-147, 148, 148m and 149. The equilibrium concentrations is targeted for both the decay chains in order to find an actual steady-state solution of the coupled system.

The thermal-hydraulic model is FLICA4 [START_REF] Toumi | FLICA-4: a threedimensional two-phase flow computer code with advanced numerical methods for nuclear applications[END_REF], it is used to perform subchannel calculations with the 3D 4-equations modelling in porous media and 1D heat conduction for all the fuel rods. In Fig. 1, the radial discretizations for the neutronics, thermal-hydraulics and heat conduction are depicted.

Introduction to the damped fixed-point algorithm

Mathematical background

In the context of multiphysics calculation schemes, the fixed-point algorithm is one of the simplest and most widespread. In fact, it is sufficient to formalize the problem in the form G(x) = x, under which it often naturally falls, to create a sequence of the type of Eq. ( 2). Nevertheless, multiple ways exist to define G, corresponding to different numerical methods. In this context, x is the vector of all the variables: multi-groups neutron scalar flux, heat sources, water density, fuel Doppler temperature and fission products concentrations. The solution of the problem is x * , respecting G(x * ) = x * .

x k+1 = G(x k ) (2)
A sufficient condition for local convergence of the fixed-point method is given by the following statement. This notation is adopted for the Jacobian matrix of G, J G = ∂G ∂x and for its spectral radius (i.e. the maximum in module of the eigenvalues of the matrix), ρ(J G (x)). If G has a fixed-point, G ∈ C 1 in its proximity and ρ(J G (x * )) < 1 then there exists a neighbourhood of the solution in which any initial guess makes the sequence converge [START_REF] Quarteroni | Numerical mathematics[END_REF], page 297. Moreover, all the values of the sequence stay within this interval. In respect of the speed of convergence, it could be demonstrated (e.g. [START_REF] Quarteroni | Numerical mathematics[END_REF], page 262 for the scalar case) that, except in case of null spectral radius, the fixed point converges linearly 1and its rate of convergence is equal to the spectral radius, following Eq. (3). Such a variable is almost never computed in practical applications, because to derive it, the solution of the problem should be known and even to compute the Jacobian in a given point is generally prohibitive. However, this analysis is important to better understand the fixed-point method and the relaxation (or damping) technique.

lim k→∞ x k+1 -x * x k -x * = ρ(J G (x * )) (3)
As stated before there is not a unique way of writing G. Starting from the problem in the residual form (F (x) = 0), any formulation of the type of Eq. ( 4) would be valid under the condition that H is a homogeneous and continuous operator.

G H (x) = x -H(F (x)) (4)
Homogeneous meaning that H respects H(t * x, t * y) = t n * H(x, y) for a constant n. For this reason, a simple choice of H would be a multiplicative constant, which is referred to as the relaxation factor (α). In this way, a new sequence G α is built, whose relation with G (G ≡ G I ) is specified in Eq. [START_REF] Brent | Algorithms for minimization without derivatives[END_REF].

x k+1 = G α (x k ) = x k -α * F (x k ) = α * G I (x k ) + (1 -α) * x k (5)
It is easy to demonstrate that the eigenvalues of G α are linearly related to those of G, as written in Eq. [START_REF] Keyes | Multiphysics simulations: Challenges and opportunities[END_REF]. Hence, such a simple approach can transform a non convergent sequence into a convergent one defined by Eq. [START_REF] Zhang | The comparison between nonlinear and linear preconditioning JFNK method for transient neutronics/thermal-hydraulics coupling problem[END_REF].

λ i (α) = α * (λ i -1) + 1 (6) max(λ i ) < 1 ⇒ ∃ α > 0 : |λ i (α)| < 1, ∀i (7)
Moreover, it can also accelerate its convergence by reducing the spectral radius, as suggested by Eq. ( 3). The relaxation does not impact all the eigenvalues in the same way. Considering non-relaxed eigenvalues smaller than one, in order to nullify a positive eigenvalue, α should be greater than the unity, while for any negative eigenvalue α should be comprised between zero and one.

The more the maximum and the minimum eigenvalues are distant, which can be expressed by the dominance ratio, the less room for optimization there is.

In particular, to minimize the spectral radius, it appears that the optimal α is where |λ M (α)| = |λ m (α)|, respectively being the minimum and the maximum, which brings to Eq. [START_REF] Hamilton | An assessment of coupling algorithms for nuclear reactor core physics simulations[END_REF].

α opt = 1 1 -λ M +λm 2 (8)
The entire derivation of these equations and a wider analysis of numerical methods for multiphysics coupling can be found in [START_REF] Patricot | Multi-physics couplings: methodology impact evaluation for neutron transport /heat transfer /mechanics coupling simulations[END_REF]. In Fig. 2, it is possible to visualize an example of this optimization process for arbitrary eigenvalues.

Algorithm definition

The damped fixed-point algorithm is available in Algorithm 1. It is built directly from the system of Eq. (1), applying the relaxation only on the neutronic operator. It should be noticed that the damping could also be performed on other variables, for instance in [START_REF] Hamilton | An assessment of coupling algorithms for nuclear reactor core physics simulations[END_REF] and [START_REF] Toth | Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system[END_REF], the damping parameter is applied on the thermal-hydraulic solution. The convergence is separately tested on the L2norm of the absolute residual for each variable vector x i (where x is the linear concatenation of all the x i ). In order to obtain a convergence criterion independent from the relaxation factor, if the relaxation is applied, the convergence criteria are multiplied by α.

To link the current damped fixed-point iteration function to the mathematical background provided in the previous section, it is possible to think of the 0.00 to α = 1 (dash-dotted blue line, Undamped). λ tending to zero means optimal rate of convergence (horizontal dashed blue line, Opt abs). In order to minimize the spectral radius the optimum relaxation factor is where the maximum and minimum eigenvalues are equal in modules. In the example it is 0.77 (green dash-dotted line, Opt rel).

global functional as a composition of the operators using the neutronic variables as the only unknown. In fact, this algorithm could be seen as a nonlinear elimination of the "non-neutronic" unknowns (the moderator density, fuel effective temperature and fission products equilibrium concentrations) within a damped neutronic solver. The exact iteration function associated to the fixed-point defined by Algorithm 1 is the one given in Eq. ( 9). Where

x k := [φ k , q k w , q k f ]. F P (x 0
) is the concentration of fission products as stored in the cross-sections.

G α (x k ) := α * N ([T H(x k )] ρ w , [T H(x k )] T f , F P (x k-1 )) + (1 -α) * x k (9)
From this point of view, the convergence tests on all the variables not contained in x k become additional checks on internal variables of interest.

In order to evaluate the performance of the algorithm, a common estimator Algorithm 1 Damped fixed-point.

Flat power distribution

while ∆x k i > ∆ i * α, for any i do (ρ k+1 w , T k+1 f ) = T H(q k w , q k f ) (φ k+1 , q k+1 w , q k+1 f ) = α * N (ρ k+1 w , T k+1 f , C k F P ) + (1 -α) * (φ k , q k w , q k f ) C k+1 F P = F P (φ k+1 ) end while
of the type of Eq. ( 10) is defined.

t eqv := I N * t 1N + I T H * t 1T H + I F P * (t F P + t ex. ) (10)
This estimator defines an equivalent calculation time (t eqv ) based on the number of single-solver iterations multiplied by their respective average computing time, estimated once for all the applications. This choice is based on the fact that the iterations have an almost constant time cost and they carry most of the computing time. Moreover, it allows to compare the results obtained from calculations carried out on machines with different characteristics. I N is the number of neutronic power iterations, I T H the number of thermal-hydraulic time steps for the pseudo-transient associated to the steady-state resolution and I F P the number of depletion calculation. I F P is also equivalent to the number of multiphysics iterations. For this reason, it is multiplied by an average time for the full depletion calculation summed to an estimation of the fixed time cost for data exchanges and manipulations. In all the paper, this estimator is normalized to the best performing algorithm, as shown in Eq. ( 11) for this reason it is called relative equivalent calculation time (t eqv,r ).

t eqv,r := t eqv,i t eqv,min (11)

Case study specifications

The simulated domain is a mini-core made of 5x5 PWR assemblies, with radial and axial reflectors included. The radial reflector contains 95% of stainless steel and 5% of water. All the fuel assemblies are loaded with 4% enriched urania, their 2D burnup distribution is given in Fig. 3. The operating conditions are nominal: the average linear power is 160W/cm and the mass flux is 3900kg/m 2 /s. The fuel-cladding gap heat transfer coefficient is set constant in all the core to the value of 5000W/m 2 /K, corresponding to a value typical of the beginning of cycle and standard correlations are used for the conductivities and the water heat transfer coefficient. The core size entails a power peaking factor higher than the standard value for PWRs, which implies a higher maximum power and, therefore, a stronger local coupling. Such a property is desirable to test the numerical schemes under potentially more complicated conditions.

Moreover, the small size facilitates the systematic study, which can be carried with faster calculations. 6. Application of partial-convergences to the fixed-point method

Damped fixed-point application

The damped fixed-point algorithm defined in the previous subsection is tested for a typical range of relaxation factors. The plot of the performance versus the relaxation factor is available in Fig. 4. When the algorithm converges, the associated equivalent calculation time is reported. Otherwise, if the L2-norm of the residuals do not decrease enough over the iterations, "NC" is marked on top of the relaxation factor, whereas, if the L2-norms of the residuals increase, "DIV" is used. Without damping, divergent axial oscillations are found for the neutron scalar flux. Due to the low value of the total heat transfer coefficient for the extraction of the heat from the fuel, also the fuel temperature follows similar oscillations. With respect to the density, the variations are prevalently axial.

The axial position of the peak of heat flux significantly impacts the density axial profile and the convergence speed. The divergent behaviour makes the solution exit the range of validity of the parameters, hence the iterations stop. These results are in line with what found in literature, for instance refer to [START_REF] Toth | Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system[END_REF] (single fuel assembly case study). When the oscillations are observed, the eigenvalues are expected to be negative and, due to the divergent behaviour, some of them are supposed to be inferior to minus one. Hence, the optimal relaxation factor is expected to be between zero and one.

Considering Fig. 4, a range of relaxation factors that makes the scheme converge is found, it contains the values between 0.08 and 0.39. From a theoretical standpoint, values comprised in α ∈ ]0, 0.08[ should bring to convergence in a decreasingly effective way approaching zero. Obviously α = 0 is not considered, because it means no update of the solution. Nevertheless, small values of α might bring to convergence criteria too small compared to the internal precision of the solvers. Indeed, for α = 0.05, a plateau is observed just before reaching the convergence for some physics. So even if this case is labelled as "NC", the scheme is actually converging to a result, but with a precision lower than the imposed one. To solve this problem, the internal convergence criteria should be refined, but, in this work, this option is not considered. On the contrary, for α ≥ 0.4, an interval of non-convergence is observed, in the sense that the solution is trapped in periodic oscillations within a set of values. More insight on this behaviour is given in the next sub-section. For even larger relaxation factors, the module of the residuals increases with the iterations. Given the small range of acceptable relaxation factors, the damped fixed-point algorithm could be seen as not suitable for the current application. For such a strongly coupled problem, tighter algorithms should be used. Moreover, the performance is very dependent on the value of the relaxation factor, which is very inconvenient, because optimal values are generally not available or determined empirically.

Generalised fixed-point with partial-convergences 6.2.1. Introduction

The partial convergence technique consists in chasing a progressive convergence of some of the individual physics solvers. Within this approach, the internal convergence of the neutronic and thermal-hydraulic operators is not asked at every multiphysics iteration, but it becomes just a necessary criterion for the global convergence of the scheme. The neutronic operator solves the k -eigenvalue problem associated to the Boltzmann equation with the power iteration method. For this reason, the partial convergence can be imposed in two ways: by limiting the number of power iterations to be done at each multiphysics call or by asking the solver for a precision on the solution that increases with the multiphysics iterations. These two options are not equivalent at all, depending on the initialization of the solver, a widely varying number of iterations would be needed to obtain a given precision. The thermal-hydraulic operator searches the asymptotic solution of the pseudo-transient associated to the steady-state problem. For this reason, in a similar way, the partial convergence can be controlled by imposing a maximum number of time steps or by asking an increasing precision on the solution. In this work, a constant time step of 0.01 s is considered and it is tested to limit the number of power iterations and the number of time steps (both being also referred to as single-physics iterations). Since the fission products equilibrium research takes a much lower computing time and its slope of convergence may be rather stiff and not monotonous, full convergence is always asked for this calculation. Within this approach, a key role is played by the solver initialization at each call. The solver restarts from the results computed at the last call. In this way, in case of stable algorithm, the number of iterations needed to meet the convergence criteria reduces with the multiphysics calls and after a certain point, the internal convergence becomes reachable within a number of single-physics iterations lower than the limits.

The main idea behind this technique is to avoid the extra single-solver iterations needed to very precisely compute every physics while away from the solution. Instead, a progressive convergence of all the variables is targeted. Another justification to the use of this technique could be found in its similarities with the under-relaxation (α < 1). When the convergence of the internal solvers is monotonous, the solution with partial convergence is "contained in between" the result from the previous iteration and the fully converged new guess. For this reason, when the under-relaxation is desirable, the partial convergence allows to achieve similar savings in terms of multiphysics iterations, but also to reduce the number of single-solver iterations, by avoiding the unnecessary iterations to full convergence when far from the solution. On the other hand, it should be noticed that the exact operator N is substituted by an operator N k nint. that depends on the maximum number of power iterations "n int. " and on the initialization at the iteration "k". If the internal convergence is met, the outputs of N k nint. and N must coincide within the fixed precision and regardless of the parameters. For this reason, as stated before, the internal convergence is necessary for the global convergence. The same observations are valid for the thermal-hydraulics ("m int. " is used for the limit on the time steps) and for the global operators which is described in Eq. [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF]. In general, given the nature of the fixed-point algorithm and of the solvers, this loss of consistency should not cause any problem. Algorithm 2 represents the generalised fixed-point with partial-convergences, using the previous notation for the partial operators.

G k α (x k ) := α * N k ([T H k (x k )] ρ w , [T H k (x k )] T f , F P (x k-1 )) + (1 -α) * x k (12)
Algorithm 2 Generalised fixed-point with partial-convergences.

Flat power distribution

while ∆x k i > ∆ i * α, for any i do (ρ k+1 w , T k+1 f ) = T H k mint. (q k w , q k f ) (φ k+1 , q k+1 w , q k+1 f ) = α * N k nint. (ρ k+1 w , T k+1 f , C k F P ) + (1 -α) * (φ k , q k w , q k f ) C k+1
F P = F P (φ k+1 ) end while

Systematic study of the parameters choice

This algorithm is applied to the case study previously defined, for different tuples of maximum single-solver iterations ("N N -N T H ") and for a range of relaxation factors. The objective is to show the interplay between these parameters and the potential benefits this technique could bring2 .

The first test concerns a symmetrical reduction of both the neutronic and the thermal-hydraulic iterations limits with respect to the total number of itera-tions required to solve the neutronics for isothermal conditions and the thermalhydraulics for flat power distribution.

The number of multiphysics iterations required to solve the full problem, normalised to the lowest obtained value, is plotted in Fig. 5 for the considered tuples of limitations. From these results, it clearly appears that lowering the limits on the singlephysics iterations leads to a significant reduction of sensibility on the damping factor. In the extreme case, in which only one neutronic iteration and a thermal-hydraulic one are performed during each solver call ("1-1"), the algorithm converges over the entire range of relaxation factor and the number of required multiphysics iterations is almost constant. Moreover, when limiting the number of single-physics iterations, except for "1-1", the number of multiphysics iterations needed to fully solve the problem does not increase much. For this reason, savings on the total computing time are expected. Therefore, the same data are analysed in terms of relative equivalent calculation time, the plot is available in Fig. 6. In this figure, "inf-inf" corresponds to no limits on the iterations, so to the same curve of Fig. 4, but with the current normalisation. Thanks to the significantly lower time cost of each multiphysics iteration, for all the explored settings the iterations actually converge faster than the full convergence fixed-point. The equivalent calculation time can be up to 15 times smaller than the best result obtained with the standard algorithm. Also in terms of performance, it clearly appears that the dependency on the relaxation factor is drastically reduced. Even if it might seem that the problem is only transported from the choice of a relaxation factor to that of the tuple of maximum singlesolver iterations, the dependency on such a parameter is much smaller for any choice within few tens of iterations.

The comparison between "1-1","10-10" and "20-40" shows that alternating the resolution of the two physics at each iteration is a suitable choice. However, the considered time cost for the data manipulation and exchange during a multiphysics iteration favours the choice of few tens of single-physics iterations (e.g.

"10-10", "20-40"). Moreover, imposing the "1-1" setting makes impossible to test the internal convergence on more than one iterate and the dedicated inter-nal acceleration techniques, which are usually very effective, are not applicable.

For instance, in the resolution of neutronic power iterations, Chebyshev acceleration might be desirable, but the internal solver would have access to only one iterate, hence no acceleration technique would be available in this case.

The second test deals with the asymmetrical reduction of the limits on the number of single-physics iterations. The results in terms of normalised multiphysics iterations appear in Fig. 7. The results show that limiting the number of neutronic iterations has a more stabilizing effect than imposing a lower maximum of thermal-hydraulic iterations. In fact, the abrupt loss of performance and stability, observed for relaxation factors slightly larger than the optimum value, only appears when the limit on neutronic iterations is set to 100 or more. Furthermore, in most of the cases, limiting the neutronic iterations brings to larger α opt . This could be explained with the reduction in module of the most negative eigenvalues of the iteration function. The extreme cases "1-inf" and "inf-1" fully confirm this trend. In particular, "inf-1" diverges for almost the entire range of standard relaxation factors and, on the contrary, "1-inf" converges for the entire range of damping factors. To compare the performance corresponding to these settings, they are plotted in terms of relative equivalent calculation time in Fig. 8. The asymmetrical reduction of the limits on the single-physics iterations seems to be slightly less effective in terms of performance improvement. The extreme settings ("1-inf","10-inf","inf-10" and "inf-1") confirm that they are not particularly attractive: "1-inf" and "10-inf", which appear as the best compromises out of the four, are more robust than the full convergence fixed-point, but they do not offer large time savings.

The results reported so far are obtained assuming that the time for the data exchange and manipulation is in the order of 2 s. In order to test their sensitivity to the manipulation time, the curves are plotted again in Fig. 9 for a data exchange and manipulation time of 200 s. The shape of the curves does not drastically change, which makes the study slightly more general, in terms of small sensitivity of the analysis to the exact estimation of this time cost.

However, in this case, algorithms with very small limits on the single-physics iterations are strongly disadvantaged. In particular, the main difference is that the "1-1" setting is no longer competitive in terms of equivalent calculation time.

Ultimately, this result underlines the importance of optimising the data manipulations and exchanges when using the partial-convergences, as this technique may potentially lead to a larger number of multiphysics iterations. An analysis, similar to that reported in this section, has been published by [START_REF] Clarno | High fidelity modeling of pellet-clad interaction using the CASL virtual 40 environment for reactor applications[END_REF]. A single rod case study is adopted and only the limitation of the neutronic iterations is analysed for a given damping factor. Therefore, the relevant settings that should be compared are "1-inf", "10-inf" and "inf-inf". Considering that a different combination of solvers is employed, the results seem to be in rather good agreement.

Analysing the number of neutronic and thermal-hydraulic iterations separately, it is possible to get more insight on the efficiency of limiting the maximum number of iterations. It should be noticed, that in this kind of analysis, the time for the data manipulation and exchange does not have any influence. In Fig. 10a and 10c, the ratio of the total number of neutronic iterations to a reference number is reported. This reference corresponds to the amount of iterations needed to fully solve the neutronics considering the fission products concentrations as stored in the cross-sections and a given temperature and density profiles, which correspond to a flat power distribution. Similar plots are produced for the thermal-hydraulic iterations, see Fig. 10b and 10d, using as reference value the number of time steps required to fully solve the thermal-hydraulics for the flat power distribution. It should be noticed that the internal convergence of the thermal-hydraulics is tested on five consecutive iterations, in every case but for "1-1". For this setting, the convergence is only tested on one iteration, therefore, this method has an unjustly looser convergence criterion, which partially explains its outstanding performance.

From these results, it appears that the total number of neutronic and thermalhydraulic iterations to solve the coupled system with partial-convergences can be comparable or even lower than corresponding reference number ("20-40").

In case of asymmetrical reductions, it is found that, when N N is reduced, it is possible to make savings on the neutronic iterations at the cost of increasing the thermal-hydraulic ones and vice-versa when controlling N T H .

The fact that solving the whole coupled system could cost less single-physics iterations than solving the corresponding decoupled problem for the initial conditions is a rather unexpected result. It should be noticed that the neutronic problem defined by the initial conditions is different from that corresponding to the converged ones and, at least for this case study, the number of required iterations in the second case is significantly lower. In particular, the equilibrium fission products concentrations and the thermal-hydraulics tend to flatten the power distribution, making the neutronics easier to be solved. Similarly, the thermal-hydraulics often requires less iterations for lower power peaking factors. Hence, for similar cases and initialisations, the partial-convergences methods appears as very competitive.

In order to verify that all the schemes converge to the same solution, the discrepancies that occur on the main variables are measured. The RMS and the maximum of the absolute relative discrepancy are computed for all the converging settings using "inf-inf" with a relaxation factor of 0.3 as reference. This process is repeated for the integrated power in the fuel, the moderator density, the fuel effective temperature and the xenon-135 equilibrium concentration. The histograms relative to frequency of each discrepancy are reported for these four variables in Fig. 11 3 .

The results show that all the convergent settings lead to the same solution for all the interest variables within the 8 pcm of maximum absolute relative The examined variables are the integrated power in the fuel, the moderator density, the fuel effective temperature and the xenon-135 equilibrium concentration. discrepancy on the pin-cells, which are directly linked to the convergence criteria.

Referring again to [START_REF] Clarno | High fidelity modeling of pellet-clad interaction using the CASL virtual 40 environment for reactor applications[END_REF], it appears that limiting the neutronic power iterations (for a given damping factor) has potential to lead to false convergence on all the interest variables. In particular, it is observed that the discrepancies against the full-convergences results consistently increase when lowering the limit on the neutronic iterations. In the current work this is not found. Even if the fullconvergences fixed-point, with a relaxation factor similar to the reference (0.275 or 0.325), calculates the solution with the lowest RMS and maximum of the absolute relative discrepancy on some of the considered variables, their predictions are not always the closest to the reference (e.g. (α = 0.1, N N = 20, N T H = 40) computes the solution with the lowest RMS on the moderator density). Moreover, the methods with the smallest limitations on single-physics iterations ("1-1" and "10-10") do not even appear in the group of settings leading to the first percentile of largest discrepancies. There is not one setting producing results with the largest discrepancies on all the interest variables. These results could potentially derive from different internal controls on the internal variables of the specialised solvers.

Fixed-point bifurcations

The previously mentioned periodic oscillations of the solution observed in the non-convergence interval (acknowledged for every adopted fixed-point algorithm) correspond to the fixed-point bifurcations, which have been already studied by [START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF]. This behaviour is easier to visualize for a scalar quantity like the effective multiplication factor. An example of oscillations of period four is given in Fig. 12 4 .

It should be noticed that, like in this example, the four values periodically repeating may not contain at all the converged multiplication factor. Figure 12: Different convergence behaviours of the multiplication factor for two values of the relaxation parameter, here the "100-40" setting is considered. On the left-hand side, monotonous convergence is observed, while for slightly larger relaxation-factor, the previous mentioned non-converging oscillations of period four are reported.

Application of partial-convergences to Anderson acceleration

Introduction to the standard Anderson algorithm

The Anderson algorithm [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF] is a low-degree generalized secant method, which could also be seen as an acceleration method for the fixed-point method.

Indeed, it aims to exploit to a greater extent the results obtained during the previous iterations in order to create a faster converging sequence. More in details, this acceleration is based on the linear combination of the latest M (arbitrary parameter) iterations, using weights obtained from the minimization of a quadratic residual. This process actually corresponds to the formation of a rank-M approximation of the inverse Jacobian of F(x) as demonstrated by [START_REF] Eyert | A comparative study on methods for convergence acceleration of iterative vector sequences[END_REF].

There are multiple good reasons to prefer, among the quasi-Newton methods, a low order one. Since the evaluation of F(x) is very expensive and its solution has a precision limited by the internal solvers iterative processes, all the finite differences schemes, like JFNK, appear less competitive because of their lower tolerance to the noise. A more detailed analysis is available in [START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF]. This flexibility to a certain level of approximation in the function evaluation is a key feature for the applicability of partial-convergences. Furthermore, in the analysis presented by [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF] and [START_REF] Eyert | A comparative study on methods for convergence acceleration of iterative vector sequences[END_REF], it is argued that, approaching convergence, linear dependencies arise among the iterates, implying an increasing ill-conditioning of the matrix to be solved to obtain the M weights. Low-order approximations suffer less of this ill-conditioning.

For this work, the Anderson algorithm is imported from the Scipy-Optimize library [START_REF] Jones | Scipy: Open source scientific tools for python[END_REF]. This version of the algorithm offers more than the "simple" Anderson iterative scheme. Two main extensions are available, the possibility to apply a regularization parameter and the line-search. The regularization parameter reduces the ill-conditioning problem in change of a minor deceleration of the iteration process. Hence, it could be used to select higher M -parameter with a reduced loss of accuracy in the weights determination. The second extension allows to dynamically determine the weight used to combine the results at the previous and the current iterations in order to meet a criterion that links the descent slope and the equivalent step's size to be chosen. The criterion is given either by Wolfe's [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF][START_REF] Wolfe | Convergence conditions for ascent methods. ii: Some corrections[END_REF] or by Armijo's [START_REF] Armijo | Minimization of functions having lipschitz continuous first partial derivatives[END_REF] rules. Both the techniques may require several function evaluations before the condition is met. Due to the associated computation costs, this feature is not of particular interest in the current work, but it is still applied in some applications. Even if the relaxation technique is not totally overlapping with the line-search, in every application presented here, Anderson algorithm is used without relaxation. In fact, even Anderson suggested to not apply relaxation unless empirical experience is available [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF]. In order to use this library, it is only necessary to switch to the residual formalism. In addition to expressing the problem in terms of the residual, Eqs. ( 13)-( 15) introduce a change in notation, namely the use of normalized quantities and operators, both identified with a circumflex accent.

F (x k ) =      T H(q k w , qk f ) -(ρ k w , T k f ) N (ρ k+1 w , T k+1 f , Ĉk F P ) -( φk , qk w , qk f ) F P ( φk+1 ) - Ĉk F P      (13)
x := [q w * 1/c q w , ..., ρ w * 1/c ρ w , ...] = [q w , ..., ρw , ...]

T H( qw , qf ) := ([T H(q w , q

f )] ρ * 1/c ρ w , [T H(q w , q f )] T * 1/c T f ) (15)
This notation change derives from the desire to have an unknown vector containing quantities of the same order of magnitude. The scaling factors are the individual convergence criteria (c i ) for each variable vector. Following the definitions given in Eqs. ( 16) and ( 17), Algorithm 3 is built.

xk := xk + M k j=1 θ k j * (x k-j -xk ) (16)
F k := F (x k ) + M k j=1 θ k j * ( F (x k-j ) -F (x k )) (17)
Algorithm 3 Simple Anderson acceleration: no relaxation, no line-search and no regularization parameter.

x0 = Ĝ(q f lat ) k = 1 while ∆x k i > 1, for any i do M k = min(M, k) minimize R k for θ k j : M k j=1 ( F (x) k -F (x) k-l ) T • ( F (x) k -F (x) k-j )θ k j = ( F (x) k -F (x) k-l ) T • ( F (x) k ) ∀ l=1,...,M k xk+1 = xk + F k k = k + 1 end while
The convergence criterion is still tested on the L2-norm of the absolute difference of the variable for two consecutive iterations. The weights θ k j are determined by the minimization of the residual as defined in Eq. [START_REF] Santandrea | Accelerated polynomial axial expansions for full 3D neutron transport moc in the APOLLO3 ® code system as applied to the ASTRID fast breeder reactor[END_REF].

R k = 1 2 ( F k ) T • F k (18)

Comparison against the generalised fixed-point with partial-convergences

The benefits of the Anderson acceleration are expected both in terms of stability and convergence rate. In Fig. 13-a, the convergence of the power integrated in the water is reported for the fixed-point and for the Anderson method for a range of M -parameters. The standard fixed-point and the extended Anderson (regularization parameter and Armijo's line-search) algorithms are considered, both with no damping.

As shown in the previous section, the simple fixed-point scheme diverges after few iterations. Anderson with M=0, which is equivalent to the fixed-point method plus line-search, is slightly more stable, but it also diverges after few more iterations. On the contrary, for any tested value of M greater than zero, the convergence is reached, showing the greater robustness of the method. Us- f-pref, ρ = 0. 54, q = 1. 0 ANDref, ρ = 0. 39, q = 1. 0 (b) Rate of convergence comparison Figure 13: Convergence slopes of the residual of the power vector for different methods. On the left, the fixed-point method is compared to the Anderson acceleration for different M values (appearing in the legend). In every case, full convergence is imposed and the relaxation is not applied. On the right, the Anderson method with M =5 is compared to the best performing fixed-point algorithm with partial-convergences (α = 0.6, N N = 20, N T H = 40); ρ is the average factor by which the residual decreases and q represents the order of convergence.

ing the Anderson algorithm for the current application, the choice of the relaxation parameter is not anymore crucial for the stability of the method. In Fig. 13-b, the convergence slope of the best fixed-point with partial-convergences (α = 0.6, N N = 20, N T H = 40) is compared to the Anderson method with M=5 and no limitation on the internal iterations. In the x-axis, global iterations are reported, they differ from the multiphysics ones only because of the extraevaluations eventually performed for the line-search. Both the methods show a linear convergence. In the figure, it also appears ρ, an estimate of the rate of convergence. This quantity represents the spectral radius of the equivalent fixed-point function, which for the fixed-point trivially is G α . Anderson convergence rate is better, it needs less global iterations and less multiphysics ones as well. However, as highlighted in the previous sections, this does not necessarily imply a lower computing time. Indeed, even considering the large data manipulation and exchange time (200 s), which significantly overestimates this time cost favouring methods performing less multiphysics iterations, the Anderson algorithm with optimal M -parameter is 8.5 times slower than the best fixed-point with partial-convergences and only 5 % faster than the best fixed-point with full convergence. A further improvement can be obtained by switching off the Scipy extensions, achieving a relative equivalent calculation time of 7.2. To resume, in case of full convergence, the Anderson method has proven the expected superior performance both in terms of robustness and convergence rate (coherently with what found in [START_REF] Toth | Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system[END_REF] modelling a single fuel rod). Nevertheless, if it is compared to most of the partial-convergence fixed-point schemes, its equivalent calculation time is not totally satisfactory.

Partial-convergences implementation

Given the remarkable impact of the partial-convergences on the fixed-point method and the great robustness of the Anderson acceleration, in this subsection, two strategies to apply partial-convergences to the Anderson algorithm are presented.

Algorithm based on the control of single-solver iterations

The algorithm described here is a simple extension of what has been done for the fixed-point algorithm. The partial-convergences are imposed via the limitation of the internal iterations per solver call and to alleviate the inconsistency problem, which may arise because of the different initializations of each iterate, low values of the M -parameter are tested. Indeed, smaller groups of consecutive iterations are expected to be less inconsistent thanks to the more similar initialization they share and to the fact that lower order approximations of the inverse of the Jacobian should be more tolerant to fluctuation of the solutions. This approach is tested for some of the iterations limits of the previous section and for different values of the M -parameter. In all the cases, the Scipy extensions (line-search and regularization parameter) are not activated and the relaxation is not applied. The convergences of the residuals are plotted in Fig. 14 for "20-40" and "100-150".

The results from "20-40" seem to confirm the idea that an M equals to 1 or 2 could be combined with the considered iterations limits, while higher values of M would lead to an unstable algorithm. In respect of "100-150", a different behaviour is observed: it converges only for M >1 and the convergence rate improves up to M =4. In addition, large oscillations of the module of the residual are observed even when the scheme converges. It should be noticed that while the undamped fixed-point with "20-40" converges, "100-150" does not. As written in [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF], the convergence of the Anderson iterations is linked to that of the corresponding fixed-point ones, hence, this element should be kept in mind for further analysis. Another key point is the range of acceptable relaxation factor for the given setting, which could be somehow linked to the tolerance on the correct estimation of the weights. For this reason and to limit the evolution of the solvers initialisations at each multiphysics iteration, the setting "1-1" is tested. The results show that the algorithm does not converge for any value of the M -parameter, which underlines the necessity of having a minimum precision in the function evaluations. With small limits on the singlephysics iterations, the solvers initialisations might hinder the stability of the method. In fact, a potential cause of the oscillations could be found in the fact that the evolution of the solution returned by the Anderson algorithm might be too fast as compared to that of the solver initialisation. For this reason, a different initialisation strategy should be considered to improve the stability of this algorithm.

The time for the data manipulation and exchange is approximately the same of the standard Anderson algorithm. Hence, to simplify the comparisons, a data manipulation and exchange time of 200 s is assumed. The results in terms of relative equivalent calculation time are reported in Table 1. Overall, when the iterations converge, the equivalent calculation time is slightly larger than the equivalent fixed-point with partial-convergences with optimal damping. On the other hand, depending on a non-trivial choice of the M -parameter the iterations may converge or not. Therefore, on this case study, the generalised fixed-point with partial-convergence offers a better compromise.

Algorithm based on the control of the internal tolerances

The considered algorithm deals with the application of the partial-convergences in the Anderson algorithm via the control of the internal convergence criteria.

Modulating this parameter, instead of the single-physics iterations, allows to significantly decrease the dependency on the particular solvers initialisations.

Moreover, in this way, it is easier to bound the error introduced by the convergence limitation. With respect to the literature, [START_REF] Toth | Local improvement results for Anderson acceleration with inaccurate function evaluations[END_REF] provides a mathematical analysis of the convergence of the Anderson algorithm when introducing a bounded error in the iterative function. More into details, in this reference, a similar multiphysics problem is treated (neutronics, thermal-hydraulics and fuel performance on a single fuel rod under steady-state conditions) and the conver-gence of the multiphysics residual norm is studied for different levels of tolerance on the convergence of the internal variables of the specialised solvers. Instead of fixing constant tolerances all along the process, in the proposed method, the convergence process is divided into a sequence of groups of Anderson calculations with increasingly refined imposed precisions on the internal variables. A switching criterion is defined in order to move from one block to the following one, i.e. starting new Anderson iterations with refined internal tolerances. This criterion deals with a check on the minimum number of neutronic iterations per solver call. It should be noticed that this quantity is a measure of the "distance" of the neutronic solution between two consecutive solver calls. Hence, it can be used to avoid to push too far the multiphysic convergence for the degraded internal criteria. Analogously, a criterion could be implemented to loose the tolerance in case of diverging iterations. In the final block, the same multiphysics and internal convergence criteria are imposed. This global scheme is represented in Fig. 15. the thermal-hydraulics and the neutronics. In this example, the algorithm switches to the next block when the number of neutronic power iterations to convergence is lower than ten. In the final block, all the convergence criteria used for the other methods apply here.

An application of this method is reported here. Also in this case, the Scipy extensions are deactivated and there is no damping. A sequence of four tuples of internal precisions is chosen, the values are simply constantly reduced by a factor ten at each step but the first, in which it is divided by a factor fifty, as reported in Table 2. This larger reduction in the first block is simply motivated by the willing to reduce even more the single-solver iterations at the beginning of the process. It should be noticed that for the neutronics the test on the precision is performed only on the last two iterations, while for the thermal-hydraulics this check is set so that it has to be satisfied five consecutive times. Block number:

ε N ε T H 1 5E-2 5E-3 2 1E-3 1E-4 3 1E-4 1E-5 4 1E-5 1E-6
The convergence slopes are available in Fig. 16 for a range of M -parameters.

The fixed-point iterations, which correspond to M =0, are diverging. It is expected that with an even coarser first tuple of convergence criteria it would be possible to obtain a converging algorithm. For M =1, the method is diverging very slowly in the second Anderson block. As proposed before, in case of diverging iterations, the tuple of internal precisions could be loosen to stabilise the method. In particular, an intermediate tuple could be added so that it is in between those required by the second and the first blocks. For larger values of M, almost no oscillations are observed in the convergence of the residuals and the rate of convergence is rather satisfactory. It is observed that the time for data manipulation and exchange does not significantly increase with this algorithm customisation. Therefore, for simplicity, also in this case, this time cost is set to 200 s. The results in terms of relative equivalent calculation time are displayed in Table 3. This algorithm with the considered sequence of convergence criteria and switching criterion converges faster than the standard Anderson method with optimal M -parameter. Anyhow, in terms of equivalent calculation time, it is still significantly slower than the customised Anderson method with optimised iteration limits and M -parameter. On the other hand, it should be noticed that, contrarily to the sequence of internal tolerances and to the switching criterion, the iterations limits are issued by an empirical optimization on this given case study. Moreover, previous studies suggest that the number of iterations strongly depends on the analysed core size and on the operating conditions. Therefore, the algorithm described in this section is potentially more robust than the considered methods based on the control of the single-physics iteration limits. To improve the efficiency of this customised Anderson algorithm, it could be interesting to work on the definition of the internal precision sequence and of the switching criterion. Ultimately, for a broader comparison, the considered numerical methods could be tested for a range of core configurations and operating conditions.

Conclusions

The context of this study is the multiphysics simulation of PWRs under steady-state nominal conditions. The coupling scheme is realised combining a set of specialised solvers. The considered models are two-steps neutronics with pin-cell homogenization, subchannel thermal-hydraulics, explicit 1D heat conduction for every fuel rod and research of the equilibrium fission products.

This coupled problem constitutes a system of nonlinear equations. A simple case study (5x5 PWR mini-core) is defined and it used as applicative test for several numerical methods. In particular, this work is centred on the analysis of the application of the partial-convergences technique to the widespread fixed- In this work, it is studied to apply different limitations on the maximum numbers of both neutronic and thermal-hydraulic iterations per solver call for a wide range of damping factors. An important result is that the application of the partial-convergences to the fixed-point method significantly improves the robustness of the algorithm by strongly reducing its dependency on the relaxation factor. Moreover, for the particular combination of solvers, it appears that also the computing time can be effectively reduced. Ultimately, even if the optimisation of the parameters introduced by this technique is definitely case dependent, the range of values leading to satisfactory results is rather wide.

The application of the Anderson method on the case study has demonstrated its superiority on the standard fixed-point both in terms of robustness and time efficiency. However, for the considered combination of solvers and the given case study, the generalised fixed-point with partial-convergences offers a significantly lower computing time for a wide range of settings. For this reason, two strategies for the integration of the partial-convergences to the Anderson method are studied. The first one, dealing with the limitation of the single-physics iterations, seems not particularly suitable. Considering the same iteration limits, the fixed-point is still faster. Moreover, the choice of the M -parameter becomes more complicated and crucial for the stability of the method. The second strategy consists in controlling the internal convergence criteria of the specialised solvers. More into details, the convergence process is divided into a set of blocks of Anderson iterations with increasingly finer tolerances on the internal variables. Preliminary tests on this customisation of the Anderson algorithm leads to promising results. Future research will deal with a wide range of case studies.

Acknowledgement

The authors are thankful to Simone Santandrea and Davide Mancusi for the numerous fruitful discussions.

Figure 1 :

 1 Figure 1: Level of refinement of the radial meshes, from left to right: pin-cell neutronics, subchannel thermal-hydraulics and heat conduction for every rod. The meshes represent the north-east quarter of a mini-core presented later in the paper, the neutronics includes also the reflector.

λ 2 =

 2 -0. 8 λ 3 = -1. 5

Figure 2 :

 2 Figure 2: Effect of the relaxation factor over the modules of the eigenvalues |λ i |. |λ(α)| > 1 means divergent behaviour (red solid line, Max conv). The original eigenvalues correspond

Figure 3 :

 3 Figure 3: North-east quarter of the loading plan, burnup values relative to the fuel assembly are expressed in MWd/kg.

Figure 4 :

 4 Figure 4: Performance of the Damped fixed-point algorithm for a range of relaxation factors. "NC" stands for Non-Convergent behaviour and "DIV" for DIVergent.

Figure 5 :

 5 Figure 5: Symmetrical reduction of the limits on single-physics iterations, appearing in the tuples format: "N N -N T H ". The partial-convergences technique leads to significant savings in the number of multiphysics iterations and it drastically reduces the dependency on the relaxation factor.

Figure 6 :

 6 Figure 6: Symmetrical reduction of the limits on single-physics iterations, appearing in the tuples format: "N N -N T H ". The partial-convergence significantly improves the performance of the algorithm and drastically reduces the dependency on the relaxation factor. "SC", meaning Slow-Convergence, appears if (teqv,r) > 50.

Figure 7 :

 7 Figure 7: Asymmetrical reduction of the limits on single-physics iterations, appearing in the tuples format: "N N -N T H ". The stability of the algorithm improves more when limiting N N .

Figure 8 :

 8 Figure 8: Asymmetrical reduction of the limits on single-physics iterations, appearing in the tuples format: "N N -N T H ". The stability of the algorithm improves more when limiting N N .

Figure 9 :

 9 Figure 9: Performance analysis of different limits on the maximum number of single-physics iterations per multiphysics call. The same analysis of Fig. 6 and Fig. 8, but assuming a data exchange and manipulation time of 200 s, instead of 2 s.

Figure 10 :

 10 Figure10: Ratio of the total number of single-physics iterations needed to fully solve the coupled problem and a reference iteration number. For the neutronics, this is the number of power iterations to completely solve the neutronics for a temperature and density profile corresponding to flat power. For the thermal-hydraulics, this is the number of iterations necessary to entirely solve the thermal-hydraulics for a flat power profile.

  [-] MAX(ABS()) * 10 3 [%] RMS() * 10 3 [%] [-] MAX(ABS()) * 10 3 [%] RMS() * 10 3 [%] (d) C Xe,eq

Figure 11 :

 11 Figure 11: Convergence proof based on the analysis of the relative discrepancy of each converging setting to the reference values (produced with "inf-inf" with relaxation factor 0.3).

  keff,ref = 1.14129 (b) α = 0.7

Figure 14 :

 14 Figure 14: Convergence of the residual of the total power in the water for Anderson for the two tuples of iterations limits per solver call and for different M -parameters. No line-search, no regularization parameter and no relaxation are applied.

Table 1 :

 1 Relative equivalent calculation times for the generalised fixed-point with partialconvergences, respectively with (α = 0.6, N N = 20, N T H = 40) and (α = 0.6, N N = 100, N T H = 150), Anderson with optimal M -parameter and Anderson with partialconvergences (M = 1, α = 1, N N = 20, N T H = 40) and (M = 4, α = 1, N N = 100, N T H = 150).

FP 0. 6 -

 6 20-40 FP 0.6-100-150 AND 5 AND 1-20-40 AND 4-100

Figure 15 :

 15 Figure 15: Anderson with partial-convergences controlled by the internal precision. A sequence of Anderson calculation blocks characterized by increasingly finer convergence criteria for both

Figure 16 :Table 3 :

 163 Figure 16: Anderson with partial-convergences controlled by the internal precision. The method is tested for a range of M -parameters appearing in the label. The star-markers correspond to the first fixed-point iteration of every Anderson block, to which no calculation of the residual is associated.

  point and Anderson algorithms. This technique deals with the limitation of the convergence process of the thermal-hydraulic and neutronic solvers during each multiphysics iteration. To control the convergence degree two strategies are individuated, one dealing with the limitation of the maximum number of singlephysics iteration per solver call and the other consisting in the modulation of the internal convergence criteria of the specialised solvers. The limitation of the single-physics iteration within the fixed-point algorithm is a common practice for industrial calculations, but systematic studies on this technique are rarely found in literature.

Table 2 :

 2 Sequence defining the progressive refinement of the precision on neutronic and thermal-hydraulic variables. ε N and ε T H respectively refer to the neutron flux and the moderator density convergence criteria.

The order of converge is defined as the value q respecting the following relation lim k→∞x k+1 -x *x k -x * q < M . Linear convergence means q = 1, while quadratic q =

The results appearing in this section are attached to the article under the Mendeley format.

The RMS and maximum of the absolute relative discrepancies associated to each setting are attached to the article as Mendeley Data.

This and more data related to the observed fixed-point bifurcations are attached to the article as Mendeley data.

Data Availability

Datasets related to the performance of each presented fixed-point iteration scheme, the RMS and maximum of the absolute relative discrepancies associated to each fixed-point setting and the residuals convergence in case of fixed-point bifurcations are attached to the submitted article as Mendeley Data.