Jacobian-free Newton-Krylov Method Research Articles

Predictor–Corrector Nodal Integral Method for simulation of high Reynolds number fluid flow using larger time steps in Burgers ’ equation

Nodal Integral Methods (NIM), although frequently used for solving neutron transport equations, have not found wider acceptance for solving fluid flow problems. One of the pivotal reasons behind this is the lack of efficient non-linear solvers for the system of coupled algebraic equations obtained by these methods, thus making such methods prohibitively expensive for highly nonlinear (higher Reynolds number) fluid flow problems. Successive improvements in earlier attempts resulted in an improved version called Modified NIM (MNIM) scheme; following which, a further modified version appeared (called as M2NIM), that used the concept of the delayed coefficients. Upon comparing the solutions obtained using MNIM and M2NIM schemes, respectively, it was observed that although the modified scheme with delayed coefficients (M2NIM) has significantly faster convergence, it is less accurate. The convergence and accuracy are very likely to suffer to a significant extent in the high Reynolds number fluid flows, especially when a large time-step is to be considered. In order to resolve these difficulties, a new type of physics-based predictor–corrector algorithm is proposed in the present study. In the proposed computational algorithm, the linearized M2NIM scheme is used to predict the solution, whereas the MNIM scheme is utilized to improve the predictive guess. The novelty of such a hybrid numerical algorithm comes from the fact that it combines the advantage of the faster convergence of M2NIM with the accuracy of the MNIM scheme. The algorithm also benefits from the unique inclusion of Jacobian-free Newton-Krylov (JFNK) method, which helps in getting rid of the formation of large Jacobian matrices, thereby reducing unnecessary computational overhead. The proposed methodology has been applied to solve a non-linear convection–diffusion problem, represented by the Burgers’ equation. The computational results for bothone-dimensional and two-dimensional Burgers’ equation are presented to demonstrate the effectiveness of the developed novel algorithm, as well as, the advantage it offers over the existing numerical methods.

The JFNK method for the PWR's transient simulation considering neutronics, thermal hydraulics and mechanics

Abstract A new task of using the Jacobian-Free-Newton-Krylov (JFNK) method for the PWR core transient simulations involving neutronics, thermal hydraulics and mechanics is conducted. For the transient scenario of PWR, normally the Picard iteration of the coupled coarse-mesh nodal equations and parallel channel TH equations is performed to get the transient solution. In order to solve the coupled equations faster and more stable, the Newton Krylov (NK) method based on the explicit matrix was studied. However, the NK method is hard to be extended to the cases with more physics phenomenon coupled, thus the JFNK based iteration scheme is developed for the nodal method and parallel-channel TH method. The local gap conductance is sensitive to the gap width and will influence the temperature distribution in the fuel rod significantly. To further consider the local gap conductance during the transient scenario, a 1D mechanics model is coupled into the JFNK scheme to account for the fuel thermal expansion effect. To improve the efficiency, the physics-based precondition and scaling technique are developed for the JFNK iteration. Numerical tests show good convergence behavior of the iterations and demonstrate the influence of the fuel thermal expansion effect during the rod ejection problems.

The comparison between nonlinear and linear preconditioning JFNK method for transient neutronics/thermal-hydraulics coupling problem

Jacobian-free Newton Krylov (JFNK) method is an advanced numerical method for solving the transient multi-physics coupling problem in nuclear reactor, where all the coupled equations are simultaneously solved in a tightly nonlinear form. Constructing a high performance preconditioning is a key issue for JFNK method to efficiently solve this complex coupled system. Two families of preconditioning method, linear preconditioning and nonlinear preconditioning, have been proposed and developed independently. Linear preconditioning is a transformation of linear equations derived from the Newton linearization. Nonlinear preconditioning is also an attractive method because it can be easily implemented as a black box coupling. The comparison between nonlinear and linear preconditioning is made and validated by numerical tests using a two-dimension simplified reactor model in this work. The results show that the similar convergence rate is achieved for nonlinear and linear preconditioning in the transient neutronics/thermal-hydraulics coupling problem. However, the computational efficiency of linear preconditioning is always higher than that of nonlinear preconditioning, because an extra inverse of the preconditioner should be calculated per Krylov iteration in nonlinear preconditioning. When the computational cost of the inverse preconditioner is usually very expensive, like using the original engineering code as the preconditioning, linear preconditioning has an absolute advantage. While compared with the nonlinear preconditioning, a relatively large burden of code development has to be paid for linear preconditioning.

The Gauss-Siedel (GS) and the Jacobian-free Newton-Krylov (JFNK) methods applied to steady 2D buoyancy convection in vertical cylinders

The Gauss-Siedel (GS) and the Jacobian-free Newton-Krylov (JFNK) methods applied to steady 2D buoyancy convection in vertical cylinders

Reacting flow analysis of a cavity-based scramjet combustor using a Jacobian-free Newton–Krylov method

ABSTRACTThe scramjet is a rather a new technology and there are many issues related to their operation, especially when it comes to the combustion processes. Combustion in high-speed flows causes various problems such as flame instability and poor fuel–air mixing efficiency. One of the methods used to overcome these problems is to recess a cavity in the combustor wall where a secondary flow is generated. In this study, a computational fluid dynamics (CFD) code is developed to analyse the reacting flow passing through the cavity-based scramjet combustor. The developed code is based on three-dimensional coupled Navier–Stokes and finite rate chemistry equations. An ethylene-air reduced chemical reaction model is used as a fuel–air combination. The Spalart–Allmaras model is utilised for turbulence closure. The non-dimensional form of the flow and chemical reaction equations are discretised using a finite volume method. The Jacobian-Free Newton–Krylov (JFNK) method is used to solve the coupled system of non-linear equations. The JFNK is a matrix-free solution method which improves the computational cost of Newton’s method. The parameters that affect the performance of the JFNK method are studied in the analysis of a scramjet combustor. The influence of the forcing term on the convergence of the JFNK method is studied in the analysis of scramjet combustor. Different upwind flux vector splitting methods are utilised. Various flux limiter techniques are employed for the calculations of higher order flux vectors. The effects of flux vector splitting and flux limiter methods on the convergence and accuracy of the JFNK method are evaluated. Moreover, the variations of the mixing efficiency with fuel injection angles are studied.

Coupling Methods for HTR-PM Multicircuit Problem

The High-Temperature Gas-Cooled Reactor–Pebble Bed Module (HTR-PM) is a large-scale complex system that includes reactor core, steam generator, helium circulator, and other important components. When integrating these components, coupling problems such as multiphysics problem, multicircuit problem, multiscale problem, and multimodule problem arise in the numerical simulation. The HTR-PM multicircuit system comprises the primary circuit and secondary circuit, which are simulated by two independent codes and coupled by the interface in the once-through steam generator. Although time-consuming, Picard iteration is a feasible and convenient coupling method to integrate two components because oversolving in the early stages of the iteration causes strong fluctuation between circuits. To address this problem, optimization of the maximum subiteration number and convergence precision have been implemented to improve the efficiency and numerical stability of the Picard iteration. The Dynamic Residual Balance method, an improved version of the Residual Balance method, is proposed to prevent oversolving inside the subiterations. It takes into consideration fluctuation between circuits, and this method is robust in a wide range of cases. Moreover, the nonlinear preconditioned Jacobian-Free Newton-Krylov method, which has less fluctuation between circuits than Picard iteration, is a coupling scheme that updates all the solution variables from the primary circuit and the secondary circuit simultaneously. Outstanding convergence and efficiency can be obtained by implementing the proper preconditioner in this HTR-PM multicircuit problem. The downside is that it requires significant modification to the legacy codes.

Application of implicit Roe-type scheme and Jacobian-Free Newton-Krylov method to two-phase flow problems

A new implicit second-order accurate spatial scheme for steady-state thermal-hydraulic simulations of the two-phase two-fluid six-equation model is proposed. The new scheme is based on a Roe-type numerical flux that is formulated with the help of analytical approximate eigenvalues and eigenvectors of the two-phase system. Approximate eigenvalues and eigenvectors are obtained with a structured Jacobian matrix that is general for arbitrary Equation of State (EOS). Three issues are solved in this article: (1) A fully-implicit scheme using the backward Euler method is proposed to improve the stability and avoid the time step limit of the originally explicit scheme. The fully implicit scheme is solved with the Jacobian-Free Newton-Krylov method. (2) A second-order accurate spatial scheme is proposed by converting the first-order Roe-type numerical flux to a second-order one. The conversion is made by extending an existing procedure for single-phase flows to two-phase flows. (3) Phase appearance issue is treated with a numerical procedure, which requires little modification to the schemes. The new scheme is verified and validated with two numerical tests: the faucet flow problem and the BFBT benchmark.

An Assessment of Coupling Algorithms in HTR Simulator TINTE

This paper evaluates the performance of neutronic and thermal-hydraulic coupling algorithms in transient problems based on the high-temperature gas-cooled reactor simulator TINTE. In particular, the operator splitting semi-implicit (OSSI), Picard iteration, and Jacobian-free Newton-Krylov (JFNK) methods are compared by a practical engineering model. The OSSI method is employed in the original TINTE. The fully implicit algorithms TINTE-Picard and TINTE-JFNK are implemented in this study. Several special numerical technologies are discussed to improve the performance of JFNK. First, a novel JFNK variant is employed to deal with the multiscale coupling between local fuel sphere temperature and global solid porous media temperature. Second, the preconditioning strategy is determined by making a balance between performance and code burden. Finally, the scaling modifications of the Jacobian matrix and perturbation size are investigated to solve the ill-posed problem. What is more, the framework of TINTE-Picard and TINTE-JFNK is presented, and the key points of implementation are discussed. Numerical results indicate that the advanced coupling algorithms Picard and JFNK can achieve higher computational performance than the original semi-implicit coupling algorithm in TINTE due to the accuracy and stability advantage.

Jacobian-Free Newton-Krylov Nodal Expansion Methods with Physics-Based Preconditioner and Local Elimination for Three-Dimensional and Multigroup k-Eigenvalue Problems

Motivated by the high accuracy and efficiency of nodal methods on the coarse meshes and the superlinear convergence and high efficiency of Jacobian-free Newton-Krylov (JFNK) methods for large-scale nonlinear problems, a new JFNK nodal expansion method (NEM) with the physics-based preconditioner and local elimination NEM_JFNK is successfully developed to solve three-dimensional (3D) and multigroup k-eigenvalue problems by combining and integrating the NEM discrete systems into the framework of JFNK methods. A local elimination technique of NEM_JFNK is developed to eliminate some intermediate variables, expansion coefficients, and transverse leakage terms through equivalent transformation as much as possible in order to reduce the computational cost and the number of final-solving variables and residual equations constructed in NEM_JFNK. Then efficient physics-based preconditioners are successfully developed by approximating the matrices of the diffusion and removal terms, transverse leakage terms using the three-adjacent-node quadratic fitting methods, and scatter source terms, which make full use of the traditional power iteration. In addition, the Eisenstat-Walker forcing terms are used in the developed NEM_JFNK method to adaptively choose the convergence criterion of linear Krylov iteration within each Newton iteration based on the Newton residuals and to improve computational efficiency further. Finally, the NEM_JFNK code is developed for 3D and multigroup k-eigenvalue problems in neutron diffusion calculations and the detailed study of convergence, computational cost, and efficiency is carried out for several 3D problems. Numerical results show that the developed NEM_JFNK methods have faster convergence speed and are more efficient than the traditional NEM using power iteration, and the speedup ratio is greater for the higher convergence criterion.

Coupling methods for HTR-PM primary circuit

The purpose of this paper is to analyze coupling methods for High-Temperature Gas-Cooled Reactor-Pebble bed Module (HTR-PM) primary circuit, including reactor core, steam generator, and other components. Each component has individual characteristics and is simulated by an independent code package. A feasible way to integrate these independent codes is by using Picard iteration. In practice, the primary circuit is divided into multiple parts based on physical fields or spatial domains. Another coupling method is Jacobian-free Newton-Krylov (JFNK), which is frequently applied for the simulation of coupled multi-physics systems due to its efficiency and capability for large-scale nonlinear problems. To demonstrate the efficiency of coupling method for tightly coupled primary circuit systems, several Picard and JFNK methods were compared consistently in this paper, where all the models, equations and codes for each method are the same. Especially for neutronics, we use the same additional equation for JFNK as the outer iteration step in power iteration. Results show that the equations in the reactor core and outside the reactor core should be solved simultaneously in Picard iteration, otherwise additional thermal-hydraulics iterations are necessary. JFNK methods have shown better convergence behavior than Picard iteration, especially for computational efficiency. Since the nonlinear residual is evaluated directly from nonlinear equation set without iteration, linear preconditioned JFNK has shown efficient convergence, which is recommended by the authors as the first choice for HTR-PM primary circuit if the linear preconditioned JFNK is applicable and the well-performing preconditioner is available. However, for utilizing physical iteration processes, nonlinear preconditioned JFNK is a practical choice when integrating legacy codes.

A Study of Preconditioned Jacobian-Free Newton-Krylov Discontinuous Galerkin Method for Compressible Flows on 3D Hexahedral Grids

Storage requirement and computational efficiency have always been challenges for the efficient implementation of discontinuous Galerkin (DG)methods for real life applications. In this paper, a fully implicit Jacobian-Free Newton-Krylov (JFNK) method is developed in the context of DG discretizations for the three-dimensional compressible Euler and Navier-Stokes equations. Compared with the Jacobian-based methods, the Jacobian-Free approach saves the storage for the Jacobian matrix which can be of great importance for DG methods. Three types of preconditioners are investigated in which the block diagonal preconditioner requires the least storage, while the block LU-SGS and ILU0 preconditioners require more storage but are more computationally efficient. An implicit time-stepping strategy is adopted for the stability of the current solver,which is based upon a hexahedral spatialmesh and the nonlinear solver package Kinsol is used to improve the computational efficiency and robustness. Numerical results demonstrate that the preconditioned JFNK-DG solver can substantially reduce the storage requirement compared with the Jacobian based method without significantly compromising accuracy or efficiency. Furthermore, as a good compromise between efficiency and storage requirement, the ILU0 preconditioner shows the best choice of the preconditioners presented.

Efficient simultaneous solution of multi-physics multi-scale nonlinear coupled system in HTR reactor based on nonlinear elimination method

Pebble bed HTR reactor core is a globally coupled system involving neutronics and thermal-hydraulic field. At the same time, the locally coupled relations exist in fuel sphere temperature, local neutron flux and pebble bed temperature. What’s more, several types of fuel spheres with different burnup and fission power must be treated separately. Due to the essential difference between local and global coupling properties, the method, named nonlinear elimination (NlEm), is developed in this paper to solve the multi-physics/multi-scale system of HTR in a tightly coupled, two-level nonlinear form. In the NlEm framework, the local coupled variables and the global coupled variables are treated separately at two different levels to enhance the computational efficiency. For comparison, the Jacobian-free Newton Krylov (JFNK) method and the Picard method are utilized to solve the multi-physics/multi-scale system respectively. The numerical results reveal that the computational efficiency of NlEm is 6.7 times higher than that of Picard method and 3.5 times higher than that of JFNK method for five transient cases of a two-dimension simplified reactor model. And with the increase of the computational scale and time step size, the speedup ratio of NlEm is further raised.

Formulation of a fully implicit numerical scheme for simulation of two-phase flow in a vertical channel using the Drift-Flux Model

This paper introduces a new fully implicit numerical scheme to simulate two-phase flow in a vertical channel. Two-phase flow is widely associated with the nuclear systems and therefore accurate prediction of the void fraction as a key thermo-hydraulic parameter is crucial in operational and safety analysis. The Drift-Flux Model (DFM) is applicable to a wide range of thermal-hydraulic phenomena in an integral pressurized and boiling water reactors. However, to obtain an accurate and stable solution is a challenge task for successful application of the model. The trend toward the more accurate simulation requires advanced numerical methods. The efficiency of the recently develpoed Jacobian-free Newton-Krylov (JFNK) method as an advanced fully implicit solver have been demonstrated in the published literature. However, there are still unsolved problems with this method. Therefore, a much simpler, more reliable and a fast fully implicit model is developed as a variant of Newton method. The objective of the present study is to investigate the validity, accuracy and CPU usage of the proposed model. Also, good agreement was found between our results and those of the experimental data and the RELAP5 code. Comparison of our results with those of JFNK proved viable for handling thermo hydraulic transients without stability constraints. It is simpler, more accurate and faster than the JFNK method.

Improvement of computational efficiency of circular function-based gas kinetic scheme by using Jacobian-free Newton–Krylov method

The gas-kinetic BGK scheme is a promising method for simulation of inviscid and viscous flows. Different from conventional Navier–Stokes (N-S) solvers, it evaluates the inviscid and viscous fluxes simultaneously by reconstructing the local solution of BGK Boltzmann equation at the cell interface. Due to its inherently superior dissipation property, it usually gives accurate and robust numerical results. However, a notable drawback of the BGK-type scheme is the low computational efficiency. Recently, aimed at reducing the computational effort, a circular function-based BGK (CBGK) scheme was developed. Nevertheless, it is still time consuming because the original scheme used an explicit way in the time integration as in most existing BGK schemes and the convergence speed can be slow. To improve the computational efficiency, an implicit CBGK scheme is developed in this paper by incorporating the Jacobian-free Newton–Krylov (JFNK) method into the scheme. Particularly, the generalized minimal residual (GMRES) approach is employed to iteratively solve the large linear equation system. With the help of the Jacobian-free approach, a faster convergence speed can be achieved without explicitly computing and storing the flux Jacobian, which is usually a large and sparse matrix. In order to reduce the number of GMRES iterations, the preconditioning is also adopted and the Lower-Upper symmetric Gauss-Seidel (LUSGS) scheme is employed as a preconditioner. For validation of the present CBGK-JFNK method, several two-dimensional inviscid and viscous test cases are investigated. The numerical results show that the needed computational time is significantly reduced as compared with the original explicit CBGK scheme and the CBGK-LUSGS method.

A Newton solution for the Superhomogenization method: The PJFNK-SPH

This work presents two novel topics regarding the Superhomogenization method: 1) the formalism for the implementation of the method with the linear Boltzmann Transport Equation, and 2) a Newton algorithm for the solution of the nonlinear problem that arises from the method. These new ideas have been implemented in a continuous finite element discretization in the MAMMOTH reactor physics application. The traditional solution strategy for this nonlinear problem uses a Picard, fixed-point iterative process whereas the new implementation relies on MOOSE’s Preconditioned Jacobian-Free Newton Krylov method to allow for a direct solution. The PJFNK-SPH can converge problems that were either intractable or very difficult to converge with the traditional iterative approach, including geometries with reflectors and vacuum boundary conditions. This is partly due to the underlying Scalable Nonlinear Equations Solvers in PETSc, which are integral to MOOSE and offer Newton damping, line search and trust region methods. The PJFNK-SPH has been implemented and tested for various discretizations of the transport equation included in the Rattlesnake transport solver. Speedups of five times for diffusion and ten to fifteen times for transport were obtained when compared to the traditional Picard approach. The three test problems cover a wide range of applications including a standard Pressurized Water Reactor lattice with control rods, a Transient Reactor Test facility control rod supercell and a prototype fast-thermal reactor. The reference solutions and initial cross sections were obtained from the Serpent 2 Monte Carlo code. The SPH-corrected cross sections yield eigenvalues that are near exact, relative to reference solutions, for reflected geometries, even with reflector regions. In geometries with vacuum boundary conditions the accuracy is problem dependent and solutions can be within a few to a few hundred pcm. The root-mean-square error in the power distribution is below 0.8% of the reference Monte Carlo. There is little benefit from SPH-corrected transport in typical scoping calculations, but for more detailed analyses it can yield superior convergence of the solution in some of the test problems. This PJFNK-SPH approach is currently being used in the modeling of the Transient Test Reactor at Idaho National Laboratory, where full reactor core SPH-corrected cross sections are employed to reduce the homogenization errors in transient or multi-physics calculations. This base implementation of the PJFNK-SPH provides an extremely robust solver and a springboard to further improve the Superhomogenization method in order to better preserve neutron currents, one of the primary deficiencies of the method.

A fully-implicit Particle-In-Cell Monte Carlo Collision code for the simulation of inductively coupled plasmas

We present a fully-implicit electromagnetic Particle-In-Cell Monte Carlo collision code, called NINJA, written for the simulation of inductively coupled plasmas. NINJA employs a kinetic enslaved Jacobian-Free Newton Krylov method to solve self-consistently the interaction between the electromagnetic field generated by the radio-frequency coil and the plasma response. The simulated plasma includes a kinetic description of charged and neutral species as well as the collision processes between them. The algorithm allows simulations with cell sizes much larger than the Debye length and time steps in excess of the Courant–Friedrichs–Lewy condition whilst preserving the conservation of the total energy. The code is applied to the simulation of the plasma discharge of the Linac4 H− ion source at CERN. Simulation results of plasma density, temperature and EEDF are discussed and compared with optical emission spectroscopy measurements. A systematic study of the energy conservation as a function of the numerical parameters is presented.

Fractional-step θ-method for solving singularly perturbed problem in ecology

We consider a predator-prey model arising in ecology that describes a slow-fast dynamical system. The dynamics of the model is expressed by a system of nonlinear differential equations having different time scales. Designing numerical methods for solving problems exhibiting multiple time scales within a system, such as those considered in this paper, has always been a challenging task. To solve such complicated systems, we therefore use an efficient time-stepping algorithm based on fractional-step methods. To develop our algorithm, we first decouple the original system into fast and slow sub-systems, and then apply suitable sub-algorithms based on a class of θ-methods, to discretize each sub-system independently using different time-steps. Then the algorithm for the full problem is obtained by utilizing a higher-order product method by merging the sub-algorithms at each time-step. The nonlinear system resulting from the use of implicit schemes is solved by two different nonlinear solvers, namely, the Jacobian-free Newton-Krylov method and the well-known Anderson’s acceleration technique. The fractional-step θ-methods give us flexibility to use a variety of methods for each sub-system and they are able to preserve qualitative properties of the solution. We analyze these methods for stability and convergence. Several numerical results indicating the efficiency of the proposed method are presented. We also provide numerical results that confirm our theoretical investigations.

A new approach to the stabilization and convergence acceleration in coupled Monte Carlo–CFD calculations: The Newton method via Monte Carlo perturbation theory

Abstract This paper proposes the adoption of Monte Carlo perturbation theory to approximate the Jacobian matrix of coupled neutronics/thermal-hydraulics problems. The projected Jacobian is obtained from the eigenvalue decomposition of the fission matrix, and it is adopted to solve the coupled problem via the Newton method. This avoids numerical differentiations commonly adopted in Jacobian-free Newton–Krylov methods that tend to become expensive and inaccurate in the presence of Monte Carlo statistical errors in the residual. The proposed approach is presented and preliminarily demonstrated for a simple two-dimensional pressurized water reactor case study.

A Newton-based Jacobian-free approach for neutronic-Monte Carlo/thermal-hydraulic static coupled analysis

In the field of nuclear reactor analysis, multi-physics calculations that account for the bonded nature of the neutronic and thermal-hydraulic phenomena are of major importance for both reactor safety and design. So far in the context of Monte-Carlo neutronic analysis a kind of “serial” algorithm has been mainly used for coupling with thermal-hydraulics. The main motivation of this work is the interest for an algorithm that could maintain the distinct treatment of the involved fields within a tight coupling context that could be translated into higher convergence rates and more stable behaviour. This work investigates the possibility of replacing the usually used “serial” iteration with an approximate Newton algorithm. The selected algorithm, called Approximate Block Newton, is actually a version of the Jacobian-free Newton Krylov method suitably modified for coupling mono-disciplinary solvers. Within this Newton scheme the linearised system is solved with a Krylov solver in order to avoid the creation of the Jacobian matrix. A coupling algorithm between Monte-Carlo neutronics and thermal-hydraulics based on the above-mentioned methodology is developed and its performance is analysed. More specifically, OpenMC, a Monte-Carlo neutronics code and COBRA-EN, a thermal-hydraulics code for sub-channel and core analysis, are merged in a coupling scheme using the Approximate Block Newton method aiming to examine the performance of this scheme and compare with that of the “traditional” serial iterative scheme. First results show a clear improvement of the convergence especially in problems where significant coolant density gradients appear.

Anderson acceleration and application to the three-temperature energy equations

The Anderson acceleration method is an algorithm for accelerating the convergence of fixed-point iterations, including the Picard method. Anderson acceleration was first proposed in 1965 and, for some years, has been used successfully to accelerate the convergence of self-consistent field iterations in electronic-structure computations. Recently, the method has attracted growing attention in other application areas and among numerical analysts.Compared with a Newton-like method, an advantage of Anderson acceleration is that there is no need to form the Jacobian matrix. Thus the method is easy to implement. In this paper, an Anderson-accelerated Picard method is employed to solve the three-temperature energy equations, which are a type of strong nonlinear radiation-diffusion equations. Two strategies are used to improve the robustness of the Anderson acceleration method. One strategy is to adjust the iterates when necessary to satisfy the physical constraint. Another strategy is to monitor and, if necessary, reduce the matrix condition number of the least-squares problem in the Anderson-acceleration implementation so that numerical stability can be guaranteed. Numerical results show that the Anderson-accelerated Picard method can solve the three-temperature energy equations efficiently. Compared with the Picard method without acceleration, Anderson acceleration can reduce the number of iterations by at least half. A comparison between a Jacobian-free Newton–Krylov method, the Picard method, and the Anderson-accelerated Picard method is conducted in this paper.