Symmetric Positive Definite Systems Research Articles

Mixed constraint preconditioning in computational contact mechanics

Contact mechanics can be addressed numerically by finite elements using either a penalty formulation or the Lagrange multipliers. The penalty approach leads to a linearized symmetric positive definite system which can prove severely ill-conditioned, with the iterative solution to large 3D problems requiring expensive preconditioners to accelerate, or even to allow for, convergence. In the present paper a mixed constraint preconditioner (MCP) is developed and implemented for the iterative solution of contact problems based on the penalty approach. A theoretical bound for the eigenspectrum of the preconditioned matrix is provided depending on the user-specified parameters that control the preconditioner fill-in. The performance of MCP combined with the preconditioned conjugate gradient is studied in two real 3D test cases describing the geomechanics of faulted porous media and compared to a well-established preconditioner such as the incomplete Cholesky decomposition ILLT with optimal fill-in. Numerical results show that typically MCP outperforms ILLT especially in the most challenging problems where very large penalty parameters must be used. Finally, a modification aimed at simplifying the MCP implementation and increasing its robustness is developed and experimented with, showing that the modified MCP can be viewed as an efficient and viable alternative to both MCP and ILLT, especially in problems where the size of the contact surfaces is large with respect to the global domain.

Parallel MIC(0) preconditioning of 3D elliptic problems discretized by Rannacher–Turek finite elements

Novel parallel algorithms for the solution of large FEM linear systems arising from second order elliptic partial differential equations in 3D are presented. The problem is discretized by rotated trilinear nonconforming Rannacher–Turek finite elements. The resulting symmetric positive definite system of equations A x = f is solved by the preconditioned conjugate gradient algorithm. The preconditioners employed are obtained by the modified incomplete Cholesky factorization MIC(0) of two kinds of auxiliary matrices B that both are constructed as locally optimal approximations of A in the class of M -matrices. Uniform estimates for the condition number κ ( B − 1 A ) are derived. Two parallel algorithms based on the different block structures of the related matrices B are studied. The numerical tests confirm theory in that the algorithm scales as O ( N 7 / 6 ) in the matrix order N .

An efficient solver for the equations of resistive MHD with spatially-varying resistivity

We regularize the variable coefficient Helmholtz equations arising from implicit time discretizations for resistive MHD, in a way that leads to a symmetric positive-definite system uniformly in the time step. Standard centered-difference discretizations in space of the resulting PDE leads to a method that is second-order accurate, and that can be used with multigrid iteration to obtain efficient solvers.

A comparison of abstract versions of deflation, balancing and additive coarse grid correction preconditioners

AbstractIn this paper we consider various preconditioners for the conjugate gradient (CG) method to solve large linear systems of equations with symmetric positive definite system matrix. We continue the comparison between abstract versions of the deflation, balancing and additive coarse grid correction preconditioning techniques started in (SIAM J. Numer. Anal. 2004; 42:1631–1647; SIAM J. Sci. Comput. 2006; 27:1742–1759). There the deflation method is compared with the abstract additive coarse grid correction preconditioner and the abstract balancing preconditioner. Here, we close the triangle between these three methods. First of all, we show that a theoretical comparison of the condition numbers of the abstract additive coarse grid correction and the condition number of the system preconditioned by the abstract balancing preconditioner is not possible. We present a counter example, for which the condition number of the abstract additive coarse grid correction preconditioned system is below the condition number of the system preconditioned with the abstract balancing preconditioner. However, if the CG method is preconditioned by the abstract balancing preconditioner and is started with a special starting vector, the asymptotic convergence behavior of the CG method can be described by the so‐called effective condition number with respect to the starting vector. We prove that this effective condition number of the system preconditioned by the abstract balancing preconditioner is less than or equal to the condition number of the system preconditioned by the abstract additive coarse grid correction method. We also provide a short proof of the relationship between the effective condition number and the convergence of CG. Moreover, we compare the A‐norm of the errors of the iterates given by the different preconditioners and establish the orthogonal invariants of all three types of preconditioners. Copyright © 2008 John Wiley & Sons, Ltd.

Topology optimization of shell structures using adaptive inner-front (AIF) level set method

A new topology optimization using adaptive inner-front level set method is presented. In the conventional level set-based topology optimization, the optimum topology strongly depends on the initial level set due to the incapability of inner-front creation during the optimization process. In the present work, in this regard, an algorithm for inner-front creation is proposed in which the sizes, the positions, and the number of new inner-fronts during the optimization process can be globally and consistently identified. In the algorithm, the criterion of inner-front creation for compliance minimization problems of linear elastic structures is chosen as the strain energy density along with volumetric constraint. To facilitate the inner-front creation process, the inner-front creation map is constructed and used to define new level set function. In the implementation of inner-front creation algorithm, to suppress the numerical oscillation of solutions due to the sharp edges in the level set function, domain regularization is carried out by solving the edge smoothing partial differential equation (smoothing PDE). To update the level set function during the optimization process, the least-squares finite element method (LSFEM) is adopted. Through the LSFEM, a symmetric positive definite system matrix is constructed, and non-diffused and non-oscillatory solution for the hyperbolic PDE such as level set equation can be obtained. As applications, three-dimensional topology optimization of shell structures is treated. From the numerical examples, it is shown that the present method brings in much needed flexibility in topologies during the level set-based topology optimization process.

Recursive Krylov‐based multigrid cycles

AbstractWe consider multigrid (MG) cycles based on the recursive use of a two‐grid method, in which the coarse‐grid system is solved by µ⩾1 steps of a Krylov subspace iterative method. The approach is further extended by allowing such inner iterations only at the levels of given multiplicity, whereas V‐cycle formulation is used at all other levels. For symmetric positive definite systems and symmetric MG schemes, we consider a flexible (or generalized) conjugate gradient method as Krylov subspace solver for both inner and outer iterations. Then, based on some algebraic (block matrix) properties of the V‐cycle MG viewed as a preconditioner, we show that the method can have optimal convergence properties if µ is chosen to be sufficiently large. We also formulate conditions that guarantee both, optimal complexity and convergence, bounded independently of the number of levels. Our analysis shows that the method is, at least, as effective as the standard W‐cycle, whereas numerical results illustrate that it can be much faster than the latter, and actually more robust than predicted by the theory. Copyright © 2007 John Wiley & Sons, Ltd.

Boundary conditions for a divergence free velocity–pressure formulation of the Navier–Stokes equations

New sets of boundary conditions for the velocity–pressure formulation of the incompressible Navier–Stokes equations are derived. The boundary conditions have the same form on both inflow and outflow boundaries and lead to a divergence free solution. Moreover, the specific form of the boundary condition makes it possible derive a symmetric positive definite equation system for the internal pressure. Numerical experiments support the theoretical conclusions.

Analysis of block matrix preconditioners for elliptic optimal control problems

AbstractIn this paper, we describe and analyse several block matrix iterative algorithms for solving a saddle point linear system arising from the discretization of a linear‐quadratic elliptic control problem with Neumann boundary conditions. To ensure that the problem is well posed, a regularization term with a parameter α is included. The first algorithm reduces the saddle point system to a symmetric positive definite Schur complement system for the control variable and employs conjugate gradient (CG) acceleration, however, double iteration is required (except in special cases). A preconditioner yielding a rate of convergence independent of the mesh size h is described for Ω ⊂ R2 or R3, and a preconditioner independent of h and α when Ω ⊂ R2. Next, two algorithms avoiding double iteration are described using an augmented Lagrangian formulation. One of these algorithms solves the augmented saddle point system employing MINRES acceleration, while the other solves a symmetric positive definite reformulation of the augmented saddle point system employing CG acceleration. For both algorithms, a symmetric positive definite preconditioner is described yielding a rate of convergence independent of h. In addition to the above algorithms, two heuristic algorithms are described, one a projected CG algorithm, and the other an indefinite block matrix preconditioner employing GMRES acceleration. Rigorous convergence results, however, are not known for the heuristic algorithms. Copyright © 2007 John Wiley & Sons, Ltd.

Preconditioners for the conjugate gradient algorithm using Gram–Schmidt and least squares methods

This paper is devoted to the study of some preconditioners for the conjugate gradient algorithm used to solve large sparse linear and symmetric positive definite systems. The construction of a preconditioner based on the Gram–Schmidt orthogonalization process and the least squares method is presented. Some results on the condition number of the preconditioned system are provided. Finally, numerical comparisons are given for different preconditioners.

Locally Conservative Fluxes for the Continuous Galerkin Method

The standard continuous Galerkin (CG) finite element method for second order elliptic problems suffers from its inability to provide conservative flux approximations, a much needed quantity in many applications. We show how to overcome this shortcoming by using a two-step postprocessing. The first step is the computation of a numerical flux trace defined on element interfaces and is motivated by the structure of the numerical traces of discontinuous Galerkin methods. This computation is nonlocal in that it requires the solution of a symmetric positive definite system, but the system is well conditioned independently of mesh size, so it can be solved at asymptotically optimal cost. The second step is a local element-by-element postprocessing of the CG solution incorporating the result of the first step. This leads to a conservative flux approximation with continuous normal components. This postprocessing applies for the CG method in its standard form or for a hybridized version of it. We present the hybridized version since it allows easy handling of variable-degree polynomials and hanging nodes. Furthermore, we provide an a priori analysis of the error in the postprocessed flux approximation and display numerical evidence suggesting that the approximation is competitive with the approximation provided by the Raviart-Thomas mixed method of corresponding degree.

A numerical experimental study of inverse preconditioning for the parallel iterative solution to 3D finite element flow equations

Integration of the subsurface flow equation by finite elements (FE) in space and finite differences (FD) in time requires the repeated solution to sparse symmetric positive definite systems of linear equations. Iterative techniques based on preconditioned conjugate gradients (PCG) are one of the most attractive tool to solve the problem on sequential computers. A present challenge is to make PCG attractive in a parallel computing environment as well. To this aim a key factor is the development of an efficient parallel preconditioner. FSAI (factorized sparse approximate inverse) and enlarged FSAI relying on the approximate inverse of the coefficient matrix appears to be a most promising parallel preconditioner. In the present paper PCG using FSAI, diagonal and pARMS (parallel algebraic recursive multilevel solvers) preconditioners is implemented on the IBM SP4/512 and CLX/768 supercomputers with up to 32 processors to solve underground flow problems of a large size. The results show that FSAI may allow for a parallel relative efficiency E p * larger than 50% on the largest problems with p = 32 processors. Moreover, FSAI turns out to be significantly less expensive and more robust than pARMS. Finally, it is shown that E p * for p in the upper range may be much improved if PCG–FSAI is implemented on CLX.

Analysis of a novel preconditioner for a class of p‐level lower rank extracted systems

AbstractThis paper proposes and studies the performance of a preconditioner suitable for solving a class of symmetric positive definite systems, Âx=b, which we call p‐level lower rank extracted systems (p‐level LRES), by the preconditioned conjugate gradient method. The study of these systems is motivated by the numerical approximation of integral equations with convolution kernels defined on arbitrary p‐dimensional domains. This is in contrast to p‐level Toeplitz systems which only apply to rectangular domains. The coefficient matrix, Â, is a principal submatrix of a p‐level Toeplitz matrix, A, and the preconditioner for the preconditioned conjugate gradient algorithm is provided in terms of the inverse of a p‐level circulant matrix constructed from the elements of A. The preconditioner is shown to yield clustering in the spectrum of the preconditioned matrix which leads to a substantial reduction in the computational cost of solving LRE systems. Copyright © 2006 John Wiley & Sons, Ltd.

A Parallel Distributed Solver for Large Dense Symmetric Systems: Applications to Geodesy and Electromagnetism Problems

In this paper we describe the parallel distributed implementation of a linear solver for large-scale applications involving real symmetric positive definite or complex symmetric non-Hermitian dense systems. The advantage of this routine is that it performs a Cholesky factorization by requiring half the storage needed by the standard parallel libraries ScaLAPACK and PLAPACK. Our solver uses a Jvariant Cholesky algorithm and a one-dimensional blockcyclic column data distribution but gives similar Gigaflops performance when applied to problems that can be solved on moderately parallel computers with up to 32 processors. Experiments and performance comparisons with ScaLAPACK and PLAPACK on our target applications are presented. These applications arise from the Earth's gravity field recovery and computational electromagnetics.

A bound constrained conjugate gradient solution method as applied to crystallographic refinement problems

Small modifications to the conjugate gradient method for solving symmetric positive definite systems have resulted in an increase in performance over LU decomposition by a factor of around 84 for solving a dense system of 1325 unknowns. Performance is further increased in the case of applying upper- and lower-bound parameter constraints. For structure solution employing simulated annealing and the Newton–Raphson method of non-linear least squares, the overall performance gain can be a factor of four, depending on the applied constraints. In addition, the new algorithm with bounding constraints often seeks out lower minima than would otherwise be attainable without constraints. The behaviour of the new algorithm has been tested against the crystallographic problems of Pawley refinement, rigid-body and general crystal structure refinement.

A parallel least-squares spectral element solver for incompressible flow problems on unstructured grids

The parallelisation of the least-squares spectral element formulation of the Stokes problem is discussed for incompressible flow problems on unstructured grids. The method leads to a large symmetric positive definite algebraic system, that is solved iteratively by the conjugate gradient method. To improve the convergence rate, both Jacobi and Additive Schwarz preconditioners are applied. Numerical simulations have been performed to validate the scalability of the different parts of the proposed method. The experiments entailed simulating several large-scale incompressible flows on a Cray T3E and on an SGI Origin 3800.

Determination of a good value of the time step and preconditioned Krylov subspace methods for the Navier-Stokes equations

The main purpose of this paper is to develop stable versions of some Krylov subspace methods for solving the linear systems of equations Ax = b which arise in the difference solution of 2-D nonstationary Navier-Stokes equations using implicit scheme and to determine a good value of the time step. Our algorithms are based on the conjugate-gradient method with a suitable preconditioner for solving the symmetric positive definite system and preconditioned GMRES, Orthomin(K), QMR methods for solving the nonsymmetric and (in)definite system. The performance of these methods is compared. In addition, we show that by using the condition number of the first nonsymmetric coefficient matrix, it is possible to determine a good value of the time step.

Conjugate gradient method for fuzzy symmetric positive definite system of linear equations

In this paper the conjugate gradient method, for solving fuzzy symmetric positive definite system of linear equation is considered. The method in detail is discussed and followed by convergence theorem and illustrated by solving some numerical examples.

Analysis of a circulant based preconditioner for a class of lower rank extracted systems

This paper proposes and studies the performance of a preconditioner suitable for solving a class of symmetric positive definite systems, Apx=b, which we call lower rank extracted systems (LRES), by the preconditioned conjugate gradient method. These systems correspond to integral equations with convolution kernels defined on a union of many line segments in contrast to only one line segment in the case of Toeplitz systems. The p × p matrix, Ap, is shown to be a principal submatrix of a larger N × N Toeplitz matrix, AN. The preconditioner is provided in terms of the inverse of a 2N × 2N circulant matrix constructed from the elements of AN. The preconditioner is shown to yield clustering in the spectrum of the preconditioned matrix similar to the clustering results for iterative algorithms used to solve Toeplitz systems. The analysis also demonstrates that the computational expense to solve LRE systems is reduced to O(N log N). Copyright © 2004 John Wiley & Sons, Ltd.

Finite element method for time dependent scattering: nonreflecting boundary condition, adaptivity, and energy decay

An adaptive finite element method is developed for acoustic wave propagation in unbounded media. The efficiency and high accuracy of the method are achieved by combining an exact nonreflecting boundary condition [SIAM J. Appl. Math. 55 (1995) 280; J. Comput. Phys. 127 (1996) 52] with space–time adaptivity [East–West J. Numer. Math. 7(4) (1999) 263]. Hence the computational effort is concentrated where needed, while the artificial boundary can be brought as close as desired to the scatterer. Both features combined yield high accuracy and keep the number of unknowns to a minimum. An energy inequality is derived for the initial-boundary value problem at the continuous level. Together with an implicit second order time discretization it guarantees unconditional stability of the semi-discrete system. The resulting fully discrete linear system that needs to be solved every time step is unsymmetric but can be transformed into an equivalent sequence of small nonsymmetric and large symmetric positive definite systems, which are efficiently solved by conjugate gradient methods. Numerical examples illustrate the high accuracy of the method, in particular in the presence of complex geometry.

W Circulant Boundary Condition

AbstractWe use the preconditioned conjugate gradient (PCG) method to analyze the symmetric positive‐definite Toeplitz systems. Then we present a new embedding constructing preconditioner way and prove that many other ways of constructing preconditioners are generally the special cases of this method. We propose the w‐circulant boundary condition and make comparison with both the ordinary circulant boundary condition and helix one. Furthermore, we obtain a new type of boundary condition: the mixed boundary condition.