Krylov Subspace Solver Research Articles

A symmetric multigrid-preconditioned Krylov subspace solver for Stokes equations

Numerical solution of discrete PDEs corresponding to saddle point problems is highly relevant to physical systems such as Stokes flow. However, scaling up numerical solvers for such systems is often met with challenges in efficiency and convergence. Multigrid is an approach with excellent applicability to elliptic problems such as the Stokes equations, and can be a solution to such challenges of scalability and efficiency. The degree of success of such methods, however, is highly contingent on the design of key components of the multigrid scheme, including the hierarchy of discretizations, and the relaxation scheme used. Additionally, in many practical cases, it may be more effective to use a multigrid scheme as a preconditioner to an iterative Krylov subspace solver, as opposed to striving for maximum efficacy of the relaxation scheme in all foreseeable settings. In this paper, we propose an efficient symmetric multigrid preconditioner for the Stokes Equations on a staggered finite difference discretization. Our contribution is focused on crafting a preconditioner that (a) is symmetric indefinite, matching the property of the Stokes system itself, (b) is appropriate for preconditioning the SQMR iterative scheme [1], and (c) has the requisite symmetry properties to be used in this context. In addition, our design is efficient in terms of computational cost and facilitates scaling to large domains.

Hybrid CPU-GPU solution to regularized divergence-free curl-curl equations for electromagnetic inversion problems

The Curl-Curl equation is the foundation of time-harmonic electromagnetic (EM) problems in geophysics. The efficiency of its solution is key to EM simulations, accounting for over 95% of the computation cost in geophysical inversions for magnetotelluric or controlled-source EM problems. However, most published EM inversion codes are still central processing unit (CPU)-based and cannot utilize recent computational developments on the graphic processing units (GPUs). Based on a previously proposed divergence-free algorithm developed on CPUs, this study demonstrates the current limits of the CPU-based inversion procedure. To exploit the high throughput capability of GPUs, we propose a hybrid CPU-GPU framework to solve forward and adjoint problems required for EM inversions. The large sparse linear systems arising from the staggered-grid finite difference approximation of the Curl-Curl equation are solved with a mixed-precision Krylov subspace solver implemented on a GPU. The algorithm is then tested in EM forward and adjoint calculations, with real-world three-dimensional numerical examples. Test results show promising 30× kernel-level speed-ups over the conventional CPU algorithm. This approach may further take the complex frequency domain EM inversions onto the next, practical stage on small affordable GPU platforms.

Explicit and implicit large eddy simulations of variable-density low-speed flows

The application of the Non-oscillatory Forward-in-Time (NFT) class of solvers to facilitate conventional explicit Large Eddy Simulations (LES) and implicit LES (ILES) of low-speed turbulent flow encompassing substantial local density variations is explored. Schemes solving two different sets of Low Mach Number (LMN) equations are proposed with the aim of capturing variable-density effects arising from either a non-isothermal distribution in single-species flow or compositional variation occurring in a binary-species mixing under isothermal conditions. Both schemes employ the semi-implicit NFT integrators based on the Multidimensional Positive Definite Advection Transport Algorithm and a non-symmetric Krylov subspace solver for the variable-coefficient Poisson problems arising from different treatments applied to the non-zero velocity divergence term in the two sets of the LMN equations. The developments are successfully validated using three diverse test cases, differentially heated cavity flow, non-isothermal free-jets, and a helium plume. The reported simulations show that both the dynamic Smagorinsky model in LES and ILES subgrid-scale treatments yield similar accuracy in capturing turbulent effects for the considered flows.

An efficient fourth-order accurate conservative scheme for Riesz space fractional Schrödinger equation with wave operator

In this study, we investigate a high-order accurate conservative finite difference scheme by utilizing a fourth-order fractional central finite difference method for the two-dimensional Riesz space-fractional nonlinear Schrödinger wave equation. The conservation laws of the discrete difference scheme are shown. Meanwhile, the exactness, uniqueness, and prior estimate of the numerical solution are rigorously established. Then, it is proved that the proposed scheme is unconditionally convergent in the discrete L2 and Hγ/2 norm, where γ is a fractional order. Furthermore, we demonstrate that when the fractional order γ and the spatial grid number J increase, the block-Toeplitz coefficient matrix generated by the spatial discretization becomes ill-conditioned. As a result, we adopt an effective linearized iteration method for the nonlinear system, allowing it to be solved efficiently by the Krylov subspace solver with an appropriate circulant preconditioner, in which the fast Fourier transform is applied to speed up the computational cost at each iterative step. Finally, numerical experiments are presented to validate the theoretical findings and the efficiency of the fast algorithm.

A second-order ADI method for pricing options under fractional regime-switching models

<abstract><p>Fractional regime-switching option models have recently attracted much attention because they can capture the sudden state movement of the market, and deal with the non-stationary behavior. A second-order numerical scheme is proposed to solve the regime-switching option pricing models with fractional derivatives in space. The sufficient conditions of the stability and convergence of the proposed scheme are studied in details. An alternating direction implicit (ADI) method is implemented to accelerate the computation in every time layer. Numerical experiments are presented to verify the convergence and efficiency of the proposed method, compared with classical Krylov subspace solvers.</p></abstract>

Heat Conduction with Krylov Subspace Method Using FEniCSx

The study of heat transfer deals with the determination of the rate of heat energy transfer from one system to another driven by a temperature gradient. It can be observed in many natural phenomena and is often the fundamental principle behind several engineering systems. Heat transfer analysis is necessary while designing any product. The most common numerical method used to analyze heat transfer is the finite element method. This paper uses the finite element method to demonstrate steady and transient heat conduction in a three-dimensional bracket. The goal here was to determine the temperature distribution and rate of heat flow in the solid. This is crucial in designing machine elements as they are subjected to various thermal loads during operation and also due to fluctuations in the surrounding environmental conditions. The temperature significantly affects stress, displacements, and volumetric strains. Thus, to analyze thermal stresses induced in a machine element, it is necessary to find the temperature field first. The thermal analysis was performed using the open-source package FEniCSx on Python. The program was run using a preconditioned Krylov subspace method for higher-order function spaces. The Krylov subspace solver drastically reduces computational time. The time taken for the execution of each order was recorded and presented.

Mixed-Precision for Linear Solvers in Global Geophysical Flows.

Semi‐implicit (SI) time‐stepping schemes for atmosphere and ocean models require elliptic solvers that work efficiently on modern supercomputers. This paper reports our study of the potential computational savings when using mixed precision arithmetic in the elliptic solvers. Precision levels as low as half (16 bits) are used and a detailed evaluation of the impact of reduced precision on the solver convergence and the solution quality is performed. This study is conducted in the context of a novel SI shallow‐water model on the sphere, purposely designed to mimic numerical intricacies of modern all‐scale weather and climate (W&C) models. The governing algorithm of the shallow‐water model is based on the non‐oscillatory MPDATA methods for geophysical flows, whereas the resulting elliptic problem employs a strongly preconditioned non‐symmetric Krylov‐subspace Generalized Conjugated‐Residual (GCR) solver, proven in advanced atmospheric applications. The classical longitude/latitude grid is deliberately chosen to retain the stiffness of global W&C models. The analysis of the precision reduction is done on a software level, using an emulator, whereas the performance is measured on actual reduced precision hardware. The reduced‐precision experiments are conducted for established dynamical‐core test‐cases, like the Rossby‐Haurwitz wavenumber 4 and a zonal orographic flow. The study shows that selected key components of the elliptic solver, most prominently the preconditioning and the application of the linear operator, can be performed at the level of half precision. For these components, the use of half precision is found to yield a speed‐up of a factor 4 compared to double precision for a wide range of problem sizes.

Block preconditioners for linear systems in interior point methods for convex constrained optimization

In this paper, we address the preconditioned iterative solution of the saddle-point linear systems arising from the (regularized) Interior Point method applied to linear and quadratic convex programming problems, typically of large scale. Starting from the well-studied Constraint Preconditioner, we review a number of inexact variants with the aim to reduce the computational cost of the preconditioner application within the Krylov subspace solver of choice. In all cases we illustrate a spectral analysis showing the conditions under which a good clustering of the eigenvalues of the preconditioned matrix can be obtained, which foreshadows (at least in case PCG/MINRES Krylov solvers are used), a fast convergence of the iterative method. Results on a set of large size optimization problems confirm that the Inexact variants of the Constraint Preconditioner can yield efficient solution methods.

A suite of Richardson preconditioners for semi-implicit all-scale atmospheric models

The paper documents a suite of preconditioners for Krylov-subspace solvers of elliptic boundary-value problems (BVPs) that underlie semi-implicit integrations of the nonhydrostatic equations governing the dynamics of all-scale atmospheric flows. Effective preconditioning of the linear operators inherent in the semi-implicit models lies at the heart of the state-of-the-art multiscale-flow simulation. This is especially evident in simulations of global weather and climate—posed on a thin spherical shell—where some form of direct tridiagonal inversion of the operator in the vertical is crucial to relax the often enormous stiffness of the problem. The documented preconditioners stem from the Richardson's (1910) idea of augmenting an elliptic BVP with a transient diffusion equation. Exploiting this idea for mixed explicit-implicit pseudo-time-stepping schemes leads to a broad suite of stationary-iteration solvers, including the many classical algorithms. Here, the high-performance all-scale EULAG model (Smolarkiewicz et al. (2014) [58]), with a flexible three-dimensional decomposition of MPI tasks, is furnished with the preconditioners akin to the classical alternating-direction-implicit (ADI) algorithms, generalized to optional permutations of parallel tridiagonal inversions. The utility of various options is found to be problem dependent, in terms of computational accuracy as well as efficiency. The main thrust of the work is on the long-range forecasts using large anisotropic grids. The relative efficiency and/or accuracy gains attainable with the developed preconditioners are illustrated for idealized scenarios representative of atmospheric flows from planetary to a single-cloud and laboratory scales. The key insight that best encapsulates the significance and novelty of the present work is that there is no single “super-preconditioner” that will perform best in all cases, yet the suite as a whole offers substantial gains in the model performance.

Fast Solution of 3-D Eddy-Current Problems in Multiply Connected Domains by a, v-φ and t-φ Formulations With Multigrid-Based Algorithm for Cohomology Generation

The fast solution of three-dimensional eddy current problems is still an open problem, especially when real-size finite element models with millions of degrees of freedom are considered. In order to lower the number of degrees of freedom a magnetic scalar potential can be used in the insulating parts of the model. This may become difficult when the model geometry presents some conductive parts which are multiply connected. In this work a multigrid-based algoritm is proposed that allows for a calculation in linear-time of cohomology, which is needed to introduce the scalar potential without cuts. This algorithm relies on an algebraic multigrid solver for curl-curl field problems, which ensures optimal computational complexity. Numerical results show that the novel algorithm outperforms state-of-the-art methods for cohomology generation based on homological algebra. In addition, based on this algoritm, novel <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$a,v$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\varphi $ </tex-math></inline-formula> and <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$t$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\varphi $ </tex-math></inline-formula> formulations to analyze three-dimensional eddy current problems in multiply connected domains are proposed. Both formulations, after discretization by the cell method, lead to a complex symmetric system of linear equations amenable to fast iterative solution by Krylov-subspace solvers. These formulations are able to provide very accurate numerical results, with a minimum amount of degrees of freedoms to represent the eddy current model. In this way the computational performance is improved compared to the classical <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$A,V$ </tex-math></inline-formula> - <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$A$ </tex-math></inline-formula> formulation typically implemented in finite element software for electromagnetic design.

Решение анизотропных уравнений теплопроводности экспоненциальными схемами на крыловских полиномиальных подпространствах и подпространствах «сдвиг-обращение»

We assess performance of the exponential Krylov subspace methods for solving a class of parabolic problems with a strong anisotropy in coefficients. Different boundary conditions are considered, which have a direct impact on the smallest eigenvalue of the discretized operator and, hence, on the convergence behavior of the exponential Krylov subspace solvers. Restarted polynomial Krylov subspace methods and shift-and-invert Krylov subspace methods combined with algebraic multigrid are considered.

Solving anisotropic heat equations by exponential shift-and-invert and polynomial Krylov subspace methods

We assess performance of the exponential Krylov subspace methods for solving a class of parabolic problems with a strong anisotropy in coefficients. Different boundary conditions are considered, which have a direct impact on the smallest eigenvalue of the discretized operator and, hence, on the convergence behavior of the exponential Krylov subspace solvers. Restarted polynomial Krylov subspace methods and shift-and-invert Krylov subspace methods combined with algebraic multigrid are considered.

On convergence of a structure preserving difference scheme for two-dimensional space-fractional nonlinear Schrödinger equation and its fast implementation

In this paper we intend to construct a structure preserving difference scheme for two-dimensional space-fractional nonlinear Schrödinger (2D SFNS) equation with the integral fractional Laplacian. The temporal direction is discretized by the modified Crank-Nicolson method, and the spatial variable is approximated by a novel fractional central difference method. The mass and energy conservations and the convergence are rigorously proved for the proposed scheme. For 1D SFNS equation, the convergence relies heavily on the L∞-norm boundness of the numerical solution of the proposed scheme. However, we cannot obtain the L∞-norm boundness of the numerical solution by using the similar process for the 2D SFNS equation. One of the major significance of this paper is that we first obtain the L∞-norm boundness of the numerical solution and L2-norm error estimate via the popular “cut-off” function for the 2D SFNS equation. Further, we reveal that the spatial discretization generates a block-Toeplitz coefficient matrix, and it will be ill-conditioned as the spatial grid mesh number M and the fractional order α increase. Thus, we exploit an linearized iteration algorithm for the nonlinear system, so that it can be efficiently solved by the Krylov subspace solver with a suitable preconditioner, where the 2D fast Fourier transform (2D FFT) is applied in the solver to accelerate the matrix-vector product, and the standard orthogonal projection approach is used to eliminate the drift of mass and energy. Extensive numerical results are reported to confirm the theoretical analysis and high efficiency of the proposed algorithm.

A Note on Parallel Preconditioning for the All-at-Once Solution of Riesz Fractional Diffusion Equations

The $p$-step backwards difference formula (BDF) for solving the system of ODEs can result in a kind of all-at-once linear systems, which are solved via the parallel-in-time preconditioned Krylov subspace solvers (see McDonald, Pestana, and Wathen [SIAM J. Sci. Comput., 40(2) (2018): A1012-A1033] and Lin and Ng [arXiv:2002.01108, 17 pages]. However, these studies ignored that the $p$-step BDF ($p\geq 2$) is not selfstarting, when they are exploited to solve time-dependent PDEs. In this note, we focus on the 2-step BDF which is often superior to the trapezoidal rule for solving the Riesz fractional diffusion equations, but its resultant all-at-once discretized system is a block triangular Toeplitz system with a low-rank perturbation. Meanwhile, we first give an estimation of the condition number of the all-at-once systems and then adapt the previous work to construct two block circulant (BC) preconditioners. Both the invertibility of these two BC preconditioners and the eigenvalue distributions of preconditioned matrices are discussed in details. The efficient implementation of these BC preconditioners is also presented especially for handling the computation of dense structured Jacobi matrices. Finally, numerical experiments involving both the one- and two-dimensional Riesz fractional diffusion equations are reported to support our theoretical findings.

Marvin: A parallel three-dimensional transport code based on the discrete ordinates method for reactor shielding calculations

Marvin is a parallel three-dimensional discrete ordinates (SN) code dedicated for reactor shielding calculations. It employs the multi-group and SN approximations to discretize the energy and angular variables, respectively. Multiple discretization methods, namely weighted diamond difference, theta-weighted diamond difference and step-characteristic, have been implemented to discretize the spatial variables. The within-group transport equation is solved by the source iteration or Krylov subspace solvers, both of which employ the Koch-Baker-Alcouffe (KBA) algorithm to parallelize the transport sweep process. The upscattering in the thermal energy range is treated with thermal upscattering iteration scheme. The code design of Marvin features the mixed programming model and a three-layer architecture design. Two benchmark problems devised from realistic reactor core configurations, namely the VENUS-3 and HBR-2 problems, were calculated to evaluate the correctness, accuracy and parallel computing performance of Marvin. The calculation results were compared against TORT solutions as well as measured data. The validation results indicate that: Marvin is properly implemented and can faithfully produce solutions to the reactor shielding problems; the computing efficiency of Marvin can be significantly enhanced via the KBA parallel transport sweep algorithm.

Robust stabilised finite element solvers for generalised Newtonian fluid flows

Various materials and solid-fluid mixtures of engineering and biomedical interest can be modelled as generalised Newtonian fluids, as their apparent viscosity depends locally on the flow field. Despite the particular features of such models, it is common practice to combine them with numerical techniques originally conceived for Newtonian fluids, which can bring several issues such as spurious pressure boundary layers, unsuitable natural boundary conditions and coupling terms spoiling the efficiency of nonlinear solvers and preconditioners. In this work, we present a finite element framework dealing with such issues while maintaining low computational cost and simple implementation. The building blocks of our algorithm are (i) an equal-order stabilisation method preserving consistency even for lowest-order discretisations, (ii) robust extrapolation of velocities in the time-dependent case to decouple the rheological law from the overall system, (iii) adaptive time step selection and (iv) a fast physics-based preconditioned Krylov subspace solver, to tackle the relevant range of discretisation parameters including highly varying viscosity. Selected numerical experiments are provided demonstrating the potential of our approach in terms of robustness, accuracy and efficiency for problems of practical interest.

Fast dissipation-preserving difference scheme for nonlinear generalized wave equations with the integral fractional Laplacian

In this paper, we construct a dissipation-preserving difference scheme for two-dimensional nonlinear generalized wave equations with the integral fractional Laplacian. We discuss the discrete dissipation property of the scheme in detail, and we also analyze the existence, uniqueness and the unconditional convergence of the proposed scheme. Further, we reveal that the spatial discretization generates a block-Toeplitz coefficient matrix, and it will be ill-conditioned as the spatial grid mesh number M and the fractional order α increase. Thus, we exploit an efficient linearized iteration algorithm for the nonlinear system, such that it can be efficiently solved by the Krylov subspace solver with a suitable preconditioner, where the two-dimensional fast Fourier transform is used in the solver to accelerate the matrix-vector product. Extensive numerical experiments are provided to verify the theoretical analysis and discrete dissipation law of the proposed scheme in long-time computations.

Comparison the iterative solvers for large sparse matrix in 3D electromagnetic forward modelling

In 3D electromagnetic (EM) forward modeling, an analytical solution is generally not available. Numerical solution is commonly applied to solve the forward modeling problems, mostly based on iterative solvers. The efficiency of EM forward modeling is critical for the development of practical inversion for EM data. The Krylov subspace solvers are widely used to solve frequency-domain EM forward modeling problems. However, these solvers converge remarkably more slowly as the operating period increases. This can be improved by the use of preconditioner and divergence correction. Multigrid (MG) solver is efficient for solving EM forward modelling problems without the use of preconditioner and divergence correction. In this paper, a MG solver is compared with Bi-Conjugate Gradients Stabilized (BCG) solvers with different preconditioners. They are compared, in terms of iteration number and computing time, indicating the MG solver is much more efficient.

Optimized Structured Preconditioner for a Class of Complex Symmetric Linear Systems

An optimized structured (OS) preconditioner for a class of complex symmetric linear systems is proposed. The OS-preconditioner results in a fast Krylov subspace solver, which is robust with respect to the mesh-size. The computational complexity of OS-preconditioner is analyzed and it is shown that all eigenvalues of the corresponding preconditioned matrix are positive real and distributed on[$\frac{1}{2}$+$\frac{\varepsilon }{{2\sqrt {1 + {\varepsilon ^2}} }}$,1]. Numerical experiments confirming the theoretical derivations are presented to verify the effectiveness and stability of the OS-preconditioner.

A statically condensed discontinuous Galerkin spectral element method on Gauss-Lobatto nodes for the compressible Navier-Stokes equations

We present a static-condensation method for time-implicit discretizations of the Discontinuous Galerkin Spectral Element Method on Gauss-Lobatto points (GL-DGSEM). We show that, when solving the compressible Navier-Stokes equations, it is possible to reorganize the linear system that results from the implicit time-integration of the GL-DGSEM as a Schur complement problem, which can be efficiently solved using static condensation. The use of static condensation reduces the linear system size and improves the condition number of the system matrix, which translates into shorter computational times when using direct and iterative solvers.The statically condensed GL-DGSEM presented here can be applied to linear and nonlinear advection-diffusion partial differential equations in conservation form with the standard and the split-form variants of the DGSEM. To test it we solve the compressible Navier-Stokes equations with direct and Krylov subspace solvers, and we show for a selected problem that using the statically condensed GL-DGSEM leads to speed-ups of up to 200 when compared to the time-explicit GL-DGSEM, and speed-ups of up to three when compared with the time-implicit GL-DGSEM that solves the global system.The GL-DGSEM has gained increasing popularity in recent years because it satisfies the summation-by-parts property, which enables the construction of provably entropy stable schemes, and because it is computationally very efficient. In this paper, we show that the GL-DGSEM has an additional advantage: It can be statically condensed.