Solution Of Large Linear Systems Research Articles

Caveats of three direct linear solvers for finite element analyses

AbstractThe solution of large linear systems of equations with sparse matrices is a critical component of finite element analyses. Three linear solvers are investigated here: MUMPS, UMFPACK, and Intel DSS (PARDISO). Often, these solvers are employed as “black boxes.” However, some caveats in their implementation must be observed. For instance, the solvers may yield incorrect results or perform extremely poorly in a multithread environment. These issues are demonstrated, and suggestions to fix them are provided. Some performance benchmarks are also presented with a focus on the multithreaded behavior.

On iterative methods based on Sherman-Morrison-Woodbury splitting

We consider a new splitting based on the Sherman-Morrison-Woodbury formula, which is particularly effective with iterative methods for the numerical solution of large linear systems and especially for systems involving matrices that are perturbations of circulant or block circulant matrices. Such matrices typically arise in the discretization of differential equations using finite element or finite difference methods. We prove the convergence of the new iteration without making any assumptions regarding the symmetry or diagonal-dominance of the matrix, which are limiting factors for most classical iterative methods.To illustrate the efficacy of the new iteration we present various applications. These include extensions of the new iteration to block matrices that arise in certain saddle point problems as well as two-dimensional finite difference discretizations. The new method exhibits fast convergence in all of the test cases we used. It has minimal storage requirements, straightforward implementation and compatibility with nearly circulant matrices via the Fast Fourier Transform. Remarkably, the new method was tested against very large matrices demonstrating extremely fast convergence. For these reasons it can be a valuable tool for the solution of various finite elements and finite differences discretizations of differential equations.

A new numerical mesoscopic scale one-domain approach solver for free fluid/porous medium interaction

A new numerical continuum one-domain approach (ODA) solver is presented for the simulation of the transfer processes between a free fluid and a porous medium. The solver is developed in the mesoscopic scale framework, where a continuous variation of the physical parameters of the porous medium (e.g., porosity and permeability) is assumed. The Navier–Stokes–Brinkman equations are solved along with the continuity equation, under the hypothesis of incompressible fluid. The porous medium is assumed to be fully saturated and can potentially be anisotropic. The domain is discretized with unstructured meshes allowing local refinements. A fractional time step procedure is applied, where one predictor and two corrector steps are solved within each time iteration. The predictor step is solved in the framework of a marching in space and time procedure, with some important numerical advantages. The two corrector steps require the solution of large linear systems, whose matrices are sparse, symmetric and positive definite, with M-matrix property over Delaunay-meshes. A fast and efficient solution is obtained using a preconditioned conjugate gradient method. The discretization adopted for the two corrector steps can be regarded as a Two-Point-Flux-Approximation (TPFA) scheme, which, unlike the standard TPFA schemes, does not require the grid mesh to be K-orthogonal, (with K the anisotropy tensor). As demonstrated with the provided test cases, the proposed scheme correctly retains the anisotropy effects within the porous medium. Furthermore, it overcomes the restrictions of existing mesoscopic scale one-domain approaches proposed in the literature.

Stein-based preconditioners for weak-constraint 4D-var

Algorithms for data assimilation try to predict the most likely state of a dynamical system by combining information from observations and prior models. Variational approaches, such as the weak-constraint four-dimensional variational data assimilation formulation considered in this problem, can ultimately be interpreted as a minimization problem. One of the main challenges of such a formulation is the solution of large linear systems of equations which arise within the inner linear step of the adopted nonlinear solver. Depending on the adopted approach, these linear algebraic problems amount to either a saddle point linear system or a symmetric positive definite (SPD) one. Both formulations can be solved by means of a Krylov method, like GMRES or CG, that needs to be preconditioned to ensure fast convergence in terms of the number of iterations. In this paper we illustrate novel, efficient preconditioning operators which involve the solution of certain Stein matrix equations. In addition to achieving better computational performance, the latter machinery allows us to derive tighter bounds for the eigenvalue distribution of the preconditioned linear system for certain problem settings. A panel of diverse numerical results displays the effectiveness of the proposed methodology compared to current state-of-the-art approaches.

Quadratic programming based partitioning for Block Cimmino with correct value representation

The block Cimmino method is successfully used for the parallel solution of large linear systems of equations due to its amenability to parallel processing. Since the convergence rate of block Cimmino depends on the orthogonality between the row blocks, advanced partitioning methods are used for faster convergence. In this work, we propose a new partitioning method that is superior to the state-of-the-art partitioning method, GRIP, in several ways. Firstly, our proposed method exploits the Mongoose partitioning library which can outperform the state-of-the-art methods by combining the advantages of classical combinatoric methods and continuous quadratic programming formulations. Secondly, the proposed method works on the numerical values in a floating-point format directly without converting them to integer format as in GRIP. This brings an additional advantage of obtaining higher quality partitionings via better representation of numerical values. Furthermore, the preprocessing time is also improved since there is no overhead in converting numerical values to integer format. Finally, we extend the Mongoose library, which originally partitions graphs into only two parts, by using the recursive bisection paradigm to partition graphs into more than two parts. Extensive experiments conducted on both shared and distributed memory architectures demonstrate the effectiveness of the proposed method for solving different types of real-world problems.

Bridging finite element and deep learning: High-resolution stress distribution prediction in structural components

Finite-element analysis (FEA) for structures has been broadly used to conduct stress analysis of various civil and mechanical engineering structures. Conventional methods, such as FEA, provide high fidelity results but require the solution of large linear systems that can be computationally intensive. Instead, Deep Learning (DL) techniques can generate results significantly faster than conventional run-time analysis. This can prove extremely valuable in real-time structural assessment applications. Our proposed method uses deep neural networks in the form of convolutional neural networks (CNN) to bypass the FEA and predict high-resolution stress distributions on loaded steel plates with variable loading and boundary conditions. The CNN was designed and trained to use the geometry, boundary conditions, and load as input to predict the stress contours. The proposed technique’s performance was compared to finite-element simulations using a partial differential equation (PDE) solver. The trained DL model can predict the stress distributions with a mean absolute error of 0.9% and an absolute peak error of 0.46% for the von Mises stress distribution. This study shows the feasibility and potential of using DL techniques to bypass FEA for stress analysis applications.

A method for computing a few eigenpairs of large generalized eigenvalue problems

Many methods for the computation of selected eigenpairs of generalized eigenproblems for matrix pairs use a shift-and-invert technique. When applied to large-scale problems, this requires the solution of large linear systems of equations. This paper proposes an application of an Arnoldi method described in Voss (2004) [21] to the computation of a few extreme eigenpairs of a matrix pair. An advantage of this approach, when compared to methods that use the shift-and-invert technique, is that no large systems of equations have to be solved. We compare this approach to using a technique for simultaneously reducing a pair of large matrices to a pair of small matrices by a generalized Arnoldi process described in Li and Ye (2003) [12] and Hoffnung et al. (2006) [6]. The latter technique does not require the solution of large linear systems of equations either. Computed examples show the proposed method to yield approximations of the desired eigenpairs of higher accuracy when using about the same amount of computer storage space.

An assessment of solvers for algebraically stabilized discretizations of convection–diffusion–reaction equations

Abstract We consider flux-corrected finite element discretizations of 3D convection-dominated transport problems and assess the computational efficiency of algorithms based on such approximations. The methods under investigation include flux-corrected transport schemes and monolithic limiters. We discretize in space using a continuous Galerkin method and ℙ1 or ℚ1 finite elements. Time integration is performed using the Crank–Nicolson method or an explicit strong stability preserving Runge–Kutta method. Nonlinear systems are solved using a fixed-point iteration method, which requires solution of large linear systems at each iteration or time step. The great variety of options in the choice of discretization methods and solver components calls for a dedicated comparative study of existing approaches. To perform such a study, we define new 3D test problems for time dependent and stationary convection–diffusion–reaction equations. The results of our numerical experiments illustrate how the limiting technique, time discretization and solver impact on the overall performance.

A new randomized vector algorithm for iterative solution of large linear systems

In this letter we suggest a new randomized scalable stochastic-matrix-based algorithms for calculation of large matrix iterations. Special focus is on positive or irreducible nonnegative class of matrices. As an application, a new randomized vector algorithm for iterative solution of large linear systems of algebraic equations governed by M-matrices is constructed. The idea behind these stochastic methods is in a randomized vector representation of matrix iterations. The iterations are performed by sampling random columns only, thus avoiding not only matrix but also matrix vector multiplications. As a result, the algorithm is highly efficient for solving linear equations of high dimension, its computational cost depends linearly on the dimension. Extensions of the suggested randomized iteration method to general classes of matrices are also discussed.

A new solver for incompressible non-isothermal flows in natural and mixed convection over unstructured grids

In the present paper we propose a new numerical methodology for the solution of 2D non-isothermal incompressible flows for natural and mixed convection in irregular geometries. The governing equations are the Incompressible Navier-Stokes equations and the energy conservation equation. Fluid velocity and temperature are coupled in the buoyancy term of the momentum equations according to the Oberbeck-Boussinesq approximation. The governing equations are discretized over unstructured triangular meshes satisfying the Delaunay property. Thanks to the Oberbeck-Boussinesq hypothesis, the flow and energy problems are solved in an uncoupled way, and two fractional time step procedures are sequentially applied to solve each problem. The prediction steps of both procedures are solved applying a Marching in Space and Time (MAST) numerical Eulerian scheme, which explicitly handles the non-linear terms in the momentum equations of the fluid problem, and allows numerical stability for Courant numbers greater than one. An analytical solution is applied for the prediction thermal problem. The correction steps of the two fractional time step procedures involve the solution of large linear systems, whose matrices are sparse and symmetric and have the “M” property if the mesh satisfies the Delaunay condition. This allows a well performing condition number. The matrix coefficients are constant in time, so that they are calculated and factorized only once, before the simulations loop starts, saving a lot of computational time. We present four numerical applications. For the first test, an analytical solution is available, and this makes it possible to analyze the spatial convergence order of the numerical solver. In the other applications, we investigate the capability of the proposed algorithm to handle irregular geometries, and we compare the computed results with experimental data and the outputs of literature models. A study of the required computational costs is also presented.

Second Order Fully Semi-Lagrangian Discretizations of Advection-Diffusion-Reaction Systems

We propose a second order, fully semi-Lagrangian method for the numerical solution of systems of advection-diffusion-reaction equations, which is based on a semi-Lagrangian approach to approximate in time both the advective and the diffusive terms. The proposed method allows to use large time steps, while avoiding the solution of large linear systems, which would be required by implicit time discretization techniques. Standard interpolation procedures are used for the space discretization on structured and unstructured meshes. A novel extrapolation technique is proposed to enforce second-order accurate Dirichlet boundary conditions. We include a theoretical analysis of the scheme, along with numerical experiments which demonstrate the effectiveness of the proposed approach and its superior efficiency with respect to more conventional explicit and implicit time discretizations.

Fully implicit multiple graphics processing units’ schemes for hypersonic flows with lower upper symmetric Gauss-Seidel preconditioner on unstructured non-orthogonal grids

The paper discusses the simulation of the hypersonic flow problems on unstructured grids using fully implicit scheme and the analysis of the orthogonal corrections on the quality of heat fluxes on the surfaces of the bodies subject to hypersonic flow. The governing equations are the Navier-Stokes equations for viscous calorically perfect gas. Second order finite volume discretization is applied to the problem on arbitrary unstructured grids with implicit temporal treatment. The problem is solved using the Newton-Raphson method. Such methods require the solution of large linear systems with iterative solvers and require a preconditioning operator to converge. In this paper a well known Lower-Upper Symmetric Gauss Seidel (LU-SGS) preconditioner is applied. The new method is based on the tricky reordering for the factored Jacobi matrix that allows one to execute block triangular matrix solvers on GPUs in parallel without the loss of algebraic properties of the original non-factored operator. Good convergence properties are demonstrated for large scale problems of external aerodynamics with implicit Courant number around 1000 - 10000 for flows with Mach number around 13-25. The influence of the non-orthogonal corretions on the fluxes for highly scewed anisotropic grids is demonstrated at Mach numbers around 20. It is shown that the non-orthogonality corrections are of the second infinitesimal order compared to the effects of the shock fix influence on the heat fluxes on the body.

Mixed Kirchhoff stress–displacement–pressure formulations for incompressible hyperelasticity

The numerical approximation of hyperelasticity must address nonlinear constitutive laws, geometric nonlinearities associated with large strains and deformations, the imposition of the incompressibility of the solid, and the solution of large linear systems arising from the discretisation of 3D problems in complex geometries. We adapt the three-field formulation for nearly incompressible hyperelasticity introduced in Chavan et al. (2007) to the fully incompressible case. The mixed formulation is of Hu–Washizu type and it differs from other approaches in that we use the Kirchhoff stress, displacement, and pressure as principal unknowns. We also discuss the solvability of the linearised problem restricted to neo-Hookean materials, illustrating the interplay between the coupling blocks. We construct a family of mixed finite element schemes (with different polynomial degrees) for simplicial meshes and verify its error decay through computational tests. We also propose a new augmented Lagrangian preconditioner that improves convergence properties of iterative solvers. The numerical performance of the family of mixed methods is assessed with benchmark solutions, and the applicability of the formulation is further tested in a model of cardiac biomechanics using orthotropic strain energy densities. The proposed methods are advantageous in terms of physical fidelity (as the Kirchhoff stress can be approximated with arbitrary accuracy and no locking is observed) and convergence (the discretisation and the preconditioners are robust and computationally efficient, and they compare favourably at least with respect to classical displacement–pressure schemes).

Randomized Kaczmarz with averaging

The randomized Kaczmarz (RK) method is an iterative method for approximating the least-squares solution of large linear systems of equations. The standard RK method uses sequential updates, making parallel computation difficult. Here, we study a parallel version of RK where a weighted average of independent updates is used. We analyze the convergence of RK with averaging and demonstrate its performance empirically. We show that as the number of threads increases, the rate of convergence improves and the convergence horizon for inconsistent systems decreases.

A fast matrix-free elasto-plastic solver for predicting residual stresses in additive manufacturing

Process planning for additive manufacturing (AM) today relies heavily on multi-physics, multi-scale simulation. The focus of this paper is on one aspect of AM simulation, namely, part-level elasto-plastic simulation for residual stress and distortion predictions. This is one of the crucial steps in AM process optimization, but is computationally expensive, often requiring the use of large computer clusters.The primary bottleneck in elasto-plastic simulation is the repeated solution of large linear systems of equations. While there is a wide range of linear solvers, most cannot exploit the unique structured nature of the mesh underlying AM simulation. Here, we revisit a specific matrix-free solver, namely rigid-body deflated solver that has been successfully deployed for solving large linear elastic problems in such scenarios. The salient feature of this solver is that the stiffness matrix is never assembled, thereby reducing the memory requirements significantly, leading to large computational gains.The objective of this paper is to extend the above solver to elasto-plasticity by efficiently updating the element tangent stiffness matrices, and the corresponding deflation matrix. The performance of the proposed method is evaluated on a benchmark problem using multi-core CPU and GPU architectures, and compared against ANSYS. Then, part-level residual stresses and distortion are predicted using the proposed solver. The present work is restricted to associative plasticity with von-Mises yield criteria, but can be extended to other plasticity models.

A Scalable Multigrid Reduction Framework for Multiphase Poromechanics of Heterogeneous Media

Simulation of multiphase poromechanics involves solving a multi-physics problem in which multiphase flow and transport are tightly coupled with the porous medium deformation. To capture this dynamic interplay, fully implicit methods, also known as monolithic approaches, are usually preferred. The main bottleneck of a monolithic approach is that it requires solution of large linear systems that result from the discretization and linearization of the governing balance equations. Because such systems are non-symmetric, indefinite, and highly ill-conditioned, preconditioning is critical for fast convergence. Recently, most efforts in designing efficient preconditioners for multiphase poromechanics have been dominated by physics-based strategies. Current state-of-the-art "black-box" solvers such as algebraic multigrid (AMG) are ineffective because they cannot effectively capture the strong coupling between the mechanics and the flow sub-problems, as well as the coupling inherent in the multiphase flow and transport process. In this work, we develop an algebraic framework based on multigrid reduction (MGR) that is suited for tightly coupled systems of PDEs. Using this framework, the decoupling between the equations is done algebraically through defining appropriate interpolation and restriction operators. One can then employ existing solvers for each of the decoupled blocks or design a new solver based on knowledge of the physics. We demonstrate the applicability of our framework when used as a "black-box" solver for multiphase poromechanics. We show that the framework is flexible to accommodate a wide range of scenarios, as well as efficient and scalable for large problems.

A Bayesian Conjugate Gradient Method (with Discussion)

A fundamental task in numerical computation is the solution of large linear systems. The conjugate gradient method is an iterative method which offers rapid convergence to the solution, particularly when an effective preconditioner is employed. However, for more challenging systems a substantial error can be present even after many iterations have been performed. The estimates obtained in this case are of little value unless further information can be provided about, for example, the magnitude of the error. In this paper we propose a novel statistical model for this error, set in a Bayesian framework. Our approach is a strict generalisation of the conjugate gradient method, which is recovered as the posterior mean for a particular choice of prior. The estimates obtained are analysed with Krylov subspace methods and a contraction result for the posterior is presented. The method is then analysed in a simulation study as well as being applied to a challenging problem in medical imaging.

Arnoldi decomposition, GMRES, and preconditioning for linear discrete ill-posed problems

GMRES is one of the most popular iterative methods for the solution of large linear systems of equations that arise from the discretization of linear well-posed problems, such as boundary value problems for elliptic partial differential equations. The method is also applied to the iterative solution of linear systems of equations that are obtained by discretizing linear ill-posed problems, such as many inverse problems. However, GMRES does not always perform well when applied to the latter kind of problems. This paper seeks to shed some light on reasons for the poor performance of GMRES in certain situations, and discusses some remedies based on specific kinds of preconditioning. The standard implementation of GMRES is based on the Arnoldi process, which also can be used to define a solution subspace for Tikhonov or TSVD regularization, giving rise to the Arnoldi–Tikhonov and Arnoldi-TSVD methods, respectively. The performance of the GMRES, the Arnoldi–Tikhonov, and the Arnoldi-TSVD methods is discussed. Numerical examples illustrate properties of these methods.

Weighted and deflated global GMRES algorithms for solving large Sylvester matrix equations

The solution of a large-scale Sylvester matrix equation plays an important role in control and large scientific computations. In this paper, we are interested in the large Sylvester matrix equation with large dimensionA and small dimension B, and a popular approach is to use the global Krylov subspace method. In this paper, we propose three new algorithms for this problem. We first consider the global GMRES algorithm with weighting strategy, which can be viewed as a precondition method. We present three new schemes to update the weighting matrix during iterations. Due to the growth of memory requirements and computational cost, it is necessary to restart the algorithm effectively. The deflation strategy is efficient for the solution of large linear systems and large eigenvalue problems; to the best of our knowledge, little work is done on applying deflation to the (weighted) global GMRES algorithm for large Sylvester matrix equations. We then consider how to combine the weighting strategy with deflated restarting, and propose a weighted global GMRES algorithm with deflation for solving large Sylvester matrix equations. In particular, we are interested in the global GMRES algorithm with deflation, which can be viewed as a special case when the weighted matrix is chosen as the identity. Theoretical analysis is given to show rationality of the new algorithms. Numerical experiments illustrate the numerical behavior of the proposed algorithms.

A greedy algorithm for nonlinear inverse problems with an application to nonlinear inverse gravimetry

Based on the Regularized Functional Matching Pursuit (RFMP) algorithm for linear inverse problems, we present an analogous iterative greedy algorithm for nonlinear inverse problems, called RFMP_NL. In comparison to established methods for nonlinear inverse problems, the algorithm is able to combine very diverse types of basis functions, for example, localized and global functions. This is important, in particular, in geoscientific applications, where global structures have to be distinguished from local anomalies. Furthermore, in contrast to other methods, the algorithm does not require the solution of large linear systems. We apply the RFMP_NL to the nonlinear inverse problem of gravimetry, where gravitational data are inverted for the shape of the surface or inner layer boundaries of planetary bodies. This inverse problem is described by a nonlinear integral operator, for which we additionally provide the Frechet derivative. Finally, we present two synthetic numerical examples to show that it is beneficial to apply the presented method to inverse gravimetric problems.