Jacobian-free Newton-Krylov Solver Research Articles

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.

A hyperbolic Poisson solver for tetrahedral grids

This paper presents a robust and efficient Poisson solver that can produce accurate solutions and gradients (e.g., heat flux) on unstructured tetrahedral grids. The solver is constructed based on the hyperbolic method for diffusion, where the Laplacian operator is discretized in the form of a hyperbolic system with solution gradients introduced as additional unknown variables. A practical formula for defining a reference length is proposed, which is needed to properly scale the relaxation length associated with the hyperbolic formulation for scale-invariant computations. The resulting system of residual equations is efficiently solved by a Jacobian-Free Newton-Krylov solver with an implicit defect-correction solver used as a variable-preconditioner. Robustness and superior gradient accuracy are demonstrated for linear and nonlinear Poisson equations through a series of numerical tests with unstructured tetrahedral grids.

Effects of high-frequency damping on iterative convergence of implicit viscous solver

This paper discusses effects of high-frequency damping on iterative convergence of an implicit defect-correction solver for viscous problems. The study targets a finite-volume discretization with a one parameter family of damped viscous schemes. The parameter α controls high-frequency damping: zero damping with α=0, and larger damping for larger α(>0). Convergence rates are predicted for a model diffusion equation by a Fourier analysis over a practical range of α. It is shown that the convergence rate attains its minimum at α=1 on regular quadrilateral grids, and deteriorates for larger values of α. A similar behavior is observed for regular triangular grids. In both quadrilateral and triangular grids, the solver is predicted to diverge for α smaller than approximately 0.5. Numerical results are shown for the diffusion equation and the Navier–Stokes equations on regular and irregular grids. The study suggests that α=1 and 4/3 are suitable values for robust and efficient computations, and α=4/3 is recommended for the diffusion equation, which achieves higher-order accuracy on regular quadrilateral grids. Finally, a Jacobian-Free Newton–Krylov solver with the implicit solver (a low-order Jacobian approximately inverted by a multi-color Gauss–Seidel relaxation scheme) used as a variable preconditioner is recommended for practical computations, which provides robust and efficient convergence for a wide range of α.

A comparison of viscous-plastic sea ice solvers with and without replacement pressure

Recent developments of the explicit elastic-viscous-plastic (EVP) solvers call for a new comparison with implicit solvers for the equations of viscous-plastic sea ice dynamics. In Arctic sea ice simulations, the modified and the adaptive EVP solvers, and the implicit Jacobian–free Newton–Krylov (JFNK) solver are compared against each other. The adaptive EVP method shows convergence rates that are generally similar or even better than those of the modified EVP method, but the convergence of the EVP methods is found to depend dramatically on the use of the replacement pressure (RP). Apparently, using the RP can affect the pseudo-elastic waves in the EVP methods by introducing extra non-physical oscillations so that, in the extreme case, convergence to the VP solution can be lost altogether. The JFNK solver also suffers from higher failure rates with RP implying that with RP the momentum equations are stiffer and more difficult to solve. For practical purposes, both EVP methods can be used efficiently with an unexpectedly low number of sub-cycling steps without compromising the solutions. The differences between the RP solutions and the NoRP solutions (when the RP is not being used) can be reduced with lower thresholds of viscous regularization at the cost of increasing stiffness of the equations, and hence the computational costs of solving them.

A second-order accurate in time IMplicit–EXplicit (IMEX) integration scheme for sea ice dynamics

Current sea ice models use numerical schemes based on a splitting in time between the momentum and continuity equations. Because the ice strength is explicit when solving the momentum equation, this can create unrealistic ice stress gradients when using a large time step. As a consequence, noise develops in the numerical solution and these models can even become numerically unstable at high resolution. To resolve this issue, we have implemented an iterated IMplicit–EXplicit (IMEX) time integration method. This IMEX method was developed in the framework of an already implemented Jacobian-free Newton–Krylov solver. The basic idea of this IMEX approach is to move the explicit calculation of the sea ice thickness and concentration inside the Newton loop such that these tracers evolve during the implicit integration. To obtain second-order accuracy in time, we have also modified the explicit time integration to a second-order Runge–Kutta approach and by introducing a second-order backward difference method for the implicit integration of the momentum equation. These modifications to the code are minor and straightforward. By comparing results with a reference solution obtained with a very small time step, it is shown that the approximate solution is second-order accurate in time. The new method permits to obtain the same accuracy as the splitting in time but by using a time step that is 10 times larger. Results show that the second-order scheme is more than five times more computationally efficient than the splitting in time approach for an equivalent level of error.

A parallel Jacobian-free Newton–Krylov solver for a coupled sea ice-ocean model

The most common representation of sea ice dynamics in climate models assumes that sea ice is a quasi-continuous non-normal fluid with a viscous-plastic rheology. This rheology leads to non-linear sea ice momentum equations that are notoriously difficult to solve. Recently a Jacobian-free Newton–Krylov (JFNK) solver was shown to solve the equations accurately at moderate costs. This solver is extended for massive parallel architectures and vector computers and implemented in a coupled sea ice–ocean general circulation model for climate studies. Numerical performance is discussed along with numerical difficulties in realistic applications with up to 1920 CPUs. The parallel JFNK-solverʼs scalability competes with traditional solvers although the collective communication overhead starts to show a little earlier. When accuracy of the solution is required (i.e. reduction of the residual norm of the momentum equations of more that one or two orders of magnitude) the JFNK-solver is unrivalled in efficiency. The new implementation opens up the opportunity to explore physical mechanisms in the context of large scale sea ice models and climate models and to clearly differentiate these physical effects from numerical artifacts.

A Jacobian-Free Newton Krylov Implicit-Explicit Time Integration Method for Incompressible Flow Problems

AbstractWe have introduced a fully second order IMplicit/EXplicit (IMEX) time in-tegration technique for solving the compressible Euler equations plus nonlinear heat conduction problems (also known as the radiation hydrodynamics problems) in Kadioglu et al., J. Comp. Physics [22,24]. In this paper, we study the implications when this method is applied to the incompressible Navier-Stokes (N-S) equations. The IMEX method is applied to the incompressible flow equations in the following manner. The hyperbolic terms of the flow equations are solved explicitly exploiting the well understood explicit schemes. On the other hand, an implicit strategy is employed for the non-hyperbolic terms. The explicit part is embedded in the implicit step in such a way that it is solved as part of the non-linear function evaluation within the framework of the Jacobian-Free Newton Krylov (JFNK) method [8,29,31]. This is done to obtain a self-consistent implementation of the IMEX method that eliminates the potential order reduction in time accuracy due to the specific operator separation. We employ a simple yet quite effective fractional step projection methodology (similar to those in [11,19,21,30]) as our preconditioner inside the JFNK solver. We present results from several test calculations. For each test, we show second order time convergence. Finally, we present a study for the algorithm performance of the JFNK solver with the new projection method based preconditioner.

CFD High-order Accurate Scheme JFNK Solver

CFD finite difference method high-order accurate scheme can provide more exact flow field solution than second order accurate shceme, but its computing stencil is very border, and hard to get Jacobian matrix for the LU-SGS method, leading to the fact that the LU-SGS convergence speed is slow. A Jacobian-free Newton-Krylov (JFNK) method was developed to a seventh-order accurate scheme HDCS-E8T7, and a kind of nonlinear preconditoner was applied based on traditional 7 diagonal matrix, which was solved with LU-SGS method. In cylinder steady flow case, JFNK method was better than oringal LU-SGS method, nearly half time was saved.

A point collocation method for geometrically nonlinear membranes

This paper describes the development of a numerical model for geometrically nonlinear membranes and evaluates its performance for membranes at static equilibrium. The scheme has several features not commonly seen in structural finite element analysis: the point collocation method, group formulation, and a staggered mesh. In the point collocation finite element method, the partial differential equations are solved at each node instead of by integrating over elements. The group formulation simplifies the handling of nonlinearities by interpolating the nonlinear products of variables, as opposed to seeking the product of independently interpolated variables. The domain is discretized with a staggered mesh of linear triangles and associated polygons. Two sequential gradient recovery operations are performed: first the gradients of the linear triangles are calculated and converted to stresses; then, polygon derivative shape functions derived in this paper are used to determine the internal forces from the stress gradients. The resulting system of nonlinear equations is solved with a Jacobian-free Newton–Krylov solver. The code is first verified using the patch test and the method of manufactured solutions. Then the results are validated using experimental data and benchmark code results in the literature for the Hencky problem (a circular membrane with a fixed perimeter and uniform inflation pressure). The observed rates of convergence for both displacement and radial strain were two. For the configurations and grids used in this investigation, the scheme was suitable for accurately predicting sub-hyperelastic deformations.

A comparison of the Jacobian-free Newton–Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: A serial algorithm study

Numerical convergence properties of a recently developed Jacobian-free Newton–Krylov (JFNK) solver are compared to the ones of the widely used EVP model when solving the sea ice momentum equation with a Viscous-Plastic (VP) formulation. To do so, very accurate reference solutions are produced with an independent Picard solver with an advective time step of 10s and a tight nonlinear convergence criterion on 10, 20, 40, and 80-km grids. Approximate solutions with the JFNK and EVP solvers are obtained for advective time steps of 10, 20 and 30min. Because of an artificial elastic term, the EVP model permits an explicit time-stepping scheme with a relatively large subcycling time step. The elastic waves excited during the subcycling are intended to damp out and almost entirely disappear such that the approximate solution should be close to the VP solution. Results show that residual elastic waves cause the EVP approximate solution to have notable differences with the reference solution and that these differences get more important as the grid is refined. Compared to the reference solution, additional shear lines and zones of strong convergence/divergence are seen in the EVP approximate solution. The approximate solution obtained with the JFNK solver is very close to the reference solution for all spatial resolutions tested.

An efficient mixed-precision, hybrid CPU–GPU implementation of a nonlinearly implicit one-dimensional particle-in-cell algorithm

Recently, an implicit, nonlinearly consistent, energy- and charge-conserving one-dimensional (1D) particle-in-cell method has been proposed for multi-scale, full-f kinetic simulations [G. Chen et al., J. Comput. Phys. 230 (18) (2011)]. The method employs a Jacobian-free Newton–Krylov (JFNK) solver, capable of using very large timesteps without loss of numerical stability or accuracy. A fundamental feature of the method is the segregation of particle-orbit computations from the field solver, while remaining fully self-consistent. This paper describes a very efficient, mixed-precision hybrid CPU–GPU implementation of the 1D implicit PIC algorithm exploiting this feature. The JFNK solver is kept on the CPU in double precision (DP), while the implicit, charge-conserving, and adaptive particle mover is implemented on a GPU (graphics processing unit) using CUDA in single-precision (SP). Performance-oriented optimizations are introduced with the aid of the roofline model. The implicit particle mover algorithm is shown to achieve up to 400GOp/s on a Nvidia GeForce GTX580. This corresponds to 25% absolute GPU efficiency against the peak theoretical performance, and is about 100 times faster than an equivalent single-core CPU (Intel Xeon X5460) compiler-optimized execution. For the test case chosen, the mixed-precision hybrid CPU–GPU solver is shown to over-perform the DP CPU-only serial version by a factor of ∼100, without apparent loss of robustness or accuracy in a challenging long-timescale ion acoustic wave simulation.

The energy conserving particle-in-cell method

A new particle-in-cell (PIC) method, that conserves energy exactly, is presented. The particle equations of motion and the Maxwell’s equations are differenced implicitly in time by the midpoint rule and solved concurrently by a Jacobian-free Newton Krylov (JFNK) solver. Several tests show that the finite grid instability is eliminated in energy conserving PIC simulations, and the method correctly describes the two-stream and Weibel instabilities, conserving exactly the total energy. The computational time of the energy conserving PIC method increases linearly with the number of particles, and it is rather insensitive to the number of grid points and time step. The kinetic enslavement technique can be effectively used to reduce the problem matrix size and the number of JFNK solver iterations.

An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm

This paper discusses a novel fully implicit formulation for a one-dimensional electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov–Poisson formulation), ours is based on a nonlinearly converged Vlasov–Ampére (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (Courant–Friedrichs–Lewy) and spatial (aliasing) instabilities. It is charge- and energy-conserving to numerical round-off for arbitrary implicit time steps (unlike the earlier “energy-conserving” explicit PIC formulation, which only conserves energy in the limit of arbitrarily small time steps). While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling (a deleterious effect for long-term accuracy). The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton–Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.

Implementation of the Jacobian-free Newton–Krylov method for solving the first-order ice sheet momentum balance

We have implemented the Jacobian-free Newton–Krylov (JFNK) method for solving the first-order ice sheet momentum equation in order to improve the numerical performance of the Glimmer-Community Ice Sheet Model (Glimmer-CISM), the land ice component of the Community Earth System Model (CESM). Our JFNK implementation is based on significant re-use of existing code. For example, our physics-based preconditioner uses the original Picard linear solver in Glimmer-CISM. For several test cases spanning a range of geometries and boundary conditions, our JFNK implementation is 1.8–3.6 times more efficient than the standard Picard solver in Glimmer-CISM. Importantly, this computational gain of JFNK over the Picard solver increases when refining the grid. Global convergence of the JFNK solver has been significantly improved by rescaling the equation for the basal boundary condition and through the use of an inexact Newton method. While a diverse set of test cases show that our JFNK implementation is usually robust, for some problems it may fail to converge with increasing resolution (as does the Picard solver). Globalization through parameter continuation did not remedy this problem and future work to improve robustness will explore a combination of Picard and JFNK and the use of homotopy methods.

Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton–Krylov method

We have implemented the Jacobian-free Newton–Krylov (JFNK) method to solve the sea ice momentum equation with a viscous-plastic (VP) formulation. The JFNK method has many advantages: the system matrix (the Jacobian) does not need to be formed and stored, the method is parallelizable and the convergence can be nearly quadratic in the vicinity of the solution. The convergence rate of our JFNK implementation is characterized by two phases: an initial phase with slow convergence and a fast phase for which the residual norm decreases significantly from one Newton iteration to the next. Because of this fast phase, the computational gain of the JFNK method over the standard solver used in existing VP models increases with the required drop in the residual norm (termination criterion). The JFNK method is between 3 and 6.6 times faster (depending on the spatial resolution and termination criterion) than the standard solver using a preconditioned generalized minimum residual method. Resolutions tested in this study are 80, 40, 20 and 10 km. For a large required drop in the residual norm, both JFNK and standard solvers sometimes do not converge. The failure rate for both solvers increases as the grid is refined but stays relatively small (less than 2.3% of failures). With increasing spatial resolution, the velocity gradients (sea ice deformations) get more and more important. Nonlinear solvers such as the JFNK method tend to have difficulties when there are such sharp structures in the solution. This lack of robustness of both solvers is however a debatable problem as it mostly occurs for large required drops in the residual norm. Furthermore, when it occurs, it usually affects only a few grid cells, i.e., the residual is small for all the velocity components except in very localized regions. Globalization approaches for the JFNK solver, such as the line search method, have not yet proven to be successful. Further investigation is needed.