Development of an efficient tightly coupled method for multiphysics reactor transient analysis

Progress in Nuclear Energy 103 (2018) 33–44

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

Development of an efficient tightly coupled method for multiphysics reactor transient analysis

MARK

Jaron P. Senecal, Wei Ji∗ Department of Mechanical, Aerospace, and Nuclear Engineering, Rensselaer Polytechnic Institute, 110 8th St, Troy, NY 12180, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Over-solving Multiphysics Picard iteration

Picard Iteration is a widely used coupling method for multiphysics simulations. This method allows one to directly leverage existing and well-developed single-physics programs without re-writing large portions of the codes. In Picard Iteration, single-physics codes just iteratively pass solutions to each other as inputs until each code has reached a converged solution. However, multiphysics computation linked by Picard Iteration is susceptible to over-solving, which can make the overall computation much less efficient. Over-solving means that each single-physics code provides an accurate solution in each Picard Iteration, which is not necessary in practice. Solving the single-physics codes in an inexact manner, i.e. with relaxed termination criteria, can help avoid this problem. This work develops a modified Picard Iteration coupling method with adaptive, inexact termination criteria for the underlying single-physics codes. Also, nested within the inexact Picard Iteration, inexact Newton methods were applied in the single-physics codes. The effect on the overall computation efficiency due to the inexact (relaxed) termination criteria at both levels is investigated by applying them to solve reactor transient problems. A reactor dynamics problem with temperature feedback in one-dimensional slab geometry is used to scope the behavior of nested inexact solvers. Then these methods are applied to a larger twodimensional Boiling Water Reactor (BWR) problem. Computational time savings reach 55% for the two-dimensional problem. Additionally, applying an inexact termination criterion (inexact Newton method) to each single-physics code results in a further time savings of up to 18%.

1. Introduction Multiphysics modeling and analysis of nuclear reactor core and system designs is one of the current thrusts in nuclear energy research. New computational frameworks and codes have been developed to execute multiphysics simulations incorporating coupled neutronic, thermal-hydraulic, and mechanical behaviors in nuclear reactors (Gaston et al., 2015; Palmtag et al., 2014; Siegel et al., 2007; Magedanz et al., 2015; Chanaron et al., 2015). A multiphysics analysis that accurately incorporates feedbacks from each physics process allows for the best predictions of realistic system behaviors. This is in contrast to conventional single-physics analysis where other physical processes are approximated with fixed input parameters. For the analysis of nuclear reactor transients, which are inherently “multi-physical,” the motivation for developing multiphysics codes is obvious. Rather than using simplified thermal models and linear feedback coefficients, sophisticated high-fidelity codes can be linked together for detailed analysis. Much work is being done along these lines to enable multiphysics reactor simulations. Readers can find many examples in literature (Ivanov

∗

and Avramova, 2007; Bennett et al., 2016; Zerkak et al., 2015; Yilmaz et al., 2017; Kochunas et al., 2017; Mahadevan et al., 2014; Leppänen et al., 2015; Ellis et al., 2017; Herman et al., 2015; Mylonakis et al., 2014). Conceptually, multiphysics coupling can be implemented in two ways. One is the monolithic approach. All the equations that represent different physics are formulated into a single solution scheme. The different single-physics problems are treated as a single problem and are solved simultaneously (Keyes et al., 2012). The coupling of different physical models is implicitly accounted for in the solution scheme. The monolithic approach is suitable for solving strongly coupled multiphysics problems. The other coupling method is called the partitioned approach. Multiple solvers are used and each tackles a different singlephysics problem. These solvers explicitly communicate their answers to each other until the final converged solutions are obtained for all the solvers (Tautges et al., 2011). Each solver has its own solution scheme and is linked to other solvers by a coupling scheme (method). The partitioned approach is much more popular because it allows the reuse of legacy codes (Keyes et al., 2012; Ganine et al., 2013). This is

Corresponding author. E-mail address: [email protected] (W. Ji).

https://doi.org/10.1016/j.pnucene.2017.10.012 Received 10 May 2017; Received in revised form 21 September 2017; Accepted 28 October 2017 0149-1970/ © 2017 Elsevier Ltd. All rights reserved.

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

zero as the simulation progresses (Birken, 2015). He also provided numerical results to support this conclusion; and in all of our example problems the inexact methods arrive at the same solution as the exact methods. In our previous work (Senecal and Ji, 2017, 2015, 2016), several methods were introduced that reduce over-solving. In addition to applying the best of these methods (the Residual Balance method) to nuclear reactor transient problems, further improvements have been developed in order to make it more broadly applicable. A further novel contribution of the present work is to demonstrate the combined benefit of removing over-solving on multiple levels (within both the multiphysics and single-physics solvers). Before proceeding further, several terms are defined for the sake of clarity. “Global iterations” refer to solving each of the single-physics problems and performing the associated data mapping operations. “Constituent iterations” denote the iterations performed to solve a single-physics problem. “Linear iterations” solve for the update step in the Newton method. Throughout this paper the norm operator, ⋅ , refers to the Euclidean norm.

important because the separate legacy codes use specialized methods (Herman et al., 2015) and meshes (Hansen and Owen, 2008), which may differ greatly between the single-physics problems. Using different meshes for each single-physics problem greatly complicates (but does not prevent (Mahadevan et al., 2012)) the utilization of the monolithic approach. Furthermore, not all single-physics problems are strongly coupled to each other, in which case the monolithic approach is superfluous (Gaston et al., 2015). Of the partitioned methods, we focus on Picard Iteration (also referred to as fixed-point iteration (Birken, 2015) or Block Gauss-Seidel (Hamilton et al., 2016)), which iterates between the single-physics problems until the coupled problem has converged. But because each single-physics problem is solved independently and repeatedly, Picard Iteration may suffer from poor performance (Keyes et al., 2012). The poor performance stems from two aspects. One is that Picard Iteration only attains a linear convergence rate. Acceleration schemes can be applied to improve the multiphysics convergence rate and stability. Several schemes have been developed (Macleod, 1986; Walker and Ni, 2011) to address the issue of slow linear convergence rates. The other primary cause of poor performance is generally referred to as over-solving, which is the focus of, and is addressed in, this paper. Over-solving is working to obtain a precise solution to an imprecise problem. For example, there is no need to find an extremely precise temperature distribution until the flux/power distribution has been solved to a commensurate level of precision. However in standard Picard Iteration, each single-physics problem is solved to full precision at every iteration. This is not necessary because the solution keeps changing as a function of the feedback from the other solver(s). Only approximate solutions are needed until all of the solvers approach their final converged solutions (in the last few Picard iterations). Over-solving is well-known in Newton methods where it is not necessary to solve the linear system for the next update to great precision when the current guess is still far from the exact solution (Dembo et al., 1982). Over-solving in Newton solvers has been addressed with a variety of strategies (Dembo et al., 1982; Eisenstat and Walker, 1996; Cai et al., 1994; An et al., 2007; Ter, 2007), referred to as inexact Newton methods. However, over-solving has not been thoroughly treated for multiphysics problems and it presents more complexity than what a single Newton solver faces. In multiphysics simulations using Picard Iteration, there is a plurality of single-physics solvers interacting with each other by iteratively exchanging solutions. Within these Picard iterations, each constituent single-physics solver, depending on its solution method, may have multiple levels of nested iterations—as is the case with Newton-based methods. If so, over-solving may occur within the single-physics solvers as well. Thus over-solving should be dealt with at each level to attain better performance. Over-solving in multiphysics is often addressed only as an afterthought in multiphysics problems (Lipnikov et al., 2013; Clarno et al., 2015; Jareteg et al., 2013), until recently (Birken, 2015). Our previous work (Senecal and Ji, 2017) focused solely on methods for reducing over-solving in partitioned multiphysics problems, here we build upon that progress. The present study focuses on multiphysics problems involving only two single-physics solvers, each implementing a Newton-based method. The basic strategy is to employ inexact methods to relax the termination criteria at each level of iteration. Inexact methods (methods that do not solve the inner level of a nested problem to a tight numerical tolerance) are acceptable because they still arrive at the correct solution. Assuming that the coupled problem converges with standard termination criteria, Birken has shown that the exact answer can still be obtained by using relaxed tolerances in the constituent solvers, provided the tolerance approaches

2. Numerical methods A brief background on Picard Iteration is provided before describing the Residual Balance method. Afterwards, inexact Newton methods are discussed as a means to reduce over-solving in the single-physics solvers. 2.1. Picard Iteration Picard Iteration (PI) is a commonly used method for coupling multiple physics codes. At each time step (steady-state problems have only one “time step”), the algorithm executes the solvers sequentially and iteratively. Because the constituent single-physics codes are executed sequentially, they can be independent codes that have been developed previously (Mylonakis et al., 2014). Fig. 1 portrays the general solution scheme for a Picard Iteration-based multiphysics solver involving two single-physics solvers. Generally, weakly coupled problems require few global iterations to resolve the feedbacks between the

Fig. 1. Schematic diagram of a partitioned multiphysics solver. The “Global” loop refers to the Picard Iteration scheme that connects the single-physics codes. Each single-physics solver is executed in a series of “Constituent” iterations. Finally the “Linear” iterations are nested within the constituent iterations.

34

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Picard Iteration Neutronics Heat Transfer

10 0

Residual Norm

The goal of the Residual Balance (RB) method is to eliminate oversolving without degrading the convergence rate of the coupled problem (under-solving). Its strategy for accomplishing this is to control the solvers such that they keep pace with each other. This requires that the constituent solvers share their statuses with each other. It is assumed that each single-physics solver has a residual vector that describes how accurately the problem's governing equations have been solved. The norm of this vector is a useful metric for determining the solution progress. The residual norm of two constituent solvers (S1 and S2) will be denoted by RiS, j1 and RiS, j2 . The indices i and j refer to the global and constituent iteration, respectively. The Residual Balance method and its recent improvements are described next.

10 -5

10 -10

N 0

HT 2

4

N 6

HT 8

10

N 12

14

HT

N

HT

16

18

20

2.2.1. Residual Balance termination criterion The logic that the Residual Balance method uses to mitigate oversolving is to reduce the residual norm of the constituent solver only just below the initial residual norm of the other solver. In this way each solver will catch up to the progress made by the other without going too far beyond. The termination criterion for the first and second solvers (S1 and S2) are, respectively

Total Constituent Iterations Fig. 2. Residual norm progression during the solution of a time step of the one-dimensional reactor prompt excursion problem described in Section 3.2.2. Each data marker denotes a constituent iteration. Increasing residual norms are due to the feedback from the other single-physics solver. The active solver is denoted along the bottom edge of the figure.

RiS, j1 <

S1 ⎧ a R 0,0 ⎨ a RiS−21,0 ⎩

for i = 0 for i > 0

(1a)

RiS, j2 < a RiS,01 .

(1b)

A user-specified value of a = 0.1 generally works quite well (Lipnikov et al., 2013; Birken, 2015; Senecal and Ji, 2017). This termination criterion is applied alongside standard absolute and relative termination criteria. The standard absolute termination criterion is important because it prevents the solver from stagnating due to roundoff errors. The development leading to the Residual Balance termination criterion is explained further in (Senecal and Ji, 2015, 2017). Balancing the residual norms of the constituent single-physics solvers helps ensure that they work effectively. The solver with the larger residual norm catches up to the other. And if one solver already has a lower residual norm, it waits for the other by doing only a limited amount of work. In this way, Equation (1) helps each solver avoid needless work.

constituent single-physics solvers. As the coupling becomes stronger, the computational costs increase because the feedbacks require more global iterations to resolve. Usually reactor dynamics problems involve strong coupling between neutronics and thermal-hydraulics. The Picard Iteration method is well-understood and converges linearly. With the addition of numerical damping (Hamilton et al., 2013; Toth et al., 2015) or acceleration (e.g. Anderson Acceleration (Hamilton et al., 2016)) schemes, the convergence rate and stability of Picard Iteration can be greatly improved. A notable weakness of Picard Iteration is an inefficiency known as over-solving. The overall accuracy of the simulation is initially quite poor, therefore it is superfluous to solve to a tight tolerance in the early global iterations. Only as the global iterations progress and the multiphysics feedbacks converge is it necessary for the constituent solvers to attain high precision via tight tolerances. Over-solving can be seen in Fig. 2 which shows the solution process for a single time step of the onedimensional reactor prompt reactivity excursion described in Section 3.2.2. The feedback from one solver causes the residual norm of the other to rebound to a higher level. One may surmise that the iterations well below the rebound are unneeded. Thus some inexact method should be used to eliminate these unneeded iterations.

2.2.2. Normalized residual norms One difficulty with multiphysics problems is that the constituent problems may have residuals with differing scales due to differences in physical units or discretization. If the scales differ significantly, the problem with the larger residual will be solved more precisely than the other problem. In order to treat both problems equally, they must be normalized by their respective numerical tolerances. Most iterative numerical solvers use absolute and relative termination criteria to control the magnitude of the residual norm. The standard absolute termination criterion for solver S1 is

2.2. The Residual Balance method Inexact methods, which modify the termination criteria of the constituent solvers, can be used as alternatives to standard Picard Iteration. Previous research (Senecal and Ji, 2017) compared the performance of several inexact methods. Simple methods work well for some problems, but perform poorly for others. For example, the Alternating Nonlinear method (where each single-physics solver does only a single constituent iteration per global iteration) excels at strongly coupled problems but its performance deteriorates for weakly coupled problems. The Residual Balance method is the focus of this work because it is more versatile.

RiS, j1 < εaS1 ,

(2)

and the relative termination criterion is

RiS, j1 < εrS1 RiS,01 .

(3)

εaS1

and the relative The user sets the values of the absolute tolerance tolerance εrS1 (or accepts the program defaults). The norm of Solver 1 at global iteration i and constituent iteration j, RiS, j1 , is normalized to

35

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

S1 Rˆ i, j =

RiS, j1 ε S1

Residual Balance

,

(4)

S1 ε S1 = max(εaS1, εrS1 R 0,0 ).

(5)

S1 ⎧ a Rˆ 0,0

for i = 0

⎨ a Rˆ S 2 i − 1,0 ⎩

for i > 0

Residual Norm

Similar expressions are used to normalize the residual norm of Solver 2. With this normalization, Rˆ i, j is greater than 1.0 until the solver has reached its standard termination criterion. Equation (1) now becomes S1 Rˆ i, j <

10 -10

(6b)

With this normalization, residuals from both solvers can simply be added to compute the total residual norm S1 S2 Ritot = Rˆ i, j + Rˆ i, j .

0

for i > 0 .

for i = 0

⎨ ρi Rˆ S 2 ⎩ 2 i − 1,0

for i > 0

ε S1

⎧a ⎪ < ⎨ ρi ⎪2 ⎩

S1 R 0,0

ε S1 RiS−21,0 ε S2

(9)

. (10a)

The corresponding termination criterion for the second solver (S2)

ε S2

⎧a ⎪ < ⎨ ρi ⎪2 ⎩

S1 R 0,0

ε S1 RiS,01 ε S1

6

7

8

9

10

In this section we turn our focus from over-solving in the global (multiphysics) level to over-solving in the constituent (single-physics) solvers. The constituent solvers in this work use the Jacobian-Free Newton-Krylov (JFNK) method, but the following discussion applies to most Newton-based methods. Newton-based solvers are appealing for single-physics problems because of their asymptotic quadratic convergence rate. Exact Newton methods are rarely used in practice (An et al., 2007; Rheinboldt, 1998). Instead inexact Newton methods are used to mitigate over-solving, caused by a tight termination criterion for the linear system. This results in significantly reduced computational times because fewer linear iterations are required. Newton methods iteratively search for an approximate solution x j that solves a (vector) equation of the form

R (x ) = 0 .

is

RiS, j2

5

2.3. Inexact Newton methods

for i = 0 for i > 0

4

results were more sensitive to the user-supplied parameter. When the single-physics solvers use the termination criterion of Equation (10) it will be referred to as the Adaptive Residual Balance (ARB) method. A listing of the pseudo-code for the ARB method is given in Algorithm 1. Fig. 3 shows the typical behavior of ARB. In this example, ARB does about half as many constituent iterations as PI, shown in Fig. 2.

(8)

S1 ⎧ a Rˆ 0,0

3

Fig. 3. Solution progression with the Adaptive Residual Balance method for one time step of the prompt reactivity excursion described in Section 3.2.2. Note the significant reduction in constituent iterations when compared to Fig. 2.

or

RiS, j1

2

Total Constituent Iterations

The estimate of the multiphysics convergence rate changes with the global iteration i, but since Picard Iteration converges linearly, ρi should come to an approximately constant value. The convergence rate can be used to project the appropriate solver tolerance for the current global iteration. In order to maintain the convergence rate, the active solver's residual norm should be decreased by a factor of at least ρi below the other solver's initial residual norm. In this work, the reduction factor ρi /2 is used instead of a. Equation (6) is so modified to arrive at the final formulation, S1 Rˆ i, j <

1

(7)

2.2.3. Adaptive termination criterion This section describes an adaptive method for setting the value of a based on the problem's convergence rate. The convergence rate of the coupled problem can be estimated by the ratio of the current and previous total residual norm,

Ritot Ritot −1

10 -5

(6a)

S2 S1 Rˆ i, j < a Rˆ i,0 .

ρi =

Neutronics Heat Transfer

10 0

where the effective tolerance is defined as

(11)

The approximate solution is updated by the equation

for i = 0

x j + 1 = x j + Δx j

. for i > 0

for j = 0,1,2… ,

(12)

until R (x j ) is sufficiently close to zero. The update Δx j is found by solving

(10b)

After one global iteration, the estimate for the convergence rate ρi is available. Therefore the user-provided value of a is only needed for the first global iteration. The logic of Equation (10) selects a good termination criterion without need for further user input. This is a notable improvement from previous work (Senecal and Ji, 2015), where the

R'(x j )Δx j = −R (x j ) .

(13)

R'(x j ) is the Jacobian of the residual/function. In JFNK this linear system is solved iteratively by a Krylov method such as GMRES or BiCGSTAB (Knoll and Keyes, 2004). Further details of the mathematical

36

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Gropp, Keyes, and Tidriri selected η = 10−4 . When it is necessary for clarity, the constant forcing term strategy will be denoted by CGKT, as it was in (Ter, 2007). Unless otherwise noted, CGKT is used to solve all of the demonstration problems. An alternative inexact Newton method is Eisenstat and Walker's (EW) method, which specifies the following termination criterion (Eisenstat and Walker, 1996)

procedures used in the JFNK method are available elsewhere (Knoll and Keyes, 2004).

R (x j ) ⎞ ηj = γ ⎛⎜ ⎟ R ⎝ (x j − 1) ⎠

α

for j = 1,2,3… (16)

Because R (x j − 1) is not available on the first iteration ( j = 0 ), the user must specify η0 ∈ (0,1) . The user must also provide γ ∈ [0,1] and α ∈ (1,2]. In this work, the EW method is used with η0 = 0.1, γ = 1, and 1+ 5

(see (Balay et al., 2016; Eisenstat and Walker, 1996)). α= 2 Equation (8) is the multiphysics analogue of Equation (16) with γ = 1 and α = 1. The pseudo-code describing Inexact Newton methods is listed in Algorithm 2. This pseudo-code is a more detailed description of the constituent solver (“For j = 0,1,2, …”) loops in Algorithm 1.

When x j is far from the solution, Δx j will not generally lead to an accurate solution. Therefore Equation (13) does not need to be solved to a tight tolerance in those cases. So to avoid over-solving, inexact Newton methods accept the update step Δx j (or, more precisely, (Δx j )k ) when

R (x j ) + R′ (x j )(Δx j )k ≤ ηj R (x j ) ,

In total, six methods have been presented, they are summarized in Table 1. They can be categorized the termination criterion used in each iteration level. At both the constituent and linear solver levels the termination criteria can be either exact, inexact and fixed, or inexact and adaptive. In the next section, several reactor transient example problems are solved to investigate whether applying inexact methods to multiple levels of a nested solver (RB at the multiphysics level, and EW at the single-physics level) is an efficient scheme. The Adaptive Residual Balance method and the example problems have been implemented in the Multiphysics Object Oriented Simulation Environment (MOOSE) (Gaston et al., 2015). MOOSE is built on top of PETSc (Balay et al., 2016), which implements numerical algorithms such as EW.

(14)

where ηj is known as the “forcing term” (Eisenstat and Walker, 1996). The simplest inexact Newton method is to specify a constant forcing term

ηj = η ,

(15)

where η ∈ (0,1) . In one example (Cai et al., 1994) of this approach, Cai, Table 1 Summary of methods. Iteration Level

Exact

Inexact

Adaptive

Constituent Tolerance

Picard Iteration (PI) based on Equations (2) and (3) Newton based on Equation 13

Residual Balance (RB) based on Equation (6)

Adaptive Residual Balance (ARB) based on Equation (10) Eisenstat and Walker (EW) based on Equation (16)

Linear Tolerance

Inexact Newton (CGKT) based on Equation (14)

37

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

3. Example problems and results

temperature feedback due to Doppler broadening in the one group diffusion equation

We begin with a simplified one-dimensional reactor dynamics problem to demonstrate the behavior of the Residual Balance method. The second example problem is a more realistic two-dimensional reactor. Both problems show significant performance improvement when oversolving is eliminated.

T + 273. 15 − T0 + 273. 15 ⎞ Σa (T ) = Σa (T0) + 0. 01Σa (T0 ) ⎜⎛ ⎟ T0 + 273. 15 ⎝ ⎠ and

D (T ) =

3.1. Slab reactor with fuel temperature feedback

k = 0.1 ⎛ ⎝

cp = 1 +

∂T (x , t ) −∇⋅k (T ) ∇T (x , t )−Ef Σf ϕ (x , t ) = 0. ∂t

(18)

(19)

The slab reactor is 7.5 cm thick with symmetry boundary conditions on one end. At the other end there is no incoming neutron flux and the temperature is held fixed at 100 °C.

T (0) = 100 and ∂ϕ

∂T ∂x

(7. 5) = 0 ∂ϕ (7. ∂x

1

D (T ) ∂x (0) = 2 ϕ (0) and

(20)

5) = 0

(21)

The problem starts from a steady state configuration with a total power of 32.5 W (average power of 2.17 W/cm3 ). The average temperature begins at 410°C. The steady state solution was solved for by progressing the transient problem until the response arrived at a fixed value. The variables and problem parameters are listed in Table 2. In order to create a delayed supercritical (k = 1.00040) system, Σa (T0) = 0.0457766 cm−1 was selected, simulating a control rod movement. The following relations were developed to approximate

Value

Units

Description

ϕ

a

cm−2s−1

Neutron flux

C

a

T v D Σa β ν Σf

a

220000 1/3 a 0.0458a 0.007 2.43 0.024

cm−3 °C cm s-1 cm-1 cm – – cm-1

Temperature Neutron velocity Neutron Diffusion coefficient Macroscopic absorption cross section Delayed neutron fraction Neutrons per fission Macroscopic fission cross section

λ k ρ cp

0.08 0.1a 5a 1a

s-1 W/cmK g/cm3 J/gK

Delayed neutron precursor decay constant Thermal conductivity Density Specific heat

Ef

3.08E-11

J

Energy per fission

a

T − T0  J/gK 500

(24) (25) (26)

3.2. One-dimensional reactor results From steady state, a reactivity of 40 pcm, or 0.057 $, is inserted by changing the absorption cross section to Σa (T0) = 0.0457766 cm−1. The system becomes supercritical instantaneously so there is a 6.1% prompt jump in the flux, followed by a delayed-critical rise. As the flux increases, the fission heating rate increases and Doppler broadening occurs. After 275 s , the flux peaks and the temperature feedback causes the reactor power to decrease. The simulation continues to 700 s in 90 non-uniform time steps. Fig. 4 depicts the average behavior of the slab reactor and Fig. 5 shows the temperature and flux spatial distributions at various times. It is important to note that each solution technique arrives at the same answer. This is because the global termination criteria have the final authority over the precision. That is, the solution for the time step will not be accepted until one of the global termination criteria is met, regardless of what the relaxed constituent solver tolerances may be. The solutions will therefore agree up to the precision of the global tolerance, which is, after all, how the user defines an acceptable error.

Table 2 Problem parameters. Variable

(23)

3.1.1. Implementations notes The Backwards Euler method was selected for time integration. The problem was discretized with linear finite elements. Although the problem is one-dimensional, a two-dimensional 150 × 5 mesh was used in order to simplify visualization and integration over the boundaries. Convergence studies in time and space show that the values of the peak flux and total energy deposition change by less than 0.9% when the time or space grids are refined by a factor of 2. The Residual Balance parameter a was set to 0.1 (Equations (6) and (10)), which usually results in alternating iterations of the constituent solvers. Both solvers used a relative tolerance of εrS1 = εrS 2 = 10−8 . The neutronics solver used an absolute tolerance of εaS1 = 10−12 , whereas the heat transfer solver had an absolute tolerance of εaS2 = 5 × 10−12 . The relative tolerance of the underlying linear solver was η = 10−5 , unless otherwise specified. The solution to a time step was considered converged if both constituent solvers satisfied their termination criteria simultaneously.

(17)

The governing equation for the thermal system is

ρ (T ) cp (T )

T + 273.15 0.41 ⎞  W/cmK 373.15 ⎠

ρ = 5(1 − 22.5 × 10−6 (T − T0 )) g/cm3

1 ∂ϕ (x , t ) − ∇⋅D (T ) ∇ϕ (x , t ) + Σa (T ) ϕ (x , t ) − (1 − β ) ν Σf ϕ (x , t ) v ∂t

∂C (x , t ) −βνΣf ϕ (x , t ) + λC (x , t ) = 0 ∂t

1 1 ⎛ T + 273.15 − T0 + 273.15 ⎞ − ⎜ ⎟ . 3 300 ⎝ T + 273.15 ⎠

The baseline temperature is T0 = 100∘C. This simple model varies with T (compare to Section 3.3 and (Yesilyurt et al., 2012)). The thermal properties k, ρ, and cp each vary with the temperature.

This example problem is a homogeneous, one dimensional reactor with fuel temperature (Doppler) feedback. The system initially has a positive reactivity which is then countered by a negative temperature feedback. A variety of timescales are involved, ranging from the prompt neutron jump to the diffusion of heat through the reactor. This relatively simple problem highlights the value of the Residual Balance method in the context of reactor dynamics. The neutron transport problem is described in simplified terms by the diffusion equation and a single group of delayed neutrons.

− λC (x , t ) = 0

(22)

Delayed neutron precursor concentration

3.2.1. Nested inexact solvers This section compares the performance of the Residual Balance method to Picard Iteration. Then both methods are used with inexact Newton methods. The results listed in Table 3 show that the Adaptive Residual Balance method achieves a time savings of about 35% (1.5x speedup) compared to Picard Iteration in the reactor transient problem. This is

Value varies.

38

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Table 4 Parameter study for the 0.057 $ reactivity insertion problem. Computational time is given relative to Picard Iteration in Table 3. Method Global

Constituent

Table 3 Performance comparison for a 0.057 $ step reactivity insertion. The GMRES linear solver has a default relative tolerance of η = 10−5 . For EW, it begins with a relative tolerance of η0 = 0.1. Values are accumulated over all 90 time steps. Iterations

Computational Time

Global

Constituent

Linear

Picard Iteration ARB

309 335

1108 670

26846 16586

1 0.647

PI + EW ARB + EW

309 516

1746 1366

30920 23284

1.18 0.735

Linear

531

1483

14902

0.733

ARB + CGKT, η = 10−2

346

928

13705

0.601

ARB + CGKT, η = 10−3

343

693

13199

0.557

ARB + CGKT, η = 10−4

335

670

14388

0.579

ARB + CGKT, η = 10−5

335

670

16586

0.647

ARB + EW, η0 = 10−1

516

1366

15730

0.735

ARB + EW, η0 = 10−2

346

929

15038

0.634

ARB + EW, η0 = 10−3

343

691

13183

0.554

ARB + EW, η0 = 10−4

335

670

14388

0.586

ARB + EW, η0 = 10−5

335

670

16586

0.647

accomplished by substantially reducing the number of constituent solver iterations, despite a moderate increase in the number of global iterations. In the second section of Table 3 Eisenstat and Walker's method (EW) is used in an attempt to reduce over-solving in the linear solvers, which use η = 10−5 by default. Surprisingly, PI and ARB both suffer markedly deteriorated performance due to EW with these parameters. The undersolving stems from η0 being too large, so that the linear solver loop terminates before much progress has been made. Because of this, more constituent iterations are required to reach the constituent-level termination criterion. These additional iterations are the primary cause for the extra linear iterations. Another cause is that in later constituent iterations, ηj becomes too small and surplus linear residuals are executed. Eisenstat and Walker used η0 = 1/2 , but it appears that for this problem, where the Newton solvers are wrapped in Picard Iteration, such large values of η0 are not ideal. Since EW did not perform very well, the parameters for the inexact Newton methods will now be systematically varied to explore their effects. In the first part of Table 4, ARB is used in conjunction with CGKT and the value of η is varied between 10−1 and 10−5 (results so far have used η = 10−5). The smallest values of η result in extra linear iterations (over-solving) and slow performance. The large values of η result in excessive constituent iterations (under-solving). But the minimal progress achieved by the imprecisely solved linear steps also interferes with the global convergence rate: large values of ρi in Equation (10) cause under-solving and excessive global iterations. The optimal time is achieved with η = 10−3, because it minimizes the number of linear iterations, which are the primary cost in this problem. Fig. 6 depicts the relationship between linear, constituent, and global iterations and computational time. The initial term of the EW forcing term sequence is varied in the second section of Table 4. Other sources (Eisenstat and Walker, 1996; An et al., 2007; Balay et al., 2016) suggest using η0 on the order of 10−1. However, it seems that in the context of inexact Picard Iteration a lower value should be used. The default of η0 = 10−1, used also in Table 3, is the worst value listed. The best of these options is η0 = 10−3 , which is 1.17x faster than ARB with the baseline value of η = 10−5 . Again, we find that the cost of linear iterations outweighs the cost of constituent iterations (where extra residual evaluations are needed to compute the Jacobian) and global iterations (where data must be mapped between solvers). Both selecting a constant η and using EW achieve similar results with the ARB method in this case, because each constituent (singlephysics) solver is usually doing only one constituent iteration per global iteration. Thus only the first term in the sequence of η’s is ever used, and we are led to the conclusion that the development of a good method for selecting η0 will prove to be worthwhile in future work. Some algorithm

Fig. 5. Temperature and flux distributions at 0, 150, and 275 s.

Method

Computational Time

10−1

ARB + CGKT, η =

Fig. 4. Average power density and average temperature as a function of time.

Iterations

39

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Table 5 Performance comparison for a super-prompt critical example. EW starts with η0 = 10−3 . Method

Iterations

Computational Time

Global

Constituent

Linear

Picard Iteration ARB

379 448

1744 900

37998 14970

1 0.455

PI + EW ARB + EW

380 435

2295 898

31707 11749

0.893 0.377

Table 6 Parameter study for the prompt supercritical problem solved with ARB + EW. The number of linear iterations are listed for various combinations of a and η0 . The best pair is η0 = 10−3 and a = 10−1 .

Fig. 6. Iterations and computer time as a function of η for the 0.057 $ reactivity insertion problem.

that can adaptively find the appropriate value for each solver is desirable because they may have differing behavior. 3.2.2. Prompt reactivity excursion For the next variation of this example problem, 1.04 $ of reactivity (728 pcm) was inserted by setting Σa (T0) = 0.04538 cm−1. Fig. 7 depicts the progression of the transient. The simulation was run to 2.5 s with 100 uniform time steps, which covered the rapid transient period of interest. The performance results are summarized in Table 5. In this strongly coupled problem the Residual Balance method does very well, running 2.2x faster than Picard Iteration. The effect of EW (η0 = 10−3 ) is shown in the second portion of Table 5. This time Picard Iteration does improve somewhat with the use of EW. But for this problem it is evident that over-solving within the multiphysics level is a larger issue than over-solving in the singlephysics solvers, that is, ARB has a larger impact than EW. EW is also quite helpful when applied in combination with ARB. For this strongly coupled problem, EW actually helps decrease constituent and global iterations somewhat. This is likely due to taking smaller steps in the “wrong direction” (Eisenstat and Walker, 1996), which can be an issue when two physics codes are strongly coupled. Adding EW to ARB results in a 1.21x speedup, similar to the benefit seen in the previous problem. Next, the parameters a and η0 are varied to search for insights on the selection of optimal parameters. In Table 6, the number of linear iterations are listed as a proxy for computational time. For this problem, as in the slow transient case, the best value of η0 is again near 10−3 . The best value for a is 10−1, independent of η0 . This suggests that a can be selected first—the user need not worry about the effect of η0 .

η0\\a

10−1

10−2

10−3

10−1

12591

12637

13281

10−2

12904

13861

15194

10−3

11749

12380

13129

10−4

13440

13898

14586

10−5

14764

16032

17258

The results so far can be used to inform a user's selection of inexact methods and parameters. A software user may wish to minimize a certain type of iteration depending on the characteristics of the code and the problem. In order to minimize global iterations standard Picard Iteration is the surest method. But to minimize lower levels, the story is not so simple: under-solving or over-solving can both result in excess iterations. ARB with large a and small η can be selected to keep constituent iterations low. But it may result in under-solving at the multiphysics level (excess global iterations) and over-solving within the single-physics solver (excess linear iterations). Finally, to minimize linear iterations an inexact Newton method such as EW should be applied to the method with approximately the fewest constituent iterations. It appears that the most likely way to get the best overall results is to use ARB + EW with a good selection for η0 . While the majority of the improvement comes from the multiphysics level (ARB), removing oversolving at the single-physics level (using inexact Newton) can also be significant. 3.3. LRA BWR benchmark problem This section describes a classic reactor dynamics problem found in Argonne National Laboratory's Benchmark Problem Book (Argonne Code Center, 1977). The LRA benchmark problem models a control rod drop accident in a Boiling Water Reactor (BWR) which results in a prompt supercritical reactivity excursion. Cross sections for two group diffusion are provided, see Table 7, as well as a model for adiabatic heat up and Doppler feedback. The reactor can be modeled in either two or Table 7 Two-group diffusion cross sections given by the benchmark specification (Argonne Code Center, 1977). Region

Material

Group i

Di (cm)

Σai (cm-1)

νΣfi (cm-1)

Σ1 → 2 (cm-1)

1

Fuel 1 with rod Fuel 1 without rod Fuel 2 with rod Fuel 2 without rod Reflector

1 2 1 2 1 2 1 2 1 2

1.255 0.211 1.268 0.1902 1.259 0.2091 1.259 0.2091 1.257 0.1592

0.008252 0.1003 0.007181 0.07047 0.008002 0.08344 0.008002 0.073324 0.006034 0.01911

0.004602 0.1091 0.004609 0.08675 0.004663 0.1021 0.004663 0.1021 0 0

0.02533

2 3 4 5 Fig. 7. Average power density and temperature after a 1.04 $ step reactivity insertion.

40

0.02767 0.02617 0.02617 0.04754

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Table 9 Specifications for delayed neutron precursors (Argonne Code Center, 1977). λ1 has been corrected according to (Sutton and Aviles, 1996). Group

βi

λi (s-1)

1 2

0.0054 0.001087

0.0654 1.35

modeled with a simple adiabatic assumption

∂T − α (Σf 1ϕ1 + Σf 2 ϕ2) = 0 . ∂t

(31)

The conversion factor α is the energy conversion factor divided by ε the volumetric heat capacity, α = ρc . Doppler feedback is modeled by p

adjusting the fast absorption cross section with the following relation

Σa1 (T ) = Σa1 (T0 )[1 + γ ( T − T0 )],

(32)

where the initial temperature is T0 = 300 K . Finally, to model the ejection of a control blade in region R, the thermal absorption cross section is modified in the first 2 s according to Equation (33).

Σa2 (t ) 1−0. 0606184 t t ≤ 2 =⎧ t>2 ⎨ Σa2 (0) ⎩ 0. 8787631

Fig. 8. Layout and region assignment for the LRA reactor. One control blade is ejected from the region labeled “R”.

three dimensions, in this work we consider the two-dimensional version. Fig. 8 shows the layout of one quarter of the symmetrical reactor. The governing equations are as follows. First, the two-group diffusion equations with an adjustment for axial buckling are used to obtain the flux distribution

1 ∂ϕ1 − ∇D1 ∇ϕ1 + (Σa1 + B2D1 + Σ1 → 2) ϕ1 − ν (1 − β )[Σf 1ϕ1 + Σf 2 ϕ2] v1 ∂t (27)

− λ1 C1 − λ2 C2 = 0 ,

1 ∂ϕ2 − ∇D2 ∇ϕ2 + (Σa2 + B2D2) ϕ2 − Σ1 → 2 ϕ1 = 0 , v2 ∂t

(28)

where functional dependencies have been suppressed for brevity. Table 8 lists the values of several parameters needed for these equations. For both groups, a zero flux boundary condition is applied to the top and right boundaries. Zero current boundary conditions are applied to the bottom and left boundaries, as shown in Fig. 8. Next, the delayed neutron precursor equations for two groups are

∂Ci + λi Ci − νβi (Σf 1ϕ1 + Σf 2 ϕ2) = 0 , ∂t

i = 1,2 .

3.4. Comparison with benchmark results (29) The first test of the problem implementation is to match the eigenvalue for both rods-in and rods-out conditions. Table 10 shows that our values compare well with the values reported in the literature [42, 43, 44]. Eigenvalues are a useful comparison, but they do not expose errors in the transient parameters. The transient results match nicely with previous solutions, which can be found in Table 4 of (Sutton and Aviles, 1996). The time to peak power is 1.432 s, which falls within the other results which range from 1.421 to 1.445 s. The peak power density of 5581 W/cm3 falls within the range 5451–5734 W/cm3. Fig. 9 shows

The delayed neutron precursor constants are provided in Table 9. Next, the reactor power is calculated by

P = ε (Σf 1ϕ1 + Σf 2 ϕ2) ,

(30)

where ε is the thermal energy released per fission. The flux and delayed neutron precursors are scaled so that the initial average power is 1 × 10−6 W/cm3 . The fuel temperature used for Doppler feedback is Table 8 Other data specified in (Argonne Code Center, 1977). The corrected value of γ is listed according to (Sutton and Aviles, 1996).

B 2 = 1.0 × 10−4 ν = 2.43 v1 = 3.0 × 107 cm s-1 v 2 = 3.0 × 105 cm s-1 ε = 3.204 × 10−11 T0 = 300 K

Table 10 Eigenvalue comparison.

axial buckling for both energy groups mean number of neutrons per fission fast neutron velocity thermal neutron velocity Ws/fission energy conversion factor

α = 3.83 × 10−11K cm3

initial temperature conversion factor

γ = 3.034 × 10−3 K-1/2

feedback constant

(33)

Differences in transient parameters occur in the literature. Nearly twenty years after the benchmark was first published, Sutton and Aviles (1996) noted that the parameters used to solve the problem did not match those printed in the benchmark specification. Specifically the delayed neutron precursor decay constant should be λ1 = 0.0654 s−1 instead of 0.00654 s−1, and the Doppler feedback coefficient should be γ = 0.003034 K−1/2 rather than 0.002034 K1/2 [sic]. We implement these corrections in this work, see Tables 8 and 9. To model this problem, the heat transfer equation is solved in a separate application. This allows the neutronics solver to use a fine mesh, while the heat transfer solver uses a coarse mesh that matches the benchmark specification. It receives a power distribution from the neutronics solver and returns an assembly-homogeneous temperature distribution. The two-group diffusion solver uses 64 linear finite elements per fuel assembly for a total of 7744 elements, including the reflector region. The second order backward difference formula is used for time integration. The 3 s transient is solved in 329 time steps as in [44]. Hypre's BoomerAMG preconditioner Henson and Yang, 2002. is used by the neutronics solver.

41

Author

Rods-In

Rods-Out

Langenbuch & Werner (Argonne Code Center, 1977) Finnemann (Argonne Code Center, 1977) Smith (Smith, 1979) Sutton & Aviles (Sutton and Aviles, 1996) This work

0.99633 0.99631 0.99641 0.99637 0.99639

1.01546 1.01531 – – 1.01554

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Average Power Density, W/cm3

1.0E+4

1.0E+3

ANL-7416 This work

1.0E+2 1.0E+1 1.0E+0 1.0E-1 1.0E-2 1.0E-3 1.0E-4 1.0E-5 1.0E-6 0

0.5

1

1.5

2

Time, s

2.5

Fig. 10. Global iterations per time step. The multiphysics coupling is very weak until the power increases substantially. There is some under-solving with ARB + EW, most of which occurs in the weakly coupled time steps near 1.2 s.

Fig. 9. Core-averaged power density as a function of time.

ARB. This is achieved by a substantial reduction in the linear iterations despite a minor increase in the number of constituent and global iterations. Overall, the combination of ARB and EW gets to the correct solution 2.68 times faster than PI. Fig. 10 shows the global iterations per time step for Picard Iteration and ARB + EW. The problem starts out very weakly coupled because at such low power the feedback effects are negligible. ARB + EW incurs most of its extra global iterations between 1.0 and 1.2 s, where it does not handle the weak coupling (feedback) quite as well as PI. As the power increases, and especially at the two power peaks, the coupling (feedbacks) gets much stronger and more global iterations are required. Abrupt changes in the time step size cause peaks in the number of global iterations. This effect shows up in the following figure, which generally follows the trend of this figure. Fig. 11 demonstrates the power of the Adaptive Residual Balance method to reduce over-solving at the multiphysics level. At nearly every time step, the number of constituent iterations needed is substantially less for ARB + EW than for PI. As expected from comparing Tables 5 and 3, and the findings in (Senecal and Ji, 2017), ARB has the most significant effect while the physics are strongly coupled. A substantial reduction in the number of linear iterations due to the use of ARB + EW is shown in Fig. 12. The preponderance of the computational cost is associated with the linear iterations because they are so numerous. Fig. 12, then, is useful for visualizing where performance has improved the most. Most of this benefit corresponds with the reduction in the number of constituent iterations. Fig. 13 compares the

that the transient average power results of this work compare well with other solutions. Now that we have shown that the problem is being solved correctly, we can move on to the purpose of the example problem: comparing the performance of the various solution methods.

3.5. LRA performance results The LRA BWR benchmark problem has performance characteristics similar to the one-dimensional reactor problem. That is, the Residual Balance method can be used for a substantial performance improvement and adding EW to control the linear tolerance results in further improvement. While each level has a certain amount of overhead computational cost (cost of entering a loop, as opposed to cost per iteration), the majority of the cost is due to the linear iterations because they are so much more numerous than constituent iterations—approximately 30 linear iterations per constituent iteration. The best solution method will therefore minimize the number of linear iterations and in so doing keep global and constituent iterations low as well. Table 11 shows the results of the LRA transient problem solved with several methods chosen according to the one-dimensional reactor problem results. First, we assess the effect of adding EW to PI. Surprisingly, there is no performance benefit. In fact, the under-solving of EW causes additional constituent iterations which lead to additional linear iterations as well, which finally results in an increase in computational cost of 15%. Even when η0 is set to a better value (10−3 ), EW still does not provide any performance benefit. Next, we compare PI to ARB. For this problem, over-solving on the multiphysics level is a significant issue: removing it with ARB reduces the wall-clock time by 55% (2.2x speedup). This is achieved even with an increase in the number of global iterations. Finally, ARB + EW turns out to be 1.22x faster than Table 11 Performance comparison for the LRA benchmark problem with inexact solvers at various levels. Method

Iterations

Picard Iteration PI + EW, η0 = 10−1 PI + EW, η0 = 10−3 ARB ARB + EW, η0 =

10−3

Computational Time

Global

Constituent

Linear

843 843

3053 4184

125155 139156

843

3320

126239

1.011

937 1001

1901 2030

55655 44143

0.454 0.373

1 1.15 Fig. 11. Neutronics solver constituent iterations per time step. ARB + EW successfully reduces over-solving at the multiphysics level, resulting in many fewer constituent iterations.

42

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

(Eisenstat and Walker, 1996) resulted in an additional speedup of about 1.2x, for a total speedup of 2.68x. For the problems considered here, where the inexact Newton method is nested inside of inexact Picard Iteration, the best starting value for the forcing term sequence is on the order of η0 = 10−3 . This stands in contrast to the values on the order of 0.1 suggested elsewhere (Eisenstat and Walker, 1996; Balay et al., 2016) for single-physics applications. In order to select the best value for η0 , it will be a goal of future work to implement a method that adaptively selects the linear tolerance based on the constituent tolerance. Acknowledgments The first author was supported by a Department of Energy, Office of Nuclear Energy, Integrated University Program Graduate Fellowship, and the Nuclear Regulatory Commission Fellowship Program under the Grant NRC-HQ-13-G-38-0035.

Fig. 12. Neutronics solver linear iterations per time step. The inexact method does many fewer linear iterations and much less computational work overall.

References An, H.-B., Mo, Z.-Y., Liu, X.-P., 2007. A choice of forcing terms in inexact newton method. J. Comput. Appl. Math. 200 (1), 47–60. https://doi.org/10.1016/j.cam.2005.12.030. Argonne Code Center, 1977. Benchmark Problem Book. Tech. Rep. ANL-7416, Supplement 2. Argonne National Laboratory. Balay, S., Abhyankar, S., Adams, M.F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W.D., Kaushik, D., Knepley, M.G., McInnes, L.C., Rupp, K., Smith, B.F., Zampini, S., Zhang, H., Zhang, H., 2016. PETSc Users Manual. Tech. Rep. ANL-95/11-Revision 3.7. Argonne National Laboratory. http://www.mcs.anl. gov/petsc. Bennett, A., Avramova, M., Ivanov, K., 2016. Coupled MCNP6/CTF code: development, testing, and application. Ann. Nucl. Energy 96, 1–11. https://doi.org/10.1016/j. anucene.2016.05.008. Birken, P., 2015. Termination criteria for inexact fixed-point schemes. Numer. Linear Algebra Appl. 22 (4), 702–716. http://dx.doi.org/10.1002/nla.1982. Cai, X.-C., Gropp, W.D., Keyes, D.E., Tidriri, M.D., 1994. Newton-Krylov-Schwarz Methods in CFD. Vieweg+Teubner Verlag, Wiesbaden, pp. 17–30. http://dx.doi.org/ 10.1007/978-3-663-14007-8_3. Chanaron, B., Ahnert, C., Crouzet, N., Sanchez, V., Kolev, N., Marchand, O., Kliem, S., Papukchiev, A., 2015. Advanced multi-physics simulation for reactor safety in the framework of the NURESAFE project. Ann. Nucl. Energy 84, 166–177. https://doi. org/10.1016/j.anucene.2014.12.013. Clarno, K.T., Pawlowski, R., Evans, T.M., Kochunas, B., Turner, J., 2015. High fidelity modeling of pellet-clad interaction using the CASL Virtual Environment for Reactor Applications. In: Proceedings of the ANS Joint International Conference on Mathematics and Computation (M&C 2015), Supercomputing in Nuclear Applications (SNA) and the Monte Carlo (MC) Method, Nashville, TN, USA. Dembo, R.S., Eisenstat, S.C., Steihaug, T., 1982. Inexact Newton methods. SIAM J. Numer. analysis 19 (2), 400–408. http://dx.doi.org/10.1137/0719025. Eisenstat, S.C., Walker, H.F., 1996. Choosing the forcing terms in an inexact Newton method. SIAM J. Sci. Comput. 17, 16–32. http://dx.doi.org/10.1137/0917003. Ellis, M., Gaston, D., Forget, B., Smith, K., 2017. Preliminary coupling of the Monte Carlo code OpenMC and the multiphysics Multiphysics Object-Oriented Simulation Environment for analyzing doppler feedback in Monte Carlo simulations. Nucl. Sci. Eng. 185 (1). http://doi.org/10.13182/NSE16-26. Ganine, V., Hills, N.J., Lapworth, B.L., 2013. Nonlinear acceleration of coupled fluidstructure transient thermalproblems by anderson mixing. Int. J. Numer. Methods Fluids 71 (8), 939–959. http://dx.doi.org/10.1002/fld.3689. Gaston, D.R., Permann, C.J., Peterson, J.W., Slaughter, A.E., Andrš, D., Wang, Y., Short, M.P., Perez, D.M., Tonks, M.R., Ortensi, J., et al., 2015. Physics-based multiscale coupling for full core nuclear reactor simulation. Ann. Nucl. Energy 84, 45–54. http://dx.doi.org/10.1016/j.anucene.2014.09.060. Hamilton, S., Clarno, K., Berrill, M., Evans, T., Davidson, G., Lefebvre, R., Sampath, R., Hansel, J., Ragusa, J., Josey, C., 2013. Multiphysics simulations for LWR analysis. In: Proceedings of the 2013 International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering - M and C. Hamilton, S., Berrill, M., Clarno, K., Pawlowski, R., Toth, A., Kelley, C., Evans, T., Philip, B., 2016. An assessment of coupling algorithms for nuclear reactor core physics simulations. J. Comput. Phys. 311, 241–257. http://dx.doi.org/10.1016/j.jcp.2016.02. 012. Hansen, G., Owen, S., 2008. Mesh generation technology for nuclear reactor simulation; barriers and opportunities. Nucl. Eng. Des. 238 (10), 2590–2605. http://dx.doi.org/ 10.1016/j.nucengdes.2008.05.016. Henson, V.E., Yang, U.M., 2002. BoomerAMG: a parallel algebraic multigrid solver and preconditioner. Appl. Numer. Math. 41 (1), 155–177. https://doi.org/10.1016/ S0168-9274(01)00115-5. Herman, B.R., Forget, B., Smith, K., 2015. Progress toward Monte Carlo–thermal hydraulic coupling using low-order nonlinear diffusion acceleration methods. Ann. Nucl. Energy 84, 63–72. http://dx.doi.org/10.1016/j.anucene.2014.10.029. Ivanov, K., Avramova, M., 2007. Challenges in coupled thermal–hydraulics and neutronics simulations for LWR safety analysis. Ann. Nucl. Energy 34 (6), 501–513.

Fig. 13. Neutronics solver linear iterations per time step for ARB and ARB + EW methods. EW reduces over-solving at the single-physics level, making the simulation even more efficient.

linear iterations for ARB and ARB + EW, and the advantageous effect of using EW is clearly seen. With the parameter η0 set to 10−3 , ARB + EW requires fewer linear iterations at nearly every time step. As in the one-dimensional reactor transient problems, the performance benefit of using inexact solvers is substantial. The main benefit comes from the reduced number of constituent iterations made possible by the Residual Balance method. The 2.2x speedup can be further augmented with a 1.2x speedup by adding an inexact Newton method. 4. Conclusion Several improvements made to the Residual Balance method were introduced in Section 2.2. A special method of normalization allows it to handle constituent solvers that have greatly differing scales. Furthermore, with an estimate of the convergence rate, the performance of ARB is less sensitive to potentially poor user input. Also, ARB can adapt to better handle more weakly coupled problems. Together these improvements allow the ARB to perform well for a larger range of problems. This work has investigated ways to reduce the computational burden of tightly coupled multiphysics reactor transient simulations. The foremost conclusion is that over-solving the single-physics constituents of the model can lead to significant performance losses. The Residual Balance method described in this work avoids over-solving which results in wall-clock savings of up to 55% (about 2.2x speedup). The effect of using an inexact Newton method to control the linear tolerance was also investigated. Eisenstat and Walker's method 43

Progress in Nuclear Energy 103 (2018) 33–44

J.P. Senecal, W. Ji

Physics (PHYSOR 2014), pp. 1–26. Rheinboldt, W.C., 1998. Combinations of processes. Soc. Industrial Appl. Math. http://dx. doi.org/10.1137/1.9781611970012.ch6. Senecal, J.P., Ji, W., 2015. Comparison of novel multiphysics coupling methods in MOOSE. Trans. Am. Nucl. Soc. 113, 707–710. Senecal, J.P., Ji, W., 2016. A novel tightly coupled method for reactor transient simulations. In: Transactions of the American Nuclear Society Proceedings of the International Topical Meeting on Advances in Reactor Physics (PHYSOR 2016), Sun Valley, Idaho, May 1–5. Senecal, J.P., Ji, W., 2017. Approaches for mitigating over-solving in multiphysics simulations. Int. J. Numer. Methods Eng. 112 (6), 503–528. http://dx.doi.org/10. 1002/nme.5516. Siegel, A., Tautges, T., Caceres, A., Kaushik, D., Fischer, P., Palmiotti, G., Smith, M., Ragusa, J., 2007. Software design of SHARP. In: Joint Int'l Topical Mtg. On Mathematics & Computation and Supercomputing in Nuclear Applications, Monterey, California, April 15–19. Smith, K.S., 1979. An Analytic Nodal Method for Solving the Two-group, Multidimensional, Static and Transient Neutron Diffusion Equations. Ph.D. thesis. Massachusetts Institute of Technology. Sutton, T., Aviles, B., 1996. Diffusion theory methods for spatial kinetics calculations. Prog. Nucl. Energy 30 (2), 119–182. http://dx.doi.org/10.1016/0149-1970(95) 00082-U. Tautges, T.J., Kim, H., Caceres, A., Jain, R., 2011. Coupled multi-physics simulation frameworks for reactor simulation: a bottom-up approach. In: Int. Conf. on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2011), Rio de Janiero, Brazil, 8–12 May 2011. Ter, T.-P., 2007. A Problem-solving Environment for the Numerical Solution of Nonlinear Algebraic Equations. Ph.D. thesis. University of Saskatchewan Saskatoon. Toth, A., Kelley, C., Slattery, S., Hamilton, S., Clarno, K., Pawlowski, R., 2015. Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system. In: Joint International Conference on Mathematics and Computation (M&C), Supercomputing in Nuclear Applications (SNA), and the Monte Carlo (MC) Method. Walker, H.F., Ni, P., 2011. Anderson acceleration for fixed-point iterations. SIAM J. Numer. Analysis 49 (4), 1715–1735. http://dx.doi.org/10.1137/10078356X. Yesilyurt, G., Martin, W.R., Brown, F.B., 2012. On-the-fly doppler broadening for monte carlo codes. Nucl. Sci. Eng. 171 (3), 239–257. http://dx.doi.org/10.13182/ NSE11-67. Yilmaz, M.O., Avramova, M.N., Andersen, J.G., July 2017. Multi-physics code system with improved feedback modeling. Prog. Nucl. Energy 98, 94–108. http://dx.doi.org/ 10.1016/j.pnucene.2017.03.007. Zerkak, O., Kozlowski, T., Gajev, I., 2015. Review of multi-physics temporal coupling methods for analysis of nuclear reactors. Ann. Nucl. Energy 84, 225–233. http://dx. doi.org/10.1016/j.anucene.2015.01.019.

http://dx.doi.org/10.1016/j.anucene.2007.02.016. Jareteg, K., Vinai, P., Demazière, C., 2013. Investigation of the possibility to use a finemesh solver for resolving coupled neutronics and thermal-hydraulics. In: International Conference on Mathematics and Computational Methods Applied to Nuclear Science & Engineering (M&C 2013) Sun Valley, Idaho, USA, May 5-9, 2013, on CD-ROM, vol. 3. pp. 2002–2013. Keyes, D., McInnes, L., Woodward, C., et al., 2012. Multiphysics simulations: challenges and opportunities. Int. J. High Perform. Comput. Appl. 27, 4–83. http://dx.doi.org/ 10.1177/1094342012468181. Knoll, D.A., Keyes, D.E., 2004. Jacobian-free newton-krylov methods: a survey of approaches and applications. J. Comput. Phys. 193, 357–397. http://dx.doi.org/10. 1016/j.jcp.2003.08.010. Kochunas, B., Collins, B., Stimpson, S., Salko, R., Jabaay, D., Graham, A., Liu, Y., Kim, K.S., Wieselquist, W., Godfrey, A., Clarno, K., Palmtag, S., Downar, T., Gehin, J., 2017. VERA core simulator methodology for pressurized water reactor cycle depletion. Nuclear Science and Engineering 185 (1), 217–231. http://dx.doi.org/10. 13182/NSE16-39. Leppänen, J., Hovi, V., Ikonen, T., Kurki, J., Pusa, M., Valtavirta, V., Viitanen, T., 2015. The numerical multi-physics project (NUMPS) at VTT technical research centre of Finland. Ann. Nucl. Energy 84, 55–62. http://dx.doi.org/10.1016/j.anucene.2014. 10.014. Lipnikov, K., Svyatskiy, D., Vassilevski, Y., 2013. Anderson acceleration for nonlinear finite volume scheme for advection-diffusion problems. SIAM J. Sci. Comput. 35, A1120–A1136. http://dx.doi.org/10.1137/120867846. Macleod, A.J., 1986. Acceleration of vector sequences by multi-dimensional Δ2 methods. Commun. Appl. Numer. Methods 2 (4), 385–392. http://dx.doi.org/10.1002/cnm. 1630020409. Magedanz, J., Avramova, M., Perin, Y., Velkov, A., 2015. High-fidelity multi-physics system TORT-TD/CTF/FRAPTRAN for light water reactor analysis. Ann. Nucl. Energy 84, 234–243. http://dx.doi.org/10.1016/j.anucene.2015.01.033. Mahadevan, V.S., Ragusa, J.C., Mousseau, V.A., 2012. A verification exercise in multiphysics simulations for coupled reactor physics calculations. Prog. Nucl. Energy 55 (Suppl. C), 12–32. https://doi.org/10.1016/j.pnucene.2011.10.013. Mahadevan, V.S., Merzari, E., Tautges, T., Jain, R., Obabko, A., Smith, M., Fischer, P., 2014. High-resolution coupled physics solvers for analysing fine-scale nuclear reactor design problems. Phil. Trans. Roy. Soc. Lond. A Math. Phys. Eng. Sci. 372 (2021). http://dx.doi.org/10.1098/rsta.2013.0381. Mylonakis, A.G., Varvayanni, M., Catsaros, N., Savva, P., Grigoriadis, D.G.E., 2014. Multiphysics and multi-scale methods used in nuclear reactor analysis. Ann. Nucl. Energy 72, 104–119. http://dx.doi.org/10.1016/j.anucene.2014.05.002. Palmtag, S., Clarno, K., Davidson, G., Salko, R., Evans, T., Turner, J., Schmidt, R., 2014. Coupled neutronics and thermal-hydraulic solution of a full-core PWR using VERACS. In: Proceedings of the International Topical Meeting on Advances in Reactor

44