Newton Krylov Research Articles

An Innovative and Efficient Implementation of Matrix-Free Newton Krylov Method for Neutronics/Thermal-hydraulics Coupling Simulation

The core physical behavior of reactors is essentially the result of multi-physical fields coupling feedback. High-fidelity neutronics/thermal-hydraulics (N/TH) analysis can simulate and predict nuclear reactor core phenomena realistically, providing advanced and reliable technical means during the design and safety analysis of nuclear reactor. In this work, an efficient and robustness coupling method using power density as the coupling parameter, Matrix-Free Newton Krylov (MFNK) method, is successfully developed and innovatively implemented in HNET for high-fidelity N/TH coupling simulation. To enhance the efficiency and stability, the multi-level generalized equivalence theory-based CMFD (ML-gCMFD) iterative acceleration method and ML-gCMFD coupling acceleration method are proposed. In addition, the nonlinear preconditioning and hybrid perturbation size formula are implemented to further improve the convergence. Finally, to evaluate the numerical accuracy, convergence, efficiency and stability of MFNK method, a series of representative problems, including a three-dimensional (3D) single fuel pin problem, VERA Benchmark Problem 6, and VERA Benchmark Problem 7, are analyzed by comparing with the current N/TH coupling methods. Numerical results indicate that MFNK method can obtain strong stability, high convergence performance, and relatively high computational efficiency while ensuring high accuracy. It demonstrates that MFNK method has significant performance advantages and potential for high-fidelity N/TH coupling simulation.

An unconditionally stable scheme for the immersed boundary method with application in cardiac mechanics

The stability of the immersed boundary (IB) method is a challenge in simulating fluid–structure interaction problems, where time step constraints are significantly stricter than in pure fluid simulations. We propose a novel unconditionally stable scheme for the immersed boundary finite element (IBFE) method. The structure is handled implicitly and characterized by strain energy functions, rather than being modeled as fibers or membranes. Through energy estimate, we prove the unconditional stability of the fully discrete approximation in the absence of the convective term. In real simulations of cardiac mechanics problems, the time step is much larger, only limited by the Courant–Friedrichs–Lewy condition of the fluid. The novelty of this work lies in the combination of dual interpolation and distribution operators, the Jacobian-free Newton–Krylov method for solving nonlinear algebraic systems, and the semi-Lagrangian method for handling the convective term. To validate the effectiveness and accuracy of our approach, we present various benchmarks and conduct a quasi-static simulation of a three-dimensional real left ventricular model. We have shown that the numerical stability of our scheme is very robust even with much larger time step compared to conventional explicit IB methods. Our work paves the way for further works on efficient solvers of the IBFE method.

Bifurcation analysis in the regularized four-sided cavity flow: Equilibrium states in a D2-symmetric fluid system

In this study, we investigate a handful of equilibrium states and their hydrodynamic stability for the incompressible flow in a square cavity driven by the tangential simultaneous motion of all of its lids at the same speed. This problem has been investigated in the recent past, albeit only partially and always disregarding the singular boundary conditions at the corners. The implications of the system symmetries have not been rigorously addressed, either. In our study, we introduce a regularized version of the boundary conditions to avoid the corner singularities, in a way that preserves the main features of the original setup. In addition, the Navier–Stokes equations are discretized by means of highly accurate Chebyshev spectral methods that provide exponential convergence of all computed flows and consistently eliminate any potential source of structural instability of the bifurcation scenario. We employ Newton–Krylov solvers, implemented within continuation algorithms, to accurately compute equilibrium solutions. Linear stability analysis of both the primary symmetric base flow and the secondary asymmetric states uncovers new branches of fully asymmetric steady states. The analysis has allowed identification of six previously undetected bifurcations, all of which are associated with the disruption of either the rotational invariance or the reflection symmetry. Some of these bifurcations have been found to be quite clustered in some regions of the parameter space, which points at the underlying action of higher codimension mechanisms. Notably, all the bifurcations reported occur within the range of low to moderate Reynolds numbers, making the regularized four-sided lid-driven cavity flow a reliable benchmark for assessing different numerical schemes in the context of Navier–Stokes equivariant bifurcation theory.

FreeGSNKE: A Python-based dynamic free-boundary toroidal plasma equilibrium solver

We present a Python-based numerical solver for the two-dimensional dynamic plasma equilibrium problem. We model the time evolution of toroidally symmetric free-boundary tokamak plasma equilibria in the presence of the non-linear magnetohydrodynamic coupling with both currents in the “active” poloidal field coils, with assigned applied voltages, and eddy currents in the tokamak passive structures. FreeGSNKE (FreeGS Newton–Krylov Evolutive) builds and expands on the framework provided by the Python package FreeGS (Free boundary Grad–Shafranov). FreeGS solves the static free-boundary Grad–Shafranov (GS) problem, discretized in space using finite differences, by means of Picard iterations. FreeGSNKE introduces: (i) a solver for the static free-boundary GS problem based on the Newton–Krylov (NK) method, with improved stability and convergence properties; (ii) a solver for the linearized dynamic plasma equilibrium problem; and (iii) a solver for the non-linear dynamic problem, based on the NK method. We propose a novel “staggered” solution strategy for the non-linear problem, in which we make use of a set of equivalent formulations of the non-linear dynamic problem we derive. The alternation of NK solution steps in the currents and in the plasma flux lends this strategy an increased resilience to co-linearity and stagnation problems, resulting in favorable convergence properties. FreeGSNKE can be used for any user-defined tokamak geometry and coil configuration. FreeGSNKE's flexibility and ease of use make it a suitably robust control-oriented simulator of plasma magnetic equilibria. FreeGSNKE is entirely written in Python and easily interfaced with Python libraries, which facilitates machine learning based approaches to plasma control.

An Efficient and Robust ILU(k) Preconditioner for Steady-State Neutron Diffusion Problem Based on MOOSE

Jacobian-free Newton Krylov (JFNK) is an attractive method to solve nonlinear equations in the nuclear engineering community, and has been successfully applied to steady-state neutron diffusion k-eigenvalue problems and multi-physics coupling problems. Preconditioning technique plays an important role in the JFNK algorithm, significantly affecting its computational efficiency. The key point is how to automatically construct a high-quality preconditioning matrix that can improve the convergence rate and perform the preconditioning matrix factorization efficiently and robustly. A reordering-based ILU(k) preconditioner is proposed to achieve the above objectives. In detail, the finite difference technique combined with the coloring algorithm is utilized to automatically construct a preconditioning matrix with low computational cost. Furthermore, the reordering algorithm is employed for the ILU(k) to reduce the additional non-zero elements and pursue robust computational performance. A 2D LRA neutron steady-state benchmark problem is used to evaluate the performance of the proposed preconditioning technique, and a steady-state neutron diffusion k-eigenvalue problem with thermal-hydraulic feedback is also utilized as a supplement. The results show that coloring algorithms can automatically and efficiently construct the preconditioning matrix. The computational efficiency of the FDP with coloring could be about 60 times higher than that of the preconditioner without the coloring algorithm. The reordering-based ILU(k) preconditioner shows excellent robustness, avoiding the effect of the fill-in level k choice in incomplete LU factorization. Moreover, its performances under different fill-in levels are comparable to the optimal computational cost with natural ordering.

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.

Domain Decomposition for 3-D Nonlinear Magnetostatic Problems: Newton–Krylov–Schur Versus Schur–Newton–Krylov Methods

Domain Decomposition for 3-D Nonlinear Magnetostatic Problems: Newton–Krylov–Schur Versus Schur–Newton–Krylov Methods

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.

SedTrace 1.0: a Julia-based framework for generating and running reactive-transport models of marine sediment diagenesis specializing in trace elements and isotopes

Abstract. Trace elements and isotopes (TEIs) are important tools in studying ocean biogeochemistry. Understanding their modern ocean budgets and using their sedimentary records to reconstruct paleoceanographic conditions require a mechanistic understanding of the diagenesis of TEIs, yet the lack of appropriate modeling tools has limited our ability to perform such research. Here I introduce SedTrace, a modeling framework that can be used to generate reactive-transport code for modeling marine sediment diagenesis and assist model simulation using advanced numerical tools in Julia. SedTrace enables mechanistic TEI modeling by providing flexible tools for pH and speciation modeling, which are essential in studying TEI diagenesis. SedTrace is designed to solve one particular challenge facing users of diagenetic models: existing models are usually case-specific and not easily adaptable for new problems such that the user has to choose between modifying published code and writing their own code, both of which demand strong coding skills. To lower this barrier, SedTrace can generate diagenetic models only requiring the user to supply Excel spreadsheets containing necessary model information. The resulting code is clearly structured and readable, and it is integrated with Julia's differential equation solving ecosystems, utilizing tools such as automatic differentiation, sparse numerical methods, Newton–Krylov solvers and preconditioners. This allows efficient solution of large systems of stiff diagenetic equations. I demonstrate the capacity of SedTrace using case studies of modeling the diagenesis of pH as well as radiogenic and stable isotopes of TEIs.

A hybrid numerical method for non-linear transient heat conduction problems with temperature-dependent thermal conductivity

This paper introduces a hybrid numerical technique aimed at effectively addressing two-dimensional (2D) non-linear transient heat conduction problems, featuring temperature-dependent thermal conductivity. This scheme involves integrating the Krylov deferred correction (KDC) method in the temporal discretization and the generalized finite difference method (GFDM) in the spatial discretization. The core of the KDC method introduces a new unknown variable in the form of the first-order time derivative of temperature, leading to a non-linear equation after temporal discretization. Subsequently, the corresponding non-linear equation is solved in the spatial domain using the GFDM, supported by the Jacobian-free Newton–Krylov (JFNK) solver and fourth-order Taylor series expansion. The accuracy and stability of the hybrid approach are verified through two numerical experiments, demonstrating the strong performance of the present algorithm known as the KDC-GFDM for the simulations of the interested problems in long-time intervals.

Physics‐based preconditioning of Jacobian‐free Newton–Krylov solver for Navier–Stokes equations using nodal integral method

AbstractThe nodal integral methods (NIMs) have found widespread use in the nuclear industry for neutron transport problems due to their high accuracy. However, despite considerable development, these methods have limited acceptability among the fluid flow community. One major drawback of these methods is the lack of robust and efficient nonlinear solvers for the algebraic equations resulting from discretization. Since its inception, several modifications have been made to make NIMs more agile, efficient, and accurate. Modified nodal integral method (MNIM) and modified MNIM (M2NIM) are the two most recent and efficient versions of the NIM for fluid flow problems. M2NIM modifies the MNIM by replacing the current time convective velocity with the previous time convective velocity, leading to faster convergence albeit with reduced accuracy. This work proposes a preconditioned Jacobian‐free Newton–Krylov approach for solving the Navier–Stokes equation using MNIM. The Krylov solvers do not generally work well without an appropriate preconditioner. Therefore, M2NIM is used here as a preconditioner to accelerate the solution of MNIM. Due to pressure–velocity coupling in the Navier–Stokes equation, developing a quality preconditioner for these equations needs significant effort. The momentum equation is solved using the time‐splitting alternate direction implicit method. The velocities obtained from the solution are then used to solve the pressure Poisson equation. The computational results for the Navier–Stokes equation are presented to underscore the advantages of the developed algorithm.

Transient Simulation and Parameter Sensitivity Analysis of Godiva Experiment Based on MOOSE Platform

The fast-neutron burst reactor is a chain reactor that can operate in a prompt critical state. In order to ensure the operational safety of the fast-neutron-pulse reactor and prevent the supercritical pulse from causing physical damage to the material, it is necessary to simulate and analyze the pulse operating conditions of the fast-neutron-pulse reactor. Godiva-I is a spherical assembly of highly enriched uranium metal made during the 1950s. A prompt-critical transient in such a nuclear system impels a quick power excursion, which will cause a temperature rise and a subsequent reactivity reduction because of the metal sphere’s expansion. The overall transient lasts for a few fractions of a millisecond. Based on the point kinetics and Monte Carlo method, the temporal and spatial characteristics of transient input power were calculated, the difference of the average reactivity temperature coefficient between uniform density and non-uniform density was compared, and the transient power distribution condition was loaded into the thermal–mechanics calculation of the MOOSE platform; thus, the pulse process of Godiva-I with different initial reactivity periods was simulated. The JFNK (Jacobian–Free–Newton–Krylov) direct method and multi-app indirect method were used to analyze the transient response of the pulse dynamic process using the heat conduction module and tenor mechanics module, respectively. After considering the influence of the inertia effect and wall-reflected neutrons, the simulation results were much closer to the experimental values. Based on the stochastic tools module, the uncertainty propagation and sensitivity analysis of the Godiva-I model were carried out, the uncertainty of external surface displacement of Godiva under input disturbance of material properties and heat source amplitude factor was obtained, and the sensitivity of different input parameters to output parameters was quantified. The research results can lay a technical foundation for the thermal–mechanics coupling analysis and uncertainty quantification of the metal fast reactor.

On a fully-implicit VMS-stabilized FE formulation for low Mach number compressible resistive MHD with application to MCF

This study presents the development and evaluation of a fully-implicit variational multiscale (VMS) stabilized unstructured finite element (FE) formulation for compressible magnetohydrodynamics (MHD) model, at low Mach number regime. The model describes the dynamics of a compressible conducting fluid in the low Mach number limit in the presence of electromagnetic fields and can be used to study aspects of astrophysical phenomena, important science and technology applications, and basic plasma physics phenomena. The specific applications that motivate this study are macroscopic simulations of the longer time-scale stability and disruptions of magnetic confinement fusion (MCF) devices, specifically the ITER tokamak. The discussion considers the development of the VMS FE representation, the structure of the stabilizing terms that deal with significant convective flows, the stabilization of the nearly incompressible response of the fluid flow, and the stabilization of the constraint that enforces the solenoidal involution on the magnetic field. The nonlinear discretized system is solved with scalable preconditioned Newton–Krylov iterative methods, which employs a multiphysics block preconditioning method based on approximate block factorizations and Schur complements. The study presents an evaluation of the VMS method on a 2D cartesian tearing mode instability, and illustrates the scalability of the solvers on MCF relevant problems. A set of results are also presented for longer time-scale stability and disruptions for the ITER tokamak. These include a vertical displacement event (VDE), and a (1,1) internal kink mode. The formulation is demonstrated to be scalable and also reasonably robust with respect to the Lundquist number scaling.

Implementation of a phase field damage model with a nonlinear evolution equation in an FFT-based solver

This paper focuses on the numerical implementation of phase-field models of fracture using the Fast Fourier Transform (FFT) based numerical method. Recent studies on phase-field models focus on the discussions of the choice of the value of regularization length proposed to smear the discontinuity of the sharp crack. Some studies argue that it should be considered as a material property because it has a significant impact on the mechanical behavior of a material in some phase-field models, for instance, in the model proposed by Miehe. However, our results in this study for heterogeneous materials have shown that the choice of regularization length not only affects the macroscopic mechanical behavior but also the local crack propagation patterns. As a result, it can be challenging to select an appropriate value that produces both accurate macroscopic responses and local crack patterns for certain phase-field models, such as Miehe’s model. Thus, the phase-field model proposed by Wu, which has been proven to reduce the length sensitivity for homogeneous material, has been successfully implemented in an FFT-based solver with the application of Newton–Krylov algorithm.The length sensitivity of Wu’s phase-field model for heterogeneous materials is also investigated in this paper. Our tests show that Wu’s phase-field model has partial length sensitivity for heterogeneous materials. However, the main reason for this sensitivity is that one phase enters the damage zone of other phase, which is different from Miehe’s model. As a result, a set of criteria for accurately choosing the regularization value has been established for Wu’s model in this work. Meanwhile, it has also been found that Wu’s model may be more suitable than Miehe’s model for brittle failure due to the introduction of an elastic stage.

Dissipation and time step scaling strategies for low and high Mach number flows

Aircraft aerodynamic design often involves evaluating flow conditions that span low subsonic to transonic or even supersonic Mach numbers. Compressible flow solvers are a natural choice for such design problems, but these solvers suffer from reduced accuracy and efficiency at low Mach numbers. In addition, simulations with supersonic conditions can be challenging to converge because of large gradients in the flow field. This paper presents three contributions to address these issues in the context of an approximate Newton–Krylov solver for the Reynolds-averaged Navier–Stokes equations. First, we propose a method for scaling the artificial dissipation terms in the Jameson–Schmidt–Turkel scheme to improve its accuracy at low Mach numbers while retaining the simplicity of the original scalar dissipation formulation. Second, we show that characteristic time-stepping combined with an approximate Newton method can accelerate convergence for low Mach number flows by reducing the stiffness in the linear system for each Newton iteration. Third, we introduce a dissipation-based continuation method for flows with shocks that improves robustness and accelerates convergence without sacrificing accuracy. These methods can make compressible flow solvers more accurate and efficient across low and high Mach number regimes.

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.

FEOTS v0.0.0: a new offline code for the fast equilibration of tracers in the ocean

Abstract. In this paper we introduce a new software framework for the offline calculation of tracer transport in the ocean. The Fast Equilibration of Ocean Tracers Software (FEOTS) is an end-to-end set of tools to efficiently calculate tracer distributions on a global or regional sub-domain using transport operators diagnosed from a comprehensive ocean model. To the best of our knowledge, this is the first application of a transport matrix model to an eddying ocean state. While a Newton–Krylov-based equilibration capability is still under development and not presented here, we demonstrate in this paper the transient modeling capabilities of FEOTS in an application focused on the Argentine Basin, where intense eddy activity and the Zapiola Anticyclone lead to strong mixing of water masses. The demonstration illustrates progress in developing offline passive tracer simulation capabilities, while highlighting the challenges of the impulse response functions approach in capturing tracer transports by a non-linear advection scheme. Our future work will focus on improving the computational efficiency of the code to reduce time-to-solution, using different basis functions to better represent non-linear advection operators, applying FEOTS to a parent model with unstructured grids (Ocean Model for Prediction Across Scales, MPAS-Ocean), and fully implementing a Newton–Krylov steady-state solver.

An improvement of newton – krylov method for solution of nonlinear equations

Solving problems in practice often results in a system of nonlinear equations with a large number of equations and unknowns. Finding the exact solution to this class of equations is very difficult and almost impossible. Recently, with the development of technology, many methods and algorithms have been proposed to approximate the class of these systems of equations. Especially the third-order Newton–Krylov method has solved quite well this class of systems of equations with the third degree of convergence. In this paper, we present a new improvement of the third-order Newton-Krylov method with a quaternary convergence rate and prove the convergence of the iterative formula. In addition, the paper also presents an experimental result to demonstrate the convergence speed of the method

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.