Accelerating iterative solution methods using reduced‐order models as solution predictors

We propose the use of reduced‐order models to accelerate the solution of systems of equations using iterative solvers in time stepping schemes for large‐scale numerical simulation. The acceleration is achieved by determining an improved initial guess for the iterative process based on information in the solution vectors from previous time steps. The algorithm basically consists of two projection steps: (1) projecting the governing equations onto a subspace spanned by a low number of global empirical basis functions extracted from previous time step solutions, and (2) solving the governing equations in this reduced space and projecting the solution back on the original, high dimensional one. We applied the algorithm to numerical models for simulation of two‐phase flow through heterogeneous porous media. In particular we considered implicit‐pressure explicit‐saturation (IMPES) schemes and investigated the scope to accelerate the iterative solution of the pressure equation, which is by far the most time‐consuming part of any IMPES scheme. We achieved a substantial reduction in the number of iterations and an associated acceleration of the solution. Our largest test problem involved 93 500 variables, in which case we obtained a maximum reduction in computing time of 67%. The method is particularly attractive for problems with time‐varying parameters or source terms. Copyright © 2006 John Wiley & Sons, Ltd.     

PCG methods in transient FE analysis. Part I: First order problems

The performances of some PCG (preconditioned conjugate gradient) algorithms are evaluated in the solution of first order time dependent parabolic partial differential equations, such as the heat conduction equation, which have been spatially discretized using finite elements. ‘Consistent mass’ discretizations are preferred by the authors to ‘lumped mass’ ones and various preconditioners are then compared—diagonal, incomplete Choleski and EBE (‘element-by-element’). Recommendations are made and implications for parallel computation outlined.     

Solving three-dimensional layout optimization problems using fixed scale wavelets

 The layout optimization problem in three- dimensional elasticity is solved with a meshless, wavelet-based solution scheme. A fictitious domain approach is used to embed the design domain into a simple regular domain. The material distribution and displacement field are discretized over the fictitious domain using fixed-scale, shift-invariant wavelet expansions. A discrete form of the elasticity problem is solved using a wavelet-Galerkin technique during each iteration of the layout optimization sequence. Approximate solutions are found with an efficient preconditioned conjugate gradient (PCG) solver using non-diagonal preconditioners which lead to convergence rates that are insensitive to the level of resolution. The convergence and memory management properties of the PCG algorithm make the analysis of large-scale problems possible. Several wavelet-based layout optimization examples are included.     

Evaluating the LCD algorithm for solving linear systems of equations arising from implicit SUPG formulation of compressible flows

In this work we evaluate the performance of the left conjugate direction method recently introduced by Yuan, Golub, Plemmons and Cecílio for the solution of non‐symmetric systems of linear equations arising from the implicit semi‐discrete SUPG finite element formulation of advective–diffusive and inviscid compressible flows. We extend the original algorithm to accommodate restarts and typical element‐by‐element preconditioners. We also show how to select the first left conjugate vector to start LCD. Several problems are solved, accessing performance parameters such as number of iterations, memory requirements and CPU times, and results are compared with other algorithms, such as GMRES, TFQMR and Bi‐CGSTAB. Copyright © 2004 John Wiley & Sons, Ltd.     

Galerkin reduced‐order modeling scheme for time‐dependent randomly parametrized linear partial differential equations

In this paper, we consider the problem of constructing reduced‐order models of a class of time‐dependent randomly parametrized linear partial differential equations. Our objective is to efficiently construct a reduced basis approximation of the solution as a function of the spatial coordinates, parameter space, and time. The proposed approach involves decomposing the solution in terms of undetermined spatial and parametrized temporal basis functions. The unknown basis functions in the decomposition are estimated using an alternating iterative Galerkin projection scheme. Numerical studies on the time‐dependent randomly parametrized diffusion equation are presented to demonstrate that the proposed approach provides good accuracy at significantly lower computational cost compared with polynomial chaos‐based Galerkin projection schemes. Comparison studies are also made against Nouy's generalized spectral decomposition scheme to demonstrate that the proposed approach provides a number of computational advantages. Copyright © 2012 John Wiley & Sons, Ltd.     

Dimensional model reduction in non‐linear finite element dynamics of solids and structures

A general approach to the dimensional reduction of non‐linear finite element models of solid dynamics is presented. For the Newmark implicit time‐discretization, the computationally most expensive phase is the repeated solution of the system of linear equations for displacement increments. To deal with this, it is shown how the problem can be formulated in an approximation (Ritz) basis of much smaller dimension. Similarly, the explicit Newmark algorithm can be also written in a reduced‐dimension basis, and the computation time savings in that case follow from an increase in the stable time step length. In addition, the empirical eigenvectors are proposed as the basis in which to expand the incremental problem. This basis achieves approximation optimality by using computational data for the response of the full model in time to construct a reduced basis which reproduces the full system in a statistical sense. Because of this ‘global’ time viewpoint, the basis need not be updated as with reduced bases computed from a linearization of the full finite element model. If the dynamics of a finite element model is expressed in terms of a small number of basis vectors, the asymptotic cost of the solution with the reduced model is lowered and optimal scalability of the computational algorithm with the size of the model is achieved. At the same time, numerical experiments indicate that by using reduced models, substantial savings can be achieved even in the pre‐asymptotic range. Furthermore, the algorithm parallelizes very efficiently. The method we present is expected to become a useful tool in applications requiring a large number of repeated non‐linear solid dynamics simulations, such as convergence studies, design optimization, and design of controllers of mechanical systems. Copyright © 2001 John Wiley & Sons, Ltd.     

A fast implementation of the FETI‐DP method: FETI‐DP‐RBS‐LNA and applications on large scale problems with localized non‐linearities

As parallel and distributed computing gradually becomes the computing standard for large scale problems, the domain decomposition method (DD) has received growing attention since it provides a natural basis for splitting a large problem into many small problems, which can be submitted to individual computing nodes and processed in a parallel fashion. This approach not only provides a method to solve large scale problems that are not solvable on a single computer by using direct sparse solvers but also gives a flexible solution to deal with large scale problems with localized non‐linearities. When some parts of the structure are modified, only the corresponding subdomains and the interface equation that connects all the subdomains need to be recomputed. In this paper, the dual–primal finite element tearing and interconnecting method (FETI‐DP) is carefully investigated, and a reduced back‐substitution (RBS) algorithm is proposed to accelerate the time‐consuming preconditioned conjugate gradient (PCG) iterations involved in the interface problems. Linear–non‐linear analysis (LNA) is also adopted for large scale problems with localized non‐linearities based on subdomain linear–non‐linear identification criteria. This combined approach is named as the FETI‐DP‐RBS‐LNA algorithm and demonstrated on the mechanical analyses of a welding problem. Serial CPU costs of this algorithm are measured at each solution stage and compared with that from the IBM Watson direct sparse solver and the FETI‐DP method. The results demonstrate the effectiveness of the proposed computational approach for simulating welding problems, which is representative of a large class of three‐dimensional large scale problems with localized non‐linearities. Copyright © 2005 John Wiley & Sons, Ltd.     

A SPACE-TIME COUPLED p-VERSION LEAST SQUARES FINITE ELEMENT FORMULATION FOR UNSTEADY TWO-DIMENSIONAL NAVIER–STOKES EQUATIONS

This paper presents a p-version least squares finite element formulation for two-dimensional unsteady fluid flow described by Navier–Stokes equations where the effects of space and time are coupled. The dimensionless form of the Navier–Stokes equations are first cast into a set of first-order differential equations by introducing auxiliary variables. This permits the use of C⁰ element approximation. The element properties are derived by utilizing the p-version approximation functions in both space and time and then minimizing the error functional given by the space–time integral of the sum of squares of the errors resulting from the set of first-order differential equations. This results in a true space–time coupled least squares minimization procedure. The application of least squares minimization to the set of coupled first-order partial differential equations results in finding a solution vector {δ} which makes gradient of error functional with respect to {δ} a null vector. This is accomplished by using Newton's method with a line search. A time marching procedure is developed in which the solution for the current time step provides the initial conditions for the next time step. Equilibrium iterations are carried out for each time step until the error functional and each component of the gradient of the error functional with respect to nodal degrees of freedom are below a certain prespecified tolerance. The space–time coupled p-version approximation functions provide the ability to control truncation error which, in turn, permits very large time steps. What literally requires hundreds of time steps in uncoupled conventional time marching procedures can be accomplished in a single time step using the present space–time coupled approach. The generality, success and superiority of the present formulation procedure is demonstrated by presenting specific numerical examples for transient couette flow and transient lid driven cavity. The results are compared with the analytical solutions and those reported in the literature. The formulation presented here is ideally suited for space–time adaptive procedures. The element error functional values provide a mechanism for adaptive h, p or hp refinements.     

Parallel block-preconditioned monolithic solvers for fluid-structure interaction problems

In this work, we consider the solution of fluid-structure interaction (FSI) problems using a monolithic approach for the coupling between fluid and solid subproblems. The coupling of both equations is realized by means of the arbitrary Lagrangian-Eulerian framework and a nonlinear harmonic mesh motion model. Monolithic approaches require the solution of large ill-conditioned linear systems of algebraic equations at every Newton step. Direct solvers tend to use too much memory even for a relatively small number of degrees of freedom and, in addition, exhibit superlinear growth in arithmetic complexity. Thus, iterative solvers are the only viable option. To ensure convergence of iterative methods within a reasonable amount of iterations, good and, at the same time, cheap preconditioners have to be developed. We study physics-based block preconditioners, which are derived from the block-LDU factorization of the FSI Jacobian, and their performance on distributed memory parallel computers in terms of two- and three-dimensional test cases permitting large deformations.     

Fast image inpainting using exponential-threshold POCS plus conjugate gradient

Image inpainting can remove unwanted objects and reconstruct the missing or damaged portions of an image. The projection onto convex sets (POCS) is a classical method used in image inpainting. However, the traditional POCS converges slowly due to the linear error threshold. We propose an exponential-threshold scheme, which greatly improves the convergence of the POCS. Although the exponential-threshold POCS can recover the image in about 20 iterations, it cannot reconstruct the image details very well even with hundreds of iterations. Thus, we append the non-local restoration to the exponential-threshold POCS to further refine the image details, and then we solve this objective function using the conjugate gradient. Numerical experiments show that for each iteration, the exponential-threshold POCS and the conjugate gradient have very similar computational efficiencies. For an image with various topologies of the missing areas, our scheme can recover missing pixels simultaneously and obtain a satisfied inpainting result in only 20 iterations of the exponential-threshold POCS and 20 iterations of the conjugate gradient. The proposed method can excellently restore damaged photographs and remove superimposed text. This method has less computational cost than the conjugate gradient and has a higher resolution than the POCS.     

Accurate reanalysis of structures by a preconditioned conjugate gradient method

A preconditioned conjugate gradient (PCG) method that is most suitable for reanalysis of structures is developed. The method presented provides accurate results efficiently. It is easy to implement and can be used in a wide range of applications, including non‐linear analysis and eigenvalue problems. It is shown that the PCG method presented and the combined approximations (CA) method developed recently provide theoretically identical results. Consequently, available results from one method can be applied to the other method. Effective solution procedures developed for the CA method can be used for the PCG method, and various criteria and error bounds developed for conjugate gradient methods can be used for the CA method. Numerical examples show that the condition number of the selected preconditioned matrix is much smaller than the condition number of the original matrix. This property explains the fast convergence and accurate results achieved by the method. Copyright © 2002 John Wiley & Sons, Ltd.     

ITERATIVE METHODS IN SOLVING NAVIER–STOKES EQUATIONS BY THE BOUNDARY ELEMENT METHOD

The solution of Navier–Stokes equations of time-dependent incompressible viscous fluid flow in planar geometry by the Boundary Domain Integral Method (BDIM) is discussed. The introduction of a subdomain technique to fluid flow problems is considered and improved in order to maintain the stability of BDIM. To avoid problems with flow kinematics computation in the sudomain mesh, a segmentation technique is proposed which combines the original BDIM with its subdomain variant and preserves its numerical stability. In order to reduce the computational cost of BDIM, which greatly depends on the solution of systems of linear equations, iterative methods are used. Conjugate gradient methods, conjugate gradients squared and an improved version of the biconjugate gradient method BiCGSTAB, together with the generalized minimal residual method, are used as iterative solvers. Different types of preconditioning, from simple Jacobi to incomplete LU factorization, are carried out and the performance of chosen iterative methods and preconditioners are reported. Test examples include backward facing step flow and flow through tubular heat exchangers. Test computation results show that BDIM is an accurate approximation technique which, together with the subdomain technique and powerful iterative solvers, can exhibit some significant savings in storage and CPU time requirements.     

Numerical experiments of preconditioned Krylov subspace methods solving the dense non-symmetric systems arising from BEM

Discretization of boundary integral equations leads, in general, to fully populated non-symmetric linear systems of equations. An inherent drawback of boundary element method (BEM) is that, the non-symmetric dense linear systems must be solved. For large-scale problems, the direct methods require expensive computational cost and therefore the iterative methods are perhaps more preferable. This paper studies the comparative performances of preconditioned Krylov subspace solvers as bi-conjugate gradient (Bi-CG), generalized minimal residual (GMRES), conjugate gradient squared (CGS), quasi-minimal residual (QMR) and bi-conjugate gradient stabilized (Bi-CGStab) for the solution of dense non-symmetric systems. Several general preconditioners are also considered and assessed. The results of numerical experiments suggest that the preconditioned Krylov subspace methods are effective approaches solving the large-scale dense non-symmetric linear systems arising from BEM.     

A general conjugate‐gradient‐based predictor–corrector solver for explicit finite‐element contact

A Lagrange‐multiplier based approach is presented for the general solution of multi‐body contact within an explicit finite element framework. The technique employs an explicit predictor step to permit the detection of interpenetration and then utilizes a corrector step, whose solution is obtained with a pre‐conditioned matrix‐free conjugate gradient projection method, to determine the Lagrange multipliers necessary to eliminate the predicted penetration. The predictor–corrector algorithm is developed for deformable bodies based upon the central difference method, and for rigid bodies from momentum and energy conserving approaches. Both frictionless and Coulomb‐based frictional contact idealizations are addressed. The technique imposes no time‐step constraints and quickly mitigates velocity discontinuities across closed interfaces. Special attention is directed toward contact between rigid bodies. Algorithmic moment arms conserve the translational and angular momentums of the system in the absence of external loads. Elastic collisions are captured with a two‐phase predictor–corrector approach and a geometrically approximate velocity jump criterion. The first step solves the inelastic contact problem and identifies inactive constraints between rigid bodies, while the second step generates the necessary velocity jump condition on the active constraints. The velocity criterion is shown to algorithmically preserve the system kinetic energy for two unconstrained rigid bodies. Copyright © 1999 John Wiley & Sons, Ltd. This paper was produced under the auspices of the U.S. Government and it is therefore not subject to copyright in the U.S.     

Optimal search of distributed arrival time control

Distributed arrival time control (DATC) is a heuristic feedback control algorithm for scheduling which has been developed for a real-time distributed scheduling of heterarchical systems. It has been renowned not only for its fast solution searching algorithm but also for its flexibility to changing environment. However, the optimality of this heuristic method has not been analytically explained until recently because it has been designed to discover a near optimal solution instead of the true optimum. In this paper, we provide a novel optimal search method for the DATC scheduling problem by introducing a scalar cost function over the vector space of time and show the existence and location of true optima for the DATC scheduling problem through geometric approach. This geometrical interpretation enables us to find the true optimal by direct projection without iterations like previous DATC approaches. Based on the true optimum found, we evaluate the optimality of DATC algorithms by examining their dependency on initial conditions and explain their intrinsic causality mechanism for the discrepancy with true optimum. The implication of this study is on the new viewpoint over the vector space of DATC, which not only solves the optimality issue of DATC but also provides a new direction of direct search approach like projection method for the true optimum.     

A modified Newton‐type Koiter‐Newton method for tracing the geometrically nonlinear response of structures

The Koiter‐Newton (KN) method is a combination of local multimode polynomial approximations inspired by Koiter's initial postbuckling theory and global corrections using the standard Newton method. In the original formulation, the local polynomial approximation, called a reduced‐order model, is used to make significantly more accurate predictions compared to the standard linear prediction used in conjunction with Newton method. The correction to the exact equilibrium path relied exclusively on Newton‐Raphson method using the full model. In this paper, we proposed a modified Newton‐type KN method to trace the geometrically nonlinear response of structures. The developed predictor‐corrector strategy is applied to each predicted solution of the reduced‐order model. The reduced‐order model can be used also in the correction phase, and the exact full nonlinear model is applied only to calculate force residuals. Remainder terms in both the displacement expansion and the reduced‐order model are well considered and constantly updated during correction. The same augmented finite element model system is used for both the construction of the reduced‐order model and the iterations for correction. Hence, the developed method can be seen as a particular modified Newton method with a constant iteration matrix over the single KN step. This significantly reduces the computational cost of the method. As a side product, the method has better error control, leading to more robust step size adaptation strategies. Numerical results demonstrate the effectiveness of the method in treating nonlinear buckling problems.     

Non‐linear model reduction for uncertainty quantification in large‐scale inverse problems

We present a model reduction approach to the solution of large‐scale statistical inverse problems in a Bayesian inference setting. A key to the model reduction is an efficient representation of the non‐linear terms in the reduced model. To achieve this, we present a formulation that employs masked projection of the discrete equations; that is, we compute an approximation of the non‐linear term using a select subset of interpolation points. Further, through this formulation we show similarities among the existing techniques of gappy proper orthogonal decomposition, missing point estimation, and empirical interpolation via coefficient‐function approximation. The resulting model reduction methodology is applied to a highly non‐linear combustion problem governed by an advection–diffusion‐reaction partial differential equation (PDE). Our reduced model is used as a surrogate for a finite element discretization of the non‐linear PDE within the Markov chain Monte Carlo sampling employed by the Bayesian inference approach. In two spatial dimensions, we show that this approach yields accurate results while reducing the computational cost by several orders of magnitude. For the full three‐dimensional problem, a forward solve using a reduced model that has high fidelity over the input parameter space is more than two million times faster than the full‐order finite element model, making tractable the solution of the statistical inverse problem that would otherwise require many years of CPU time. Copyright © 2009 John Wiley & Sons, Ltd.     

A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems

Several domain decomposition methods with Lagrange multipliers have been recently designed for solving iteratively large‐scale systems of finite element equations. While these methods differ typically by implementational details, they share in most cases the same substructure based preconditioners that were originally developed for the FETI method. The success of these preconditioners is due to the fact that, for homogeneous structural mechanics problems, they ensure a computational performance that scales with the problem size. In this paper, we address the suboptimal behaviour of these preconditioners in the presence of material and/or discretization heterogeneities. We propose a simple and virtually no‐cost extension of these preconditioners that exhibits scalability even for highly heterogeneous systems of equations. We consider several intricate structural analysis problems, and demonstrate numerically the optimal performance delivered by the new preconditioners for problems with discontinuities. Copyright © 1999 John Wiley & Sons, Ltd.     

Adaptive multivlevel methods in three space dimensions

We consider the approximate solution of self-adjoint elliptic problems in three space dimensions by piecewise linear finite elements with respect to a highly non-uniform tetrahedral mesh which is generated adaptively. The arising linear systems are solved iteratively by the conjugate gradient method provided with a multilevel preconditioner. Here, the accuracy of the iterative solution is coupled with the discretization error. As the performance of hierarchical bases preconditioners deteriorates in three space dimensions, the BPX preconditioner is used, taking special care of an efficient implementation. Reliable a posteriori estimates for the discretization error are derived from a local comparison with the approximation resulting from piecewise quadratic elements. To illustrate the theoretical results, we consider a familiar model problem involving reentrant corners and a real-life problem arising from hyperthermia, a recent clinical method for cancer therapy.     

The proper generalized decomposition for the simulation of delamination using cohesive zone model

The use of cohesive zone models is an efficient way to treat the damage especially when the crack path is known a priori. It is the case in the modeling of delamination in composite laminates. However, the simulations using cohesive zone models are expensive in a computational point of view. When using implicit time integration or when solving static problems, the non‐linearity related to the cohesive model requires many iteration before reaching convergence. In explicit approaches, an important number of iterations are also needed because of the time step stability condition. In this article, a new approach based on a separated representation of the solution is proposed. The proper generalized decomposition is used to build the solution. This technique coupled with a cohesive zone model allows a significant reduction of the computational cost. The results approximated with the proper generalized decomposition are very close the ones obtained using the classical finite element approach. Copyright © 2014 John Wiley & Sons, Ltd.