EFFICIENT MULTIPHYSICS ITERATIONS IN MPACT WITH PARTIALLY CONVERGENT CMFD

Partial convergence of CMFD can help to stabilize multiphysics iteration schemes. In this paper, an efficient multiphysics iteration scheme with near-optimal partially convergent CMFD implemented in MPACT is presented. In the new scheme, the feedback intensity of the problem is automatically estimated, and the relative convergence of CMFD solver is adjusted accordingly. Numerical results show that MPACT with near-optimal partially convergent CMFD can have almost the same convergence rate in problems with feedback as those without feedback. For the problems tested here the run time may be reduced by more than 20% and up to 49% compared with that of current MPACT.


INTRODUCTION
The Coarse Mesh Finite Difference (CMFD) [1] method is a standard Nonlinear Diffusion Acceleration (NDA) method for accelerating numerical transport calculations in reactor physics [2,3,4]. However, CMFD's traditional robustness and fast convergence rate can be significantly degraded when solving large-sized problems [5]. Typical iterative methods use a simple Picard iteration scheme, usually with relaxation [6,7]. In most of these implementations partial convergence is used, but little has been done to determine how much or how little the CMFD problem should be converged in a problem dependent manner. Furthermore, when CMFD is converged more tightly, the multiphysics iteration becomes less stable. This is why relaxation is commonly used. The theoretical understanding of this was detailed in [5]. More recently, our previous research [8] has shown that partially convergent CMFD can stabilize the multiphysics iteration scheme using a relaxation-free method that determines the near-optimal partial convergence of CMFD according to the feedback intensity of the problem. This method was shown to have a theoretical spectral radius comparable with NDA in problems without feedback, and nearly the same behavior for problems with feedback. For this work, the method has been implemented in MPACT [2,9]. In this paper, the implementation details are described and numerical results comparing the performance with the standard relaxed Picard iteration presently used in MPACT for a wide range of multiphysics problems are presented.

Iteration Scheme in MPACT
The flowchart of the original multiphysics iteration scheme in MPACT is shown by the part outside the red box in Figure 1. It is a Picard iteration scheme on the transport iteration level. The pin power distributions are calculated with the latest flux iterates following a neutron transport sweep in MPACT, and sent to other physics solvers. The other physics solvers in turn calculate their corresponding problems. The solutions from the other physics are then sent back to MPACT and utilized to update the cross sections.
The CMFD is utilized to accelerate the overall convergence rate of the transport calculation. The CMFD equation is typically a multigroup (MG) diffusion equation. The default CMFD solver in MPACT solves the full space and energy multigroup problem. A more efficient multilevel CMFD solver, the Multilevel in Space and Energy Diffusion (MSED) solver, has been implemented in MPACT more recently [10]. It is a multilevel method which collapses the multigroup CMFD (MGCMFD) equations to a few-group (typically one group) CMFD operator, and solves the fewgroup CMFD equations to resolve the spatial dependence of the fission source. The one group CMFD (1GCMFD) solutions are used to approximate the fission source of MGCMFD equation, and make the overall problem more tractable. Several multigrid V-cycle iterations between the MGCMFD level and the 1GCMFD level are performed until the convergence criteria of CMFD is met.

Wielandt Shifted Power Iteration
The CMFD equation is an eigenvalue problem. In matrix form, it is written as: In reactor physics, the solutions of interest are the smallest λ and corresponding eigenvector Φ.
The most popular method is the Wielandt Shifted (WS) Power Iteration (PI) scheme. The solution process for this is written as: where l is the index of the power iteration. In practice, the CMFD solution is only partially converged. The algorithm uses a certain Wielandt Shift parameter λ s and L power iterations. It turns out that the partial convergence of CMFD induces equivalent relaxation varying with the frequency of the Fourier error modes. When the partial convergence is near-optimal, these inherent relaxation factors help to stabilize the convergence of slowly oscillating error components that are sources of instability in problem with feedback [8,5].

Near-optimal Partial Convergence
The issue remains on how to estimate the degree of partial convergence in a problem dependent way. Theoretical results from [8] show that the near-optimal partial convergence is given by the Wielandt shift ratio r defined by Equation 3. where Here it is assumed that L is predefined. This is true in practical simulations since the maximum number of power iterations per outer iteration is either specified in the input or has an appropriate default value. Though in most practical simulations, the max iteration limit may not be reached, the power iteration limit L can be enforced to be reached with few modifications in the new implementation. The near-optimal partially convergent CMFD is applicable, because λ s is a natural part in solving the CMFD equations. In Equation 3, the ω 0 is a parameter characterizing the optical thickness of the modeled problem. However, as shown in [8], a fixed ω 0 can be used without degradation of performance. The suggested ω 0 is π 200 . c is the scattering ratio. The most important parameter is γ, which is introduced to characterize the localized feedback intensity [5]. γ is problem-dependent and must be determined by numerical estimation. The expression for γ is It should be noted that the formula is derived for the NDA method. The near-optimal partial convergence determined by the formula results in the same convergence rate in the feedback problem as the NDA in the problem without feedback. In practical simulations, it turns out that the iteration scheme can become even more stable when relaxation technique (power relaxation) is applied together with the near-optimally partially convergent CMFD method.

Feedback Intensity Estimation
According to the derivation of the near-optimal partial convergence formula, the estimated γ should have the following properties. (1) Locality: γ must be calculated without changing the conditions of neighboring cells. The word local here in practical simulation stands for the range inside one coarse mesh cell. Φ here is the localized one group flux. Σ a and Σ f are the localized one group cross sections. Therefore γ is space-dependent.
(2)Non-negativity: The global γ is non-negative. A positive γ denotes the negative feedback for a real reactor, required in nuclear reactor operations. A global γ can be calculated based the localized γ of different cells.(3) Proportionality: In reality, the one group cross sections do not vary much. If only TH feedback is presented, the flux is proportional to the power of the reactor. Therefore, the γ should be nearly proportional to the power for the same model.

As shown in Equation 4
, γ is a term related with the derivative of cross section to the localized flux. The most nature way to calculate these derivatives is to use finite difference method. The idea is perturbing flux first, and then applying the feedback. Thus the derivative can be approximated by: where p stands for perturbation and 0 stands for the unperturbed state. In our method the perturbation is applied in a Red-Black fashion. The fine mesh fluxes of the pin cell with the even index are perturbed by a factor 0.5α%, and those of the pin cell with the odd index are perturbed by −0.5α%. α is a small value to make the Equation 5 hold and make the computational intensity less when calculating the other physics and applying feedback during the perturbation. Here, α = 0.5.
With the Red-Black perturbation, the derivatives of all the coarse mesh cells can be calculated at the same time. Only the γ of the coarse mesh cells where there are fissionable materials will be calculated. The maximum γ is used as the global γ after the feedback intensity of all the cells of interest are calculated. It should be noted that this approach is akin to a numerical estimate of the Jacobian of the multiphysics problem, although its important to note that the procedure results in a very simple estimate of the Jacobian that is ultimately represented by a scalar quantity.
The γ can be calculated at the end of every outer iteration, and the Wielandt shift is adjusted respectively. However, as more computational effort is required by the other physics during the perturbation, calculating γ every outer iteration will make the iteration scheme less efficient. A more practical way is estimating γ at the beginning of 4 th outer iteration if the initial flux guess is random, and 2 th if the initial flux guess is the solution from the previous state. The estimation of γ at the 4 th outer iteration will be justified at subsection 3.1. Before estimating the γ by perturbation, it is suggested that γ is initialized as 0.002 which represents a reasonably large feedback intensity. The flowchart of the overall iteration scheme with near-optimal partial convergent CMFD is shown in Figure 1.

Verification of Estimated γ
Before showing the enhancement of using partially convergent CMFD, the properties of the estimated γ are examined here. The first step is determine which point during the iterations is used to estimate γ. Previously it was stated in subsection 2.4 that the start of 4 th iteration is a reasonable choice. Two problems with thermal-hydraulics (TH) feedback are used here. One is a 3D pin cell problem by extruding the 2D pin cell model specified by the VERA Core Physics Benchmark Problem 1 [11]. The second problem is the single PWR assembly, Problem 6 from [11]. The nominal power of problem 1 is the averaged pin power in the HFP condition. It can be found for both problems simulated at different power levels, the feedback intensity estimated at the 4 th outer iteration is very close to the estimated γ at the last outer iteration. The proportionality of the estimated global γ to the thermal power is also investigated, but is omitted due to the page limit. The conclusion is that the linearly fitting results indicate the estimated global γ is proportional to the power level. This also implies that the near-optimal partial convergence formula should be applicable to these feedback problems.

Performance Comparison
In this subsection, the performance of MPACT with near-optimal partially convergent CMFD is investigated. There are two kinds of near-optimal partially convergent CMFD solvers. One is on the default method in MPACT based on a single level MGCMFD, and is denoted as MGPC. The other is based on MSED, and is denoted as MSEDPC. The performance of MPACT is compared for the default MGCMFD solver, MSED solver and MSEDL solver. MSEDL is a variant of MSED that simply limits the amount of convergence of diffusion solutions of the CMFD solver in a fixed way [10]. This variant was developed because MSED was too efficient at converging the CMFD problem, resulting in less stable multiphysics iterations. For MGPC, the desired number of power iterations is 5, and the Wielandt shift parameter is updated based on the estimated feedback. For MSEDPC, the desired number of total power iterations of 1GCMFD is 20, i.e. 4 V-cycles and 5 1GCMFD power iterations per V-cycle are set. For the MGCMFD, the default setting in MPACT is L = 20 and λ s = 2/3, the partial convergence is almost fixed for the same model. As a result, it is expected that partial convergence of MGCMFD which is the standard option is relatively loose for problems with small feedback and relatively tight in problems with the feedback.

3D Problem 1
For 3D Problem 1, when it is simulated without feedback, the total number of outer iterations needed to converge is 16. Though, it is a simple case, it will become unstable when CMFD is tightly converged. Table 1 shows that for 3D Problem 1 at different power levels, MPACT with MGPC has almost the same total number of outer iterations, indicating that the convergence rate of MPACT with MGPC in problems with feedback is similar with the convergence rate of MPACT in problems without feedback. The P1 is simple and small, therefore MPACT with MGCMFD can have the same outer iteration numbers when power level is smaller than 120%. The run time of MPACT with MGCMFD at low power cases is shorter, because the convergence of MGCMFD is more loose than the predicted convergence in MGPC, reducing the computational intensity solving CMFD equation.

Problem 6
For Problem 6, all the feedback is turn on, i.e. TH feedback and equilibrium xenon feedback is applied, and critical boron search is used. The reference outer iteration is 13 when the problem is   Table 4, it can be found that when default MGCMFD solver is used, it takes more outer iterations for the cases with power no smaller than 120% of the nominal power. MSED and MSEDL without feedback are effective at reducing the total number of outer iterations when power level is low. However, as feedback is applied, MPACT overall converges more slowly and even does not converge with MSED, meaning the CMFD is being converged too tightly. MSEDL, which is used to solve the instability by loosing the convergence criteria, still takes a significant number of outer iterations to converge when the power level is 140%. On the contrary, the near-optimal partially convergent methods (MGPC/MSEDPC) adjusts the convergence of the CMFD solver according to the problem dependent feedback intensity, and allows MPACT to have almost the same number of outer iterations to converge the problem at any power level without the relaxation technique.  The VERA benchmark problem 4 (P4) is used to investigate how MPACT with the near-optimal partially convergent CMFD performs in problems with control rods. Ony TH feedback is applied. The control rod position of the single state case is around 159cm relative to the top of the bottom core plate, or inserted about halfway. The problem is modeled with power level ranging from 20% to 120%. For the multi-state case, which is listed as 100(MS), there are in total 10 states with the control rods withdrawn from the bottom by 23 steps per state. For all the cases listed in this section, no relaxation factor is used.
It can be found that MPACT with near-optimal partially convergent MSED can have a constant 13 outer iterations for all the cases modeled and the run time can be reduced by more than 35% compared with the run time of MPACT with default CMFD solver. The run time for MSEDPC is less than that of MSEDL when power level is low because the outer iteration number for MSEDL is larger than that of MSEDPC. This indicates that in this case MSEDL is too loosely convergent for low power cases but seems reasonably tuned for intermediate power levels. For the multistate cases, The run time is reduced by 38%, close to the time reduction for the single state result.

CONCLUSIONS
In this paper, an efficient multiphysics iteration scheme with near-optimal partially convergent CMFD is implemented in MPACT and presented. In the new scheme, the feedback intensity of the problem is automatically estimated, and the relative convergence of CMFD solver is adjusted accordingly. Numerical results show that MPACT with near-optimal partially convergent CMFD can have almost the same convergence rate in feedback problems as in problems without feedback. The near-optimal partially convergent MSED helps to also solve the slow convergence issue when applying MSED in problem with feedback. With the near-optimal partially convergent CMFD, the run time can be reduced by more than 20% and up to 49% in some cases.