GMRES Solver Research Articles

Hp-Multigrid Preconditioner for a Divergence-Conforming HDG Scheme for the Incompressible Flow Problems

In this study, we present an hp-multigrid preconditioner for a divergence-conforming HDG scheme for the generalized Stokes and the Navier–Stokes equations using an augmented Lagrangian formulation. Our method relies on conforming simplicial meshes in two- and three-dimensions. The hp-multigrid algorithm is a multiplicative auxiliary space preconditioner that employs the lowest-order space as the auxiliary space, and we develop a geometric multigrid method as the auxiliary space solver. For the generalized Stokes problem, the crucial ingredient of the geometric multigrid method is the equivalence between the condensed lowest-order divergence-conforming HDG scheme and a Crouzeix–Raviart discretization with a pressure-robust treatment as introduced in Linke and Merdon (Comput Methods Appl Mech Engrg 311:304–326, 2022), which allows for the direct application of geometric multigrid theory on the Crouzeix–Raviart discretization. The numerical experiments demonstrate the robustness of the proposed hp-multigrid preconditioner with respect to mesh size and augmented Lagrangian parameter, with iteration counts insensitivity to polynomial order increase. Inspired by the works by Benzi and Olshanskii (SIAM J Sci Comput 28:2095–2113, 2006) and Farrell et al. (SIAM J Sci Comput 41:A3073–A3096, 2019), we further test the proposed preconditioner on the divergence-conforming HDG scheme for the Navier–Stokes equations. Numerical experiments show a mild increase in the iteration counts of the preconditioned GMRes solver with the rise in Reynolds number up to 103\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$10^3$$\\end{document}.

Sensitivity Analysis of Underwater Structural-Acoustic Problems Based on Coupled Finite Element Method/Fast Multipole Boundary Element Method with Non-Uniform Rational B-Splines

For the direct differentiation technique-based modeling of acoustic fluid–structure interaction and structural-acoustic sensitivity analysis, a coupling algorithm based on the finite element method (FEM) and the fast multipole boundary element method (FMBEM) is suggested. By bypassing the challenging volume parameterization process in isogeometric finite element techniques and the laborious meshing process in traditional FEM/BEM, non-uniform rational B-splines (NURBS) isogeometric analysis (IGA) is utilized to immediately perform numerical analysis on CAD models. The matrix-vector products in the boundary element analysis are accelerated using the fast multipole method (FMM). To hasten the solution of the linear system of equations, the iterative solver GMRES is used. The numerical prediction of the effects of arbitrarily shaped vibrating structures on the sound field is made feasible by the FEM/FMBEM technique. A number of numerical examples are provided to show the applicability and effectiveness of the suggested approach.

3-D induction log modelling with integral equation method and domain decomposition pre-conditioning

SUMMARY The deployment of electromagnetic (EM) induction tools while drilling is one of the standard routines for assisting the geosteering decision-making process. The conductivity distribution obtained through the inversion of the EM induction log can provide important information about the geological structure around the borehole. To image the 3-D geological structure in the subsurface, 3-D inversion of the EM induction log is required. Because the inversion process is mainly dependent on forward modelling, the use of a fast and accurate forward modelling tool is essential. In this paper, we present an improved version of the integral equation (IE) based modelling technique for general anisotropic media with domain decomposition pre-conditioning. The discretized IE after domain decomposition equals a fixed-point equation that is solved iteratively with either the block Gauss–Seidel or Jacobi pre-conditioning. Within each iteration, the inverse of the block matrix is computed using a Krylov subspace method instead of a direct solver. An additional reduction in computational time is obtained by using an adaptive relative residual stopping criterion in the iterative solver. Using this domain decomposition scheme, numerical experiments show computation time reductions by factors of 1.97–2.84 compared to solving the full-domain IE with a GMRES solver and a contraction IE pre-conditioner. Additionally, the reduction of memory requirement for covering a large area of the induction tool sensitivity enables acceleration with limited GPU memory. Hence, we conclude that the domain decomposition method is improving the efficiency of the IE method by reducing the computation time and memory requirement.

Solver algorithm for stabilized space-time formulation of advection-dominated diffusion problem

This article shows how to develop an efficient solver for a stabilized numerical space-time formulation of the advection-dominated diffusion transient equation. At the discrete space-time level, we approximate the solution using higher-order continuous B-spline basis functions in both spatial and temporal dimensions. This problem is very difficult to solve numerically using the standard Galerkin finite element method due to artificial oscillations present when the advection term dominates the diffusion term. However, a first-order constraint least-square formulation allows us to obtain numerical solutions avoiding oscillations. The advantages of space-time formulations are the use of high-order methods and the feasibility of developing space-time mesh adaptive techniques on well-defined discrete problems. We develop a solver for a least-square formulation to obtain a stabilized and symmetric problem on finite element meshes. The computational cost of our solver is bounded by the cost of the inversion of the space-time mass and stiffness (with one value fixed at a point) matrices and the cost of the GMRES solver applied for the symmetric and positive definite problem. We illustrate our findings on an advection-dominated diffusion space-time model problem and present two numerical examples: one with isogeometric analysis discretizations and the second one with an adaptive space-time finite element method.

Analytical and Numerical Results for the Diffusion-Reaction Equation When the Reaction Coefficient Depends on Simultaneously the Space and Time Coordinates

We utilize the travelling-wave Ansatz to obtain novel analytical solutions to the linear diffusion–reaction equation. The reaction term is a function of time and space simultaneously, firstly in a Lorentzian form and secondly in a cosine travelling-wave form. The new solutions contain the Heun functions in the first case and the Mathieu functions for the second case, and therefore are highly nontrivial. We use these solutions to test some non-conventional explicit and stable numerical methods against the standard explicit and implicit methods, where in the latter case the algebraic equation system is solved by the preconditioned conjugate gradient and the GMRES solvers. After this verification, we also calculate the transient temperature of a 2D surface subjected to the cooling effect of the wind, which is a function of space and time again. We obtain that the explicit stable methods can reach the accuracy of the implicit solvers in orders of magnitude shorter time.

Compressed basis GMRES on high-performance graphics processing units

Krylov methods provide a fast and highly parallel numerical tool for the iterative solution of many large-scale sparse linear systems. To a large extent, the performance of practical realizations of these methods is constrained by the communication bandwidth in current computer architectures, motivating the investigation of sophisticated techniques to avoid, reduce, and/or hide the message-passing costs (in distributed platforms) and the memory accesses (in all architectures). This article leverages Ginkgo’s memory accessor in order to integrate a communication-reduction strategy into the (Krylov) GMRES solver that decouples the storage format (i.e., the data representation in memory) of the orthogonal basis from the arithmetic precision that is employed during the operations with that basis. Given that the execution time of the GMRES solver is largely determined by the memory accesses, the cost of the datatype transforms can be mostly hidden, resulting in the acceleration of the iterative step via a decrease in the volume of bits being retrieved from memory. Together with the special properties of the orthonormal basis (whose elements are all bounded by 1), this paves the road toward the aggressive customization of the storage format, which includes some floating-point as well as fixed-point formats with mild impact on the convergence of the iterative process. We develop a high-performance implementation of the “compressed basis GMRES” solver in the Ginkgo sparse linear algebra library using a large set of test problems from the SuiteSparse Matrix Collection. We demonstrate robustness and performance advantages on a modern NVIDIA V100 graphics processing unit (GPU) of up to 50% over the standard GMRES solver that stores all data in IEEE double-precision.

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

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

Comparative study of inner–outer Krylov solvers for linear systems in structured and high-order unstructured CFD problems

Advanced Krylov subspace methods are investigated for the solution of large sparse linear systems arising from stiff adjoint-based aerodynamic shape optimization problems. A special attention is paid to the flexible inner–outer GMRES strategy combined with most relevant preconditioning and deflation techniques. The choice of this specific class of Krylov solvers for challenging problems is based on its outstanding convergence properties. Typically in our implementation the efficiency of the preconditioner is enhanced with a domain decomposition method with overlapping. However, maintaining the performance of the preconditioner may be challenging since scalability and efficiency of a preconditioning technique are properties often antagonistic to each other. In this paper we demonstrate how flexible inner–outer Krylov methods are able to overcome this critical issue. A numerical study is performed considering either a Finite Volume (FV), or a high-order Discontinuous Galerkin (DG) discretization which affect the arithmetic intensity and memory-bandwidth of the algebraic operations. We consider test cases of transonic turbulent flows with RANS modeling over the two-dimensional supercritical OAT15A airfoil and the three-dimensional ONERA M6 wing. Benefits in terms of robustness and convergence compared to standard GMRES solvers are obtained. Strong scalability analysis shows satisfactory results. Based on these representative problems a discussion of the recommended numerical practices is proposed.

Low-synch Gram–Schmidt with delayed reorthogonalization for Krylov solvers

The parallel strong-scaling of iterative methods is often determined by the number of global reductions at each iteration. Low-synch Gram–Schmidt algorithms are applied here to the Arnoldi algorithm to reduce the number of global reductions and therefore to improve the parallel strong-scaling of iterative solvers for nonsymmetric matrices such as the GMRES and the Krylov–Schur iterative methods. In the Arnoldi context, the QR factorization is “left-looking” and processes one column at a time. Among the methods for generating an orthogonal basis for the Arnoldi algorithm, the classical Gram–Schmidt algorithm, with reorthogonalization (CGS2) requires three global reductions per iteration. A new variant of CGS2 that requires only one reduction per iteration is presented and applied to the Arnoldi algorithm. Delayed CGS2 (DCGS2) employs the minimum number of global reductions per iteration (one) for a one-column at-a-time algorithm. The main idea behind the new algorithm is to group global reductions by rearranging the order of operations. DCGS2 must be carefully integrated into an Arnoldi expansion or a GMRES solver. Numerical stability experiments assess robustness for Krylov–Schur eigenvalue computations. Performance experiments on the ORNL Summit supercomputer then establish the superiority of DCGS2 over CGS2.

Flexible GMRES solver for interpolating MLPG analysis of heat conduction

This article evaluates the performance of flexible generalized minimum residual (FGMRES) solver in meshless local Petrov–Galerkin (MLPG) analysis of heat conduction. Interpolating MLPG formulation has been used for the meshless discretization, resulting in an asymmetric algebraic equation system. FGMRES is a variant of GMRES, which allows the preconditioner to vary from step to step. Also, the GMRES solver can be used as a preconditioner. The FGMRES method has been used with Jacobi, SOR(k), and GMRES preconditioners for MLPG analysis of heat conduction in two-dimensional and three-dimensional geometries. The performance of the FGMRES solver is compared with the restarted GMRES solver (and BiCGSTAB solver in 3-D complex geometry). Jacobi and SOR(k) preconditioners are used in the restarted GMRES and BICGSTAB solvers. Results demonstrate that the FGMRES solver performs better than the GMRES and BiCGSTAB solvers, except for two-dimensional complex geometry.

Automatic Penalty and Degree Continuation for Parallel Pre-Conditioned Mesh Curving on Virtual Geometry

We present a distributed parallel mesh curving method for virtual geometry. The main application is to generate large-scale curved meshes on complex geometry suitable for analysis with unstructured high-order methods. Accordingly, we devise the technique to generate geometrically accurate meshes composed of high-quality elements. To this end, we advocate for degree continuation on a penalty-based second-order optimizer that uses global tight tolerances to converge the distortion residuals. To reduce the method memory footprint, waiting time, and energy consumption, we combine three main ingredients. First, we propose a matrix-free GMRES solver pre-conditioned with successive over-relaxation by blocks to reduce the memory footprint three times. We also propose an adaptive penalty technique, to reduce the number of non-linear iterations. Third, we propose an indicator of the required linear solver tolerance to reduce the number of linear iterations. On thousands of cores, the method curves meshes composed of millions of quartic elements featuring highly stretched elements while matching a virtual topology.

Efficient Simulations of Large Scale Convective Heat Transfer Problems

We describe an approach for efficient solution of large scale convective heat transfer problems, formulated as coupled unsteady heat conduction and incompressible fluid flow equations. The original problem is discretized in time using classical implicit methods, while stabilized finite elements are used for space discretization. The algorithm employed for the discretization of the fluid flow problem uses Picard's iterations to solve the arising nonlinear equations. Both problems, heat transfer and Navier-Stokes quations, give rise to large sparse systems of linear equations. The systems are solved using iterative GMRES solver with suitable preconditioning. For the incompressible flow equations we employ a special preconditioner based on algebraic multigrid (AMG) technique. The paper presents algorithmic and implementation details of the solution procedure, which is suitably tuned, especially for ill conditioned systems arising from discretizations of incompressible Navier-Stokes equations. We describe parallel implementation of the solver using MPI and elements of PETSC library. The scalability of the solver is favourably compared with other methods such as direct solvers and standard GMRES method with ILU preconditioning.

Benchmarks of Cuda-Based GMRES Solver for Toeplitz and Hankel Matrices and Applications to Topology Optimization of Photonic Components

Benchmarks of Cuda-Based GMRES Solver for Toeplitz and Hankel Matrices and Applications to Topology Optimization of Photonic Components

Inexact Inverse Subspace Iteration with Preconditioning Applied to Quadratic Matrix Polynomials

Abstract An inexact variant of inverse subspace iteration is used to find a small invariant pair of a large quadratic matrix polynomial. It is shown that linear convergence is preserved provided the inner iteration is performed with increasing accuracy. A preconditioned block GMRES solver is employed as inner iteration. The preconditioner uses the strategy of “tuning” which prevents the inner iteration from increasing and therefore results in a substantial saving in costs. The accuracy of the computed invariant pair can be improved by the addition of a post-processing step involving very few iterations of Newton’s method. The effectiveness of the proposed approach is demonstrated by numerical experiments.

Weighted integral solvers for elastic scattering by open arcs in two dimensions

AbstractWe present new methodologies for the numerical solution of problems of elastic scattering by open arcs in two dimensions. The algorithms utilize weighted versions of the classical elastic integral operators associated with Dirichlet and Neumann boundary conditions, where the integral weight accounts for (and regularizes) the singularity of the integral‐equation solutions at the open‐arc endpoints. Crucially, the method also incorporates a certain “open‐arc elastic Calderón relation” introduced in this article, whose validity is demonstrated on the basis of numerical experiments, but whose rigorous mathematical proof is left for future work. (In fact, the aforementioned open‐arc elastic Calderón relation generalizes a corresponding elastic Calderón relation for closed surfaces, which is also introduced in this article, and for which a rigorous proof is included.) Using the open‐surface Calderón relation in conjunction with spectrally accurate quadrature rules and the Krylov‐subspace linear algebra solver GMRES, the proposed overall open‐arc elastic solver produces results of high accuracy in small number of iterations, for both low and high frequencies. A variety of numerical examples in this article demonstrate the accuracy and efficiency of the proposed methodology.

Parallel graph-grammar-based algorithm for the longest-edge refinement of triangular meshes and the pollution simulations in Lesser Poland area

In this paper, we propose parallel graph-grammar-based algorithm for the longest-edge refinements and the pollution simulations in Lesser Poland area. We introduce graph-grammar productions for Rivara’s longest-edged algorithm for the local refinement of unstructured triangular meshes. We utilize the hyper-graph to represent the computational mesh and the graph-grammar productions to express the longest-edge mesh refinement algorithm. The parallelism in the original Rivara’s longest edge refinement algorithm is obtained by processing different longest edge refinement paths in different three ads. Our graph-grammar-based algorithm allows for additional parallelization within a single longest-edge refinement path. The graph-grammar-based algorithm automatically guarantees the validity and conformity of the generated mesh; it prevents the generation of duplicated nodes and edges, elongated elements with Jacobians converging to zero, and removes all the hanging nodes automatically from the mesh. We test the algorithm on generating a surface mesh based on a topographic data of Lesser Poland area. The graph-grammar productions also generate the layers of prismatic three-dimensional elements on top of the triangular mesh, and they break each prismatic element into three tetrahedral elements. Next, we propose graph-grammar productions generating element matrices and right-hand-side vectors for each tetrahedral element. We utilize the Streamline Upwind Petrov–Galerkin (SUPG) stabilization for the pollution propagation simulations in Lesser Poland area. We use the advection–diffusion-reaction model, the Crank–Nicolson time integration scheme, and the graph-grammar-based interface to the GMRES solver.

Preconditioners for multilevel Toeplitz linear systems from steady-state and evolutionary advection-diffusion equations

In this paper, we study preconditioners for multilevel Toeplitz linear systems arising from discretization of steady-state and evolutionary advection-diffusion equations, in which upwind scheme and central difference scheme are employed to discretize first-order and second-order terms, respectively. For the steady-state case, the preconditioner is constructed by replacing each of the discrete advection terms with a square root of the negative of discrete Laplacian matrix and the so constructed preconditioner is diagonalizable by a sine transform. Due to its diagonalizability, the preconditioner can be applied in a two-sided way. We prove that the GMRES solver for the preconditioned linear system has a linear convergence rate independent of discretization step-sizes. The sum of the time discretization and the steady-state preconditioner constitutes the evolutionary preconditioner. A fast implementation is proposed for the evolutionary preconditioner. Moreover, for the evolutionary case, we prove that the modulus of the eigenvalues of the preconditioned matrix is lower and upper bounded by positive constants independent of discretization step-sizes. We test the proposed preconditioners with several Krylov subspace solvers on some advection-dominated advection-diffusion problems and compare their performance with other preconditioners to show its efficiency.

A constrained transport divergence-free finite element method for incompressible MHD equations

In this paper we study finite element method for three-dimensional incompressible resistive magnetohydrodynamic equations, in which the velocity, the current density, and the magnetic induction are divergence-free. It is desirable that the discrete solutions should also satisfy divergence-free conditions exactly especially for the momentum equations. Inspired by constrained transport method, we devise a new stable mixed finite element method that can achieve the goal. We also prove the well-posedness of the discrete solutions. To solve the resulting linear algebraic equations, we propose a GMRES solver with an augmented Lagrangian block preconditioner. By numerical experiments, we verify the theoretical results and demonstrate the quasi-optimality of the discrete solver with respect to the number of degrees of freedom. A comparison with other discretization using lid driven cavity is also given.

A Chebyshev-based rectangular-polar integral solver for scattering by geometries described by non-overlapping patches

This paper introduces a high-order-accurate strategy for integration of singular kernels and edge-singular integral densities that appear in the context of boundary integral equation formulations for the problem of acoustic scattering. In particular, the proposed method is designed for use in conjunction with geometry descriptions given by a set of arbitrary non-overlapping logically-quadrilateral patches—which makes the algorithm particularly well suited for computer-aided design (CAD) geometries. Fejér's first quadrature rule is incorporated in the algorithm, to provide a spectrally accurate method for evaluation of contributions from far integration regions, while highly-accurate precomputations of singular and near-singular integrals over certain “surface patches” together with two-dimensional Chebyshev transforms and suitable surface-varying “rectangular-polar” changes of variables, are used to obtain the contributions for singular and near-singular interactions. The overall integration method is then used in conjunction with the linear-algebra solver GMRES to produce solutions for sound-soft open- and closed-surface scattering obstacles, including an application to an aircraft described by means of a CAD representation. The approach is robust, fast, and highly accurate: use of a few points per wavelength suffices for the algorithm to produce far-field accuracies of a fraction of a percent, and slight increases in the discretization densities give rise to significant accuracy improvements.

A Discussion on the Solution of Linear Systems Appeared in Gradient Coil Optimization with the Finite Element Method

In the gradient coil optimization process, a dense linear system usually needs to be solved. When the mesh discretization method is used, the matrix will be very large, which makes direct methods very inefficient. The solution of the linear system is discussed in this paper for the optimization of gradient coils with the finite element method. Finally, a transversal planar gradient coil is designed and the efficiency of different methods are compared for the solution of the problem. Numerical results showed that the GMRES solver is the most efficient method.