Jacobian-free Newton-Krylov Method Research Articles

Multi-Passage harmonic balance method for flutter prediction in multi-row configurations

A multi-passage harmonic balance method is presented for flutter prediction in multi-row configurations, and particular focus is on simulation of more than two blade rows. The unsteady flows in such configurations contain two kinds of multiple perturbations: one refers to those who have different frequencies and the other one refers to those who have the same frequency but different inter-blade phase angles (IBPAs). Most of current harmonic balance methods only consider the former one and the applications have been restricted to two-blade-row configurations. In the few studies that consider more rows, a one-to-one correspondence assumption between frequency and IBPA is usually made putting extra limitation on the used modes. This motivates our work, and a simple but effective idea is proposed by adding passages adopted per row. The theoretical basis of this work is provided by the spinning mode analysis. The present method is developed based on the time-spectral and multi-frequential form of harmonic balance. The treatment of pitchwise and inter-row boundary conditions in the case of multi passages are particularly addressed. A recently-developed Jacobian-free Newton-Krylov method with a simplified LU-SGS preconditioner is employed to solve the resultant harmonic balance equations. For validation, two typical multi-row flutter cases are investigated. The numerical results show the validity and accuracy of the developed method.

On the Use of the Jacobian-Free Newton Krylov Method to Generate One-Group Discontinuity and Super Homogenization Factors for Full-Core Neutron Diffusion Simulations

Coupled Monte Carlo (MC) and thermal-hydraulic analysis is valuable as a design or reference tool but can be slow, especially when implemented in a Picard iteration. Previous work has developed a novel prediction block to achieve convergence with fewer MC simulations. The prediction works in two stages: (1) a surrogate-like model predicts macroscopic cross sections on the fly and (2) a reduced-order neutronic model predicts the flux response to the updated cross sections. The main challenge with the prediction block is that the reduced-order neutronic model cannot reproduce the spatial flux distribution with high fidelity. This paper investigates the well-established Jacobian-free Newton Krylov (JFNK) method to preserve equivalency between a homogeneous (nodal diffusion) solution and a high-fidelity transport (MC) solution. Instead of performing multiple computationally consuming MC simulations, the nonlinear iterative approach iterates on correction parameters, e.g., assembly discontinuity factors (ADFs) or super homogenization (SPH) factors, using unexpensive nodal solutions. The JFNK approach does not require additional overhead from the MC solver to generate flux tallies. Further, the approach iterates on diffusion solutions produced directly from a desired code, thus ensuring that the parameters are compatible with that code. The approach is applied to correct a nodal diffusion solution for a realistic three-dimensional pressurized water reactor core. The results obtained in the paper show that the method is very successful in reproducing the heterogeneous solution (up to 2.5% difference in assembly flux for SPH and 0.3% for ADFs) without needing to modify the source code of the nodal diffusion solver. In addition, the results show that ADFs yield the best agreement and are also stable (i.e., weakly varying) when thermal-hydraulic fields are perturbed.

An Innovational Jacobian-Split Newton-Krylov Method Combining the Advantages of the Jacobian-Free Newton-Krylov Method and the Finite Difference Jacobian-Based Newton-Krylov Method

The Jacobian-free Newton-Krylov (JFNK) method is a widely used and flexible numerical method for solving the neutronic/thermal-hydraulic coupling system. The main property of JFNK is that the Jacobian-vector product is evaluated approximately by finite difference, avoiding the forming and storage of Jacobian explicitly. However, the lack of an efficient preconditioner is a major bottleneck for the JFNK method, leading to poor convergence. The finite difference Jacobian-based Newton-Krylov (DJNK) method is another advanced numerical method, in which the Jacobian matrix is formed and stored explicitly. The DJNK method can provide a better preconditioner for Krylov iteration than JFNK. However, how to compute the Jacobian matrix efficiently and automatically is a key issue for the DJNK method. By fully utilizing the sparsity of the Jacobian matrix and graph coloring algorithm, the Jacobian can be computed efficiently. Unfortunately, when there are dense rows/blocks, a huge computational burden will emerge due to the lack of sparsity, resulting in the extremely poor efficiency of Jacobian computation. In this work, a Jacobian-split Newton-Krylov (JSNK) method is proposed to resolve the dense row/block problem by combining the advantages of JFNK and DJNK. The main feature of the JSNK method is to split the Jacobian matrix into sparse and dense parts. The sparse part of the Jacobian matrix is explicitly constructed using the graph coloring algorithm while for the dense part, the Jacobian-vector product is approximated by finite difference. The computational complexity of the JSNK method is analyzed and compared to the JFNK method and the DJNK method from theoretical and experimental aspects and under different meshes. A simplified two-dimensional (2-D) high-temperature gas-cooled reactor (HTR) model and a simplified 2-D pressurized water reactor model are utilized to demonstrate the superiority of the JSNK method. The numerical results show that the JSNK method successfully resolved the dense rows/blocks. More importantly, its efficiency significantly outperforms the JFNK method and the DJNK method.

A GPU-Accelerated Method for 3D Nonlinear Kelvin Ship Wake Patterns Simulation

The study of ship waves is important for ship detection, coastal erosion and wave drag. This paper proposed a highly paralleled numerical computation method for efficiently simulating three-dimensional nonlinear kelvin waves. First, a numerical model for nonlinear ship waves is established based on potential flow theory, the Jacobian-free Newton–Krylov (JFNK) method and the boundary integral method. To reduce the amount of data stored in the JFNK method and improve the computational efficiency, a banded preconditioner method is then developed by formulating the optimal bandwidth selection rule. After that, a Graphics Process Unit (GPU)-based parallel computing framework is designed, and we used the Compute Unified Device Architecture (CUDA) language to develop a GPU solution. Finally, numerical simulations of 3D nonlinear ship waves under multiple scales are performed by using the GPU and CPU solvers. Simulation results show that the proposed GPU solver is more efficient than the CPU solver with the same accuracy. More than 66% GPU memory can be saved, and the computational speed can be accelerated up to 20 times. Hence, the computation time for Kelvin ship waves simulation can be significantly reduced by applying the GPU parallel numerical scheme, which lays a solid foundation for practical ocean engineering.

Development of a Fully-Implicit Second-Order system Thermal-Hydraulics analysis code for reactor modeling and simulations

There has been continuous advancement in numerical techniques and software engineering in the field of system thermal-hydraulics over the past decade. The current state-of-the-art involves the use of fully-implicit high-resolution schemes in conjunction with the Jacobian-Free Newton-Krylov (JFNK) method for both single-phase and two-phase flows. However, the application of these advanced schemes in system thermal-hydraulics codes is not a common practice. Given the opportunity and the need to develop a new modern system thermal–hydraulic code from scratch, there is no doubt that these recent advancements should be leveraged. In this work, we introduce a new modern object-oriented system thermal-hydraulics code named RETA. RETA is characterized by several key features, including: second-order temporal and spatial discretization schemes for both incompressible and compressible flows, as well as single-phase and two-phase flows; fully-implicit solution schemes utilizing both the Newton and PJFNK methods; and an object-oriented design of major fluid components and physical models, which allows for easy extension. A series of verification and demonstration studies are conducted to verify the implementation, demonstrate the accuracy improvement, and evaluate the performance of the new code. It is concluded that results from the second-order scheme are generally more reliable. Without the concern of a code stability issue, the second-order scheme should be the preferred and default choice for a code user or an analyst.

An implicit particle code with exact energy and charge conservation for electromagnetic studies of dense plasmas

A collisional particle code based on implicit energy- and charge-conserving methods is presented. A modified version of the particle-suppressed Jacobian-Free Newton-Krylov method that can enhance the solver efficiency is introduced. Mathematically, it is shown that this new approach can be viewed as a fixed-point iteration method for the particle positions. The model can exactly conserve global energy and local charge and can efficiently use time steps larger than the plasma period. The algorithm's ability to simulate dense plasmas accurately and efficiently is quantified by simulating the dynamic compression of a plasma slab via a magnetic piston in 1D planar geometry.

A Jacobian-Free Newton–Krylov Method to Solve Tumor Growth Problems with Effective Preconditioning Strategies

A Jacobian-free Newton–Krylov (JFNK) method with effective preconditioning strategies is introduced to solve a diffusion-based tumor growth model, also known as the Fisher–Kolmogorov partial differential equation (PDE). The time discretization of the PDE is based on the backward Euler and the Crank–Nicolson methods. Second-order centered finite differencing is used for the spatial derivatives. We introduce two physics-based preconditioners associated with the first- and second-order temporal discretizations. The theoretical time and spatial accuracies of the numerical scheme are verified through convergence tables and graphs that correspond to different computational settings. We present efficiency studies with and without using the preconditioners. Our numerical findings indicate the excellent performance of the newly proposed preconditioning strategies. In other words, when we turn the preconditioners on, the average number of GMRES and the Newton iterations are significantly reduced.

A Modified JFNK for Solving the HTR Steady State Secondary Circuit Problem

A nuclear power plant is a complex coupling system, which features multi-physics coupling between reactor physics and thermal-hydraulics in the reactor core, as well as the multi-circuit coupling between the primary circuit and the secondary circuit by the shared steam generator (SG). Especially in the pebble-bed modular HTR nuclear power plant, different nuclear steam supply modules are further coupled together through the shared main steam pipes and the related equipment in the secondary circuit, since the special configuration of multiple reactor modules connects to a steam turbine. The JFNK (Jacobian-Free Newton–Krylov) method provides a promising coupling framework to solve the whole HTR nuclear power plant problem, due to its excellent convergence rate and strong robustness. In this work, the JFNK method was modified and applied to the steady-state calculation of the HTR secondary circuit, which plays an important role in simultaneous solutions for the whole HTR nuclear power plant. The main components in the secondary circuit included SG, steam turbine, condenser, feed pump, high/low-pressure heat exchanger, deaerator, as well as the extraction steam from the steam turbine. The results showed that the JFNK method can effectively solve the steady state issue of the HTR secondary circuit. Moreover, the JFNK method could converge well within a wide range of initial values, indicating its strong robustness.

Jacobian-free Newton–Krylov method for the simulation of non-thermal plasma discharges with high-order time integration and physics-based preconditioning

A preconditioning framework for the numerical simulation of non-thermal streamer discharges is developed using the Jacobian-free Newton-Krylov (JFNK) method. A reduced plasma fluid model is considered, consisting of electrons, one positive ion, one negative ion, and the electrostatic potential. The plasma kinetics model includes ionization, electron-ion recombination, electron attachment, electron detachment, and ion-ion recombination. The governing equations are made dimensionless, discretized in space with finite differences, and integrated in time with a fully implicit method based on high-order backward differentiation formulas. The preconditioning framework is based on a linearized form of the governing equations and physics-based operator splitting. The efficiency of the preconditioning strategy is assessed through two test cases: streamer propagation between parallel plates and an axisymmetric pin-to-pin discharge. The fully implicit approach overcomes traditional restrictions in the time step size due to processes such as electron drift, electron diffusion, and dielectric relaxation. Excellent performance is observed through relevant statistics of the JFNK solver, although the number of linear iterations increases for the pin-to-pin discharge when nonlinear numerical boundary conditions are imposed at the electrodes. Performance studies show scalability with O(100-1000) processors for O(10M) unknowns with ample room for optimization.

Operator split, Picard iteration and JFNK methods based on nonlinear CMFD for transient full core models in the coupling multiphysics environment

Operator Split (OS), Picard iteration and Jacobian-free Newton Krylov (JFNK) methods based on nonlinear coarse mesh finite difference (CMFD) with nodal expansion methods (NEM) are respectively introduced to solve the transient full core coupled models in the COupling Multiphysics Environment (COME). In the COME, in addition to the traditional OS and Picard coupled methods, we develop two different JFNK-based coupled methods with physics and algebraic-based hybrid preconditioners at each time step: OS_JFNK and JFNK coupled methods. Specifically, the efficient JFNK method is applied to solve CMFD discrete models of only neutronics in the OS framework for OS_JFNK coupled strategies or whole Neutronics-Thermal Hydraulic (N-TH) coupled models for JFNK coupled methods. Then nonlinear corrective coupling coefficients are updated by the two-node NEM every few Newton steps at each time step. Six cases of NEACRP 3D core transient benchmarks are analyzed and numerical solutions of OS, OS_JFNK, Picard and JFNK methods in the COME generally agree well with other published solutions. JFNK and Picard coupled methods can track the reference solutions better than OS and OS_JFNK coupled methods due to multiple iterations or simultaneous solutions between neutronics and TH models in the tight coupled form. Furthermore, OS_JFNK or JFNK methods can obtain higher computational efficiency than OS or Picard methods with or without the fission source extrapolation acceleration, respectively, which indicates the potential of the developed COME using different coupled methods for the transient full core multiphysics coupled problems.

A modified JFNK with line search method for solving k-eigenvalue neutronics problems with thermal-hydraulics feedback

The k-eigenvalue neutronics/thermal-hydraulics coupling calculation is a key issue for reactor design and analysis. Jacobian-free Newton-Krylov (JFNK) method, featured with super-linear convergence rate and high efficiency, has been attracting more and more attention to solve the multi-physics coupling problem. However, it may converge to the high-order eigenmode because of the multiple solutions nature of the k-eigenvalue form of multi-physics coupling issue. Based on our previous work, a modified JFNK with a line search method is proposed in this work, which can find the fundamental eigenmode together with thermal-hydraulics feedback in a wide range of initial values. In detail, the existing modified JFNK method is combined with the line search strategy, so that the intermediate iterative solution can avoid a sudden divergence and be adjusted into a convergence basin smoothly. Two simplified 2-D homogeneous reactor models, a PWR model, and an HTR model, are utilized to evaluate the performance of the newly proposed JFNK method. The results show that the performance of this proposed JFNK is more robust than the existing JFNK-based methods.

Jacobian-free Newton Krylov coarse mesh finite difference algorithm based on high-order nodal expansion method for three-dimensional nuclear reactor pin-by-pin multiphysics coupled models

The nuclear reactor core is a large-scale, complicated and tight-coupled nonlinear multiphysics coupled system including neutron transport, heat conduction of nuclear fuel pins, coolant flow, heat transfer and so on. To meet the demands of the advanced nuclear reactor analysis and design, there is an urgent need to reduce the computational costs and improve the convergence rate for the three-dimensional (3D) nuclear reactor core multiphysics coupled models, especially for a high-fidelity pin-by-pin (at the level of fuel pin by fuel pin) full core simulation. In this work, we develop CNJ, an efficient Jacobian free Newton Krylov (JFNK) coarse mesh finite difference (CMFD) algorithm based on the nodal expansion method (NEM) to simultaneously solve the complicated pin-by-pin neutronics-thermal hydraulic (N-TH) coupled models. The high-order NEM can improve the numerical accuracy of the CMFD method on a coarse mesh and the JFNK method can ensure the rapid convergence and high efficiency on large-scale and nonlinear coupled discrete systems. By combining the CMFD, NEM and JFNK methods and making full use of their respective advantages, the CNJ algorithm can obtain high accuracy and efficiency even on a coarse mesh to greatly reduce the number of solution variables and the computational cost of the 3D coupled models. Finally, two representative and complicated 3D pin-by-pin N-TH coupled models are analyzed to evaluate the numerical accuracy and efficiency advantage of CNJ in detailed.

On the Jacobian approximation in sea ice models with viscous-plastic rheology

Viscous-plastic rheology is widely used in the sea ice modeling community at increasingly high resolutions. Due to the high degree of non-linearity of the rheological constraints, accurate approximation of the Jacobian is required to improve the efficiency of the implicit solvers of the sea ice momentum equation in pack ice. We consider the analytical Jacobian of the ice momentum equation and assess its approximation errors in the Jacobian-free Newton–Krylov (JFNK) method and in the family of more traditional schemes which neglect the dependence of viscosity coefficients on the deformation tensor. It is shown that this dependence provides a substantial contribution to the Jacobian, especially in the regions enriched by linear kinematic features like ridges and elongated polynyas. Numerical experiments indicate that performance of the Newton solvers is also sensitive to errors associated with inexact computation of the search direction, that may be caused, in particular, by numerical approximations of the Jacobian which violate its dissipative property. Based on this analysis, an improved selective damping strategy for the Newtonian solver of the momentum equation is proposed. A series of numerical experiments conducted in simulated pack ice environment demonstrate faster convergence of the updated solver with analytical Jacobian as compared to the one based on the selectively damped JFNK method with inexact GMRES solver in the inner loop.

An effective method for modeling 1D two-phase two-fluid six-equation model with automatic differentiation approach

The modeling of the two-fluid model with plenty of realistic closure correlations has been a tough problem because of the difficulty of forming the Jacobian matrix manually. To predict the phase appearance/disappearance in two-phase flow problems, the high-order spatial and fully implicit temporal discretization typically should be adopted. However, considering the requirement of a fairly accurate evaluation of the Jacobian matrix, the implementation of high-order schemes is a challenging task. In this study, the automatic differentiation (AD) method was employed to overcome these difficulties, which seems the first demonstration of the high-order scheme and the AD method in the two-fluid model with realistic relations. With the AD method, the high-order fully implicit scheme, and the preconditioned Jacobian-free Newton-Krylov (PJFNK) method, the two-fluid code was developed to solve the two-phase flow and heat transfer problems in this study. Through simulations of a series of hydrodynamics two-phase flow problems, this code was verified and validated in both transients and boiling conditions. Meanwhile, code-to-code comparisons with well-validated code, RELAP5, were presented. The results demonstrate the promising performance of the developed code for solving 1D two-phase flow problems.

A multi-level nonlinear elimination-based JFNK method for multi-scale multi-physics coupling problem in pebble-bed HTR

In the pebble-bed high-temperature gas-cooled reactor, various physical phenomena occur at different spatial scales and are coupled together tightly in a nonlinear manner. More specifically, the local temperature profiles inside the fuel sphere and TRISO particle, which are essential in reactor design and safety analysis, are coupled with the global physical fields, such as the pebble bed temperature, helium temperature, as well as the fission heat source. Therefore, the multi-physics calculation combined with multi-scale model should be employed to describe the reactor behavior. The computational efficiency is the critical issue for the numerical solution of such global and local coupled systems. Jacobian-free Newton Krylov (JFNK) algorithm is an excellent candidate for multi-physics coupling issues due to its high convergence rate and superior computational performance, but it doesn’t consider the physical features at different scales. In this work, a multi-level nonlinear elimination is proposed within the JFNK framework to deal with the multi-scale coupled system in high-temperature gas-cooled reactor. The numerous local multi-scale temperature variables, including the fuel sphere and TRISO particle temperatures, are eliminated efficiently based on the Schur complement factorization due to its well sparse matrix structure. As a result, the number of unknowns is significantly reduced and the computational performance is improved. A steady-state pebble-bed HTR model is used in this work, and the numerical results show that multi-level nonlinear elimination-based JFNK method has high computational efficiency and numerical robustness, whose performance is 4 ∼ 6 times higher than that of the original JFNK algorithm.

Large-Eddy Simulation of Transonic Buffet Using Matrix-Free Discontinuous Galerkin Method

We present an implicit large-eddy simulation of transonic buffet over the OAT15A supercritical airfoil at Mach number 0.73, angle of attack 3.5 deg, and Reynolds number . The simulation is performed using a matrix-free discontinuous Galerkin (DG) method and a diagonally implicit Runge–Kutta scheme on graphics processor units. We propose a Jacobian-free Newton–Krylov method to solve nonlinear systems arising from the discretization of the Navier–Stokes equations. The method successfully predicts the buffet onset, the buffet frequency, and turbulence statistics owing to the high-order DG discretization and an efficient mesh refinement for the laminar and turbulent boundary layers. A number of physical phenomena present in the experiment are captured in our simulation, including periodical low-frequency oscillations of shock wave in the streamwise direction, strong shear layer detached from the shock wave due to shock-wave/boundary-layer interaction and small-scale structures broken down by the shear-layer instability in the transition region, and shock-induced flow separation. The pressure coefficient, the root mean square of the fluctuating pressure, and the streamwise range of the shock wave oscillation agree well with the experimental data. The results suggest that the proposed method can accurately predict the onset of turbulence and buffet phenomena at high Reynolds numbers without a subgrid scale model or a wall model.

An adaptive scalable fully implicit algorithm based on stabilized finite element for reduced visco-resistive MHD

The magnetohydrodynamics (MHD) equations are continuum models used in the study of a wide range of plasma physics systems, including the evolution of complex plasma dynamics in tokamak disruptions. However, efficient numerical solution methods for MHD are extremely challenging due to disparate time and length scales, strong hyperbolic phenomena, and nonlinearity. Therefore the development of scalable, implicit MHD algorithms and high-resolution adaptive mesh refinement strategies is of considerable importance. In this work, we develop a high-order stabilized finite-element algorithm for the reduced visco-resistive MHD equations based on the MFEM finite element library (mfem.org). The scheme is fully implicit, solved with the Jacobian-free Newton-Krylov (JFNK) method with a physics-based preconditioning strategy. Our preconditioning strategy is a generalization of the physics-based preconditioning methods in Chacón et al. (2002) [3] to adaptive, stabilized finite elements. Algebraic multigrid methods are used to invert sub-block operators to achieve scalability. A parallel adaptive mesh refinement scheme with dynamic load-balancing is implemented to efficiently resolve the multi-scale spatial features of the system. Our implementation uses the MFEM framework, which provides arbitrary-order polynomials and flexible adaptive conforming and non-conforming meshes capabilities. Results demonstrate the accuracy, efficiency, and scalability of the implicit scheme in the presence of large scale disparity. The potential of the AMR approach is demonstrated on an island coalescence problem in the high Lundquist-number regime (≥107) with the successful resolution of plasmoid instabilities and thin current sheets.

Jacobian-Free Newton–Krylov Method for the Simulation of Non-Thermal Plasma Discharges with High-Order Time Integration and Physics-Based Preconditioning

Jacobian-Free Newton–Krylov Method for the Simulation of Non-Thermal Plasma Discharges with High-Order Time Integration and Physics-Based Preconditioning

Tusas: A fully implicit parallel approach for coupled phase-field equations

We develop a fully-coupled, fully-implicit approach for phase-field modeling of solidification in metals and alloys. Predictive simulation of solidification in pure metals and metal alloys remains a significant challenge in the field of materials science, as microstructure formation during the solidification process plays a critical role in the properties and performance of the solid material. Our simulation approach consists of a finite element spatial discretization of the fully-coupled nonlinear system of partial differential equations at the microscale, which is treated implicitly in time with a preconditioned Jacobian-free Newton-Krylov method. The approach is algorithmically scalable as well as efficient due to an effective preconditioning strategy based on algebraic multigrid and block factorization. We implement this approach in the open-source Tusas framework, which is a general, flexible tool developed in C++ for solving coupled systems of nonlinear partial differential equations. The performance of our approach is analyzed in terms of algorithmic scalability and efficiency, while the computational performance of Tusas is presented in terms of parallel scalability and efficiency on emerging heterogeneous architectures. We demonstrate that modern algorithms, discretizations, and computational science, and heterogeneous hardware provide a robust route for predictive phase-field simulation of microstructure evolution during additive manufacturing.

Assessment of a symmetry-preserving JFNK method for atmospheric convection

Numerical simulations of nonhydrostatic atmospheric flow, based on linearly decoupled semi-implicit or fully-implicit techniques, usually solve linear systems by a pre-conditioned Krylov method without preserving the skew-symmetry of convective operators. We propose to perform atmospheric simulations in such a fully-implicit manner that the difference operators preserve both the skew-symmetry and the tightly nonlinear coupling of the differential operators. We demonstrate that a symmetry-preserving Jacobian-free Newton-Krylov (JFNK) method mimics a balance between convective transport and turbulence dissipation. We present a wavelet method as an effective symmetry preserving discretization technique. The symmetry-preserving JFNK method for solving equations of nonhydrostatic atmospheric flows has been examined using two benchmark simulations of penetrative convection – a) dry thermals rising in a neutrally stratified and stably stratified environment, and b) urban heat island circulations for effects of the surface heat flux H0 varying in the range of 25≤H0≤930 Wm−2. The results show that an eddy viscosity model provides the necessary dissipation of the subgrid-scale modes, while the symmetry-preserving JFNK method provides the conservation of mass and energy at a satisfactory level. Comparisons of the results from a laboratory experiment of heat island circulation and a field measurement of potential temperature also suggest the modelling accuracy of the present symmetry-preserving JFNK framework.