Preconditioned Krylov Subspace Methods Research Articles

A Preconditioned Krylov Subspace Method for Linear Inverse Problems with General-Form Tikhonov Regularization

A Preconditioned Krylov Subspace Method for Linear Inverse Problems with General-Form Tikhonov Regularization

A matrix-free parallel two-level deflation preconditioner for two-dimensional heterogeneous Helmholtz problems

We propose a matrix-free parallel two-level deflation method combined with the Complex Shifted Laplacian Preconditioner (CSLP) for two-dimensional heterogeneous Helmholtz problems encountered in seismic exploration, antennas, and medical imaging. These problems pose challenges in terms of accuracy and convergence due to scalability issues with numerical solvers. Motivated by the limitations imposed by excessive computational time and memory constraints when employing a sequential solver with constructed matrices, we parallelize the two-level deflation method without constructing any matrices. Our approach utilizes preconditioned Krylov subspace methods and approximates the CSLP preconditioner with a parallel geometric multigrid V-cycle. For the two-level deflation, standard inter-grid deflation vectors and further high-order deflation vectors are considered. As another main contribution, the matrix-free Galerkin coarsening approach and a novel re-discretization scheme as well as high-order finite-difference schemes on the coarse grid are studied to obtain wavenumber-independent convergence. The optimal settings for an efficient coarse-grid problem solver are investigated. Numerical experiments of model problems show that the wavenumber independence has been obtained for medium wavenumbers. The matrix-free parallel framework shows satisfactory weak and strong parallel scalability.

Approximate inner solvers for block preconditioning of the incompressible Navier–Stokes problems discretized by isogeometric analysis

AbstractWe deal with efficient numerical solution of the steady incompressible Navier–Stokes equations (NSE) using our in‐house solver based on the isogeometric analysis (IgA) approach. We are interested in the solution of the arising saddle‐point linear systems using preconditioned Krylov subspace methods. Based on our comparison of ideal versions of several state‐of‐the‐art block preconditioners for linear systems arising from the IgA discretization of the incompressible NSE, suitable candidates have been selected. In the present paper, we focus on selecting efficient approximate solvers for solving subsystems within these preconditioning methods. We investigate the impact on the convergence of the outer solver and aim to identify an effective combination. For this purpose, we compare convergence properties of the selected solution approaches for problems with different viscosity values, mesh refinement levels and discretization bases.

A new block preconditioner for weighted Toeplitz regularized least-squares problems

We introduce a new block preconditioner for the solution of weighted Toeplitz regularized least-squares problems written in augmented system form. The proposed preconditioner is obtained based on the new splitting of coefficient matrix which results in an unconditionally convergent stationary iterative method. Spectral analysis of the preconditioned matrix is investigated. In particular, we show that the preconditioned matrix has a very nice eigenvalue distribution which can lead to fast convergence of the preconditioned Krylov subspace methods such as GMRES. Numerical experiments are reported to demonstrate the performance of preconditioner used with GMRES method in the solution of augmented system form of weighted Toeplitz regularized least-squares problems.

Multiscale steady discrete unified gas kinetic scheme with macroscopic coarse mesh acceleration using preconditioned Krylov subspace method for multigroup neutron Boltzmann transport equation.

A multiscale steady discrete unified gas kinetic scheme with macroscopic coarse mesh acceleration [accelerated steady discrete unified gas kinetic scheme (SDUGKS)] is proposed to improve the convergence of the original SDUGKS for an optically thick system in solving the multigroup neutron Boltzmann transport equation (NBTE) to analyze the distribution of fission energy in the reactor core. In the accelerated SDUGKS, by solving the coarse mesh macroscopic governing equations (MGEs) derived from the moment equations of the NBTE, the numerical solutions of the NBTE on fine meshes at the mesoscopic level can be rapidly obtained from the prolongation of the coarse mesh solutions of the MGE. Furthermore, the use of the coarse mesh can greatly reduce the computational variables and improve the computational efficiency of the MGE. The biconjugate gradient stabilized Krylov subspace method with the modified incomplete LU preconditioner and the lower-upper symmetric-Gauss-Seidel sweeping method are implemented to solve the discrete systems of the macroscopic coarse mesh acceleration model and mesoscopic SDUGKS to further improve the numerical efficiency. Numerical solutions validate good numerical accuracy and high acceleration efficiency of the proposed accelerated SDUGKS for the complicated multiscale neutron transport problems.

Exploiting spatial symmetries for solving Poisson's equation

This paper presents a strategy to accelerate virtually any Poisson solver by taking advantage of s spatial reflection symmetries. More precisely, we have proved the existence of an inexpensive block diagonalisation that transforms the original Poisson equation into a set of 2s fully decoupled subsystems then solved concurrently. This block diagonalisation is identical regardless of the mesh connectivity (structured or unstructured) and the geometric complexity of the problem, therefore applying to a wide range of academic and industrial configurations. In fact, it simplifies the task of discretising complex geometries since it only requires meshing a portion of the domain that is then mirrored implicitly by the symmetries' hyperplanes. Thus, the resulting meshes naturally inherit the exploited symmetries, and their memory footprint becomes 2s times smaller. Thanks to the subsystems' better spectral properties, iterative solvers converge significantly faster. Additionally, imposing an adequate grid points' ordering allows reducing the operators' footprint and replacing the standard sparse matrix-vector products with the sparse matrix-matrix product, a higher arithmetic intensity kernel. As a result, matrix multiplications are accelerated, and massive simulations become more affordable. Finally, we include numerical experiments based on a turbulent flow simulation and making state-of-the-art solvers exploit a varying number of symmetries. On the one hand, algebraic multigrid and preconditioned Krylov subspace methods require up to 23% and 72% fewer iterations, resulting in up to 1.7x and 5.6x overall speedups, respectively. On the other, sparse direct solvers' memory footprint, setup and solution costs are reduced by up to 48%, 58% and 46%, respectively.

MGM: A Meshfree Geometric Multilevel Method for Systems Arising from Elliptic Equations on Point Cloud Surfaces

We develop a new meshfree geometric multilevel (MGM) method for solving linear systems that arise from discretizing elliptic PDEs on surfaces represented by point clouds. The method uses a Poisson disk sampling-type technique for coarsening the point clouds and new meshfree restriction/interpolation operators based on polyharmonic splines for transferring information between the coarsened point clouds. These are then combined with standard smoothing and operator coarsening methods in a V-cycle iteration. MGM is applicable to discretizations of elliptic PDEs based on various localized meshfree methods, including RBF finite differences (RBF-FD) and generalized finite differences (GFD). We test MGM both as a standalone solver and preconditioner for Krylov subspace methods on several test problems using RBF-FD and GFD and numerically analyze convergence rates, efficiency, and scaling with increasing point cloud sizes. We also perform a side-by-side comparison to algebraic multigrid methods for solving the same systems. Finally, we further demonstrate the effectiveness of MGM by applying it to three challenging applications on complicated surfaces: pattern formation, surface harmonics, and geodesic distance.

Efficient ADI schemes and preconditioning for a class of high-dimensional spatial fractional diffusion equations with variable diffusion coefficients

In this paper, alternating direction implicit (ADI) finite difference method and preconditioned Krylov subspace method are combined to solve a class of high-dimensional spatial fractional diffusion equations with variable diffusion coefficients. We prove the unconditional stability and convergence rate of the ADI finite difference method provided that the diffusion coefficients satisfy the given conditions. For the linear system in each spatial direction, we establish a circulant approximate inverse preconditioner to accelerate the Krylov subspace method. In addition, we also use matrix-free algorithms and fast Fourier transforms (FFT) to speed up the solution of linear systems. Numerical experiments show the utility of the ADI method and the preconditioner.

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.

ParGeMSLR: A parallel multilevel Schur complement low-rank preconditioning and solution package for general sparse matrices

This paper discusses parGeMSLR, a C++/MPI software library for the solution of sparse systems of linear algebraic equations via preconditioned Krylov subspace methods in distributed-memory computing environments. The preconditioner implemented in parGeMSLR is based on algebraic domain decomposition and partitions the symmetrized adjacency graph recursively into several non-overlapping partitions via a p-way vertex separator, where p is an integer multiple of the total number of MPI processes. From a numerical perspective, parGeMSLR builds a Schur complement approximate inverse preconditioner as the sum between the matrix inverse of the interface coupling matrix and a low-rank correction term. To reduce the cost associated with the computation of the approximate inverse matrices, parGeMSLR exploits a multilevel partitioning of the algebraic domain. The parGeMSLR library is implemented on top of the Message Passing Interface and can solve both real and complex linear systems. Furthermore, parGeMSLR can take advantage of hybrid computing environments with in-node access to one or more Graphics Processing Units. Finally, the parallel efficiency (weak and strong scaling) of parGeMSLR is demonstrated on a few model problems arising from discretizations of 3D Partial Differential Equations.

A stochastic Galerkin method with adaptive time-stepping for the Navier–Stokes equations

We study the time-dependent Navier–Stokes equations in the context of stochastic finite element discretizations. Specifically, we assume that the viscosity is a random field given in the form of a generalized polynomial chaos expansion, and we use the stochastic Galerkin method to extend the methodology from Kay et al. (2010) [21] into this framework. For the resulting stochastic problem, we explore the properties of the resulting stochastic solutions, and we also compare the results with that of Monte Carlo and stochastic collocation. Since the time-stepping scheme is fully implicit, we also propose strategies for efficient solution of the stochastic Galerkin linear systems using a preconditioned Krylov subspace method. The effectiveness of the stochastic Galerkin method is illustrated by numerical experiments.

Modelling a permanent magnet synchronous motor in FEniCSx for parallel high-performance simulations

There are concerns that the extreme requirements of heavy-duty vehicles and aviation will see them left behind in the electrification of the transport sector, becoming the most significant emitters of greenhouse gases. Engineers extensively use the finite element method to analyse and improve the performance of electric machines, but new highly scalable methods with a linear (or near) time complexity are required to make extreme-scale models viable. This paper introduces a three-dimensional permanent magnet synchronous motor model using FEniCSx, a finite element platform tailored for efficient computing and data handling at scale. The model demonstrates comparable magnetic flux density distributions to a verification model built in Ansys Maxwell with a maximum deviation of 7% in the motor’s static regions. Solving the largest mesh, comprising over eight million cells, displayed a speedup of 198 at 512 processes. A preconditioned Krylov subspace method was used to solve the system, requiring 92% less memory than a direct solution. It is expected that advances built on this approach will allow system-level multiphysics simulations to become feasible within electric machine development. This capability could provide the near real-world accuracy needed to bring electric propulsion systems to large vehicles. • Stable and scalable electric machine analysis is possible using iterative methods. • Open-source solvers can leverage the advantages of high-performance computing. • The parallel performance of FEniCSx exceeds an established electromagnetic code.

A Comparison of Linear Solvers for Resolving Flow in Three‐Dimensional Discrete Fracture Networks

AbstractWe compare various methods for resolving steady flow within three‐dimensional discrete fracture networks, including direct methods, Krylov subspace methods with and without preconditioning, and multi‐grid methods. We compared the performance of the methods based on compute times and scaling of the solution as a function of the number of grid nodes and log‐variance of the hydraulic aperture. The methods are applied to three test cases: (a) variable density of networks with a truncated power‐law distribution of fracture lengths, (b) a fixed network composed of monodisperse fracture sizes but varied permeability/aperture heterogeneity, (c) and a network based on field site in Nevada, US. We chose these cases to allow us to study the impact of the mesh size and flow properties, as well as to demonstrate our conclusions on a large‐scale, realistic problem (more than 40 million mesh nodes). A direct solution using Cholesky factorization outperformed other methods for every example but was closely followed in performance by some algebraic multigrid (AMG) preconditioned Krylov subspace methods. Among the Krylov methods, conjugate gradients (CG) with an AMG preconditioner performs the best. Generally, Cholesky factorization is recommended, but CG with an AMG preconditioner may be suitable for very large problems beyond 40 million nodes where the entire linear system cannot reside in memory.

A Fast Preconditioned Semi-Implicit Difference Scheme for Strongly Nonlinear Space-Fractional Diffusion Equations

In this paper, we propose a semi-implicit difference scheme for solving one-dimensional nonlinear space-fractional diffusion equations. The method is first-order accurate in time and second-order accurate in space. It uses a fractional central difference formula and the backward Euler method to approximate its space and time derivatives, respectively. Stability and convergence properties of the proposed scheme are proved with the help of a discrete Grönwall inequality. Moreover, we extend the method to the solution of two-dimensional nonlinear models. A fast matrix-free implementation based on preconditioned Krylov subspace methods is presented for solving the discretized linear systems. The resulting fast preconditioned semi-implicit difference scheme reduces the memory requirement of conventional semi-implicit difference schemes from O(Ns2) to O(Ns) and the computational complexity from O(Ns3) to O(NslogNs) in each iterative step, where Ns is the number of space grid points. Experiments with two numerical examples are shown to support the theoretical findings and to illustrate the efficiency of our proposed method.

Preconditioned iterative method for nonsymmetric saddle point linear systems

In this paper, a new preconditioned iterative method is presented to solve a class of nonsymmetric nonsingular or singular saddle point problems. The implementation of the proposed preconditioned Krylov subspace method avoids solving inverse of Schur complement and only needs to solve one linear sub-system at each step, which implies that it may save considerable costs. Theoretical convergence analysis, including the bounds of eigenvalues and eigenvectors, the degree of the minimal polynomial of the preconditioned matrix, are discussed in details. Moreover, a novel algebraic estimation technique for finding a practical iteration parameter is presented, which is very effective and practical even for large scale problems. At last, some numerical examples are carried, showing that the theoretical results are valid and convincing.

Fast second-order implicit difference schemes for time distributed-order and Riesz space fractional diffusion-wave equations

In this paper, fast numerical methods are established to solve a class of time distributed-order and Riesz space fractional diffusion-wave equations. We derive new difference schemes by a weighted and shifted Grünwald formula in time and a fractional centered difference formula in space. The unconditional stability and second-order convergence in time, space and distributed-order of the difference schemes are analyzed. In the one-dimensional case, the Gohberg-Semencul formula utilizing a preconditioned Krylov subspace method is developed to solve the Toeplitz linear system derived from the proposed difference scheme. In the two-dimensional case, we also design a global preconditioned conjugate gradient method with a truncated preconditioner to solve the resulting Sylvester matrix equations. We prove that the spectrums of the preconditioned matrices in both cases are clustered around 1, such that the proposed numerical methods with preconditioners converge very quickly. Some numerical experiments are carried out to demonstrate the effectiveness of the proposed difference schemes and show that the performances of the proposed fast solution algorithms are better than other testing methods.

Highly accurate, linear, and unconditionally energy stable large time-stepping schemes for the Functionalized Cahn–Hilliard gradient flow equation

This paper aims to develop a class of large time-stepping schemes for solving the Functionalized Cahn–Hilliard (FCH) gradient flow equation. By using the Invariant Energy Quadratization (IEQ) approach, we successfully construct several highly efficient linear schemes that allow large time step sizes to be used in computations. All the developed schemes with a fixed time step size are proved to be unconditionally uniquely solvable and unconditionally energy stable in theory. Moreover, we show that the obtained semi-discrete systems can be solved efficiently by using the preconditioned Krylov subspace methods, since all the linear operators of these systems are proved to be symmetric and positive definite. Finally, several numerical experiments are carried out to verify the accuracy, efficiency and energy stability of the developed IEQ schemes. We observe that our IEQ schemes are more accurate than the Scalar Auxiliary Variable (SAV) schemes. This can be attributed to the fact that the IEQ method can preserve the operator structures of the original equation as much as possible.

Towards accuracy and scalability: Combining Isogeometric Analysis with deflation to obtain scalable convergence for the Helmholtz equation

Finding fast yet accurate numerical solutions to the Helmholtz equation remains a challenging task. The pollution error (i.e. the discrepancy between the numerical and analytical wave number k) requires the mesh resolution to be kept fine enough to obtain accurate solutions. A recent study showed that the use of Isogeometric Analysis (IgA) for the spatial discretization significantly reduces the pollution error. However, solving the resulting linear systems by means of a direct solver remains computationally expensive when large wave numbers or multiple dimensions are considered. An alternative lies in the use of (preconditioned) Krylov subspace methods. Recently, the use of the exact Complex Shifted Laplacian Preconditioner (CSLP) with a small complex shift has shown to lead to wave number independent convergence while obtaining more accurate numerical solutions using IgA.In this paper, we propose the use of deflation techniques combined with an approximated inverse of the CSLP using a geometric multigrid method. Numerical results obtained for one- and two-dimensional model problems, including constant and non-constant wave numbers, show scalable convergence with respect to the wave number and approximation order p of the spatial discretization. Furthermore, when kh is kept constant, the proposed approach leads to a significant reduction of the computational time compared to the use of the multigrid-approximated or exact inverse of the CSLP with a small shift, in particular for three-dimensional model problems.

A parallel multilevel domain decomposition method for source identification problems governed by elliptic equations

In this paper we develop a parallel multilevel domain decomposition method for large-scale source identification problems governed by elliptic equations. A popular approach is to formulate the inverse problem as a PDE-constrained optimization problem. The minima satisfies a Karush–Kuhn–Tucker (KKT) system consisting of the state, adjoint and source equations which is rather difficult to solve on parallel computers. We propose and study a parallel method that decomposes the optimization problem on the global domain into subproblems on overlapping subdomains, each subdomain is further decomposed to form an additive Schwarz preconditioner for solving these smaller subproblems simultaneously with a preconditioned Krylov subspace method. For each subproblem, the overlapping part of the solution is discarded and the remaining non-overlapping part of the solution is put together to obtain an approximated global solution to the inverse problem. Since all the subproblems are solved independently, the multilevel domain decomposition method has the advantage of higher degree of parallelism. Numerical experiments show that the algorithm is accurate in terms of the reconstruction error and has reasonably good speedup in terms of the computing time. The efficiency and robustness of the proposed approach on a parallel computer with more than 1, 000 processors are reported.

On Chebyshev accelerated iteration methods for two-by-two block linear systems

In this paper, we consider applying the Chebyshev acceleration technique to improve the convergence behavior of some iteration methods for solving the two-by-two block linear systems with semidefinite sub-blocks. Practical problem independent approximations about the eigenvalue bounds of the preconditioned matrices facilitate parameter free implementations of the concerned Chebyshev accelerated methods. To analyze their convergence behaviors, we discuss in detail the spectral properties of the corresponding preconditioned matrices, especially their eigenvector distributions, which lead to almost problem independent iteration error bounds. Numerical experiments of the complex symmetric linear systems demonstrate the efficiency of the new Chebyshev acceleration methods compared with the recently published Chebyshev accelerated PMHSS method and some preconditioned Krylov subspace methods.