Staggered discontinuous Galerkin methods for the incompressible Navier–Stokes equations: Spectral analysis and computational results

The goal of this paper is to create a fruitful bridge between the numerical methods for approximating PDEs in fluid dynamics and the (iterative) numerical methods for dealing with the resulting large linear systems. Among the main objectives are the design of new, efficient iterative solvers and a rigorous analysis of their convergence speed. The link we have in mind is either the structure or the hidden structure that the involved coefficient matrices inherit, both from the continuous PDE and from the approximation scheme; in turn, the resulting structure is used for deducing spectral information, crucial for the conditioning and convergence analysis and for the design of more efficient solvers. As a specific problem, we consider the incompressible Navier–Stokes equations; as a numerical technique, we consider a novel family of high‐order, accurate discontinuous Galerkin methods on staggered meshes, and as tools, we use the theory of Toeplitz matrices generated by a function (in the most general block, the multilevel form) and the more recent theory of generalized locally Toeplitz matrix sequences. We arrive at a somehow complete picture of the spectral features of the underlying matrices, and this information is employed for giving a forecast of the convergence history of the conjugate gradient method, together with a discussion on new and more advanced techniques (involving preconditioning, multigrid, multi‐iterative solvers). Several numerical tests are provided and critically illustrated in order to show the validity and the potential of our analysis.     

Strengthened Cauchy–Bunyakowski–Schwarz inequality for a three‐dimensional elasticity system

The constant γ of the strengthened Cauchy–Bunyakowski–Schwarz (CBS) inequality plays a fundamental role in the convergence rate of multilevel iterative methods. The main purpose of this work is to give an estimate of the constant γ for a three‐dimensional elasticity system. The theoretical results obtained are practically important for the successful implementation of the finite element method to large‐scale modelling of complicated structures as they allow us to construct optimal order algebraic multilevel iterative solvers for a wide class of real‐life elasticity problems. Copyright © 2001 John Wiley & Sons, Ltd.     

Optimal scaling parameters for sparse grid discretizations

We apply iterative subspace correction methods to elliptic PDE problems discretized by generalized sparse grid systems. The involved subspace solvers are based on the combination of all anisotropic full grid spaces that are contained in the sparse grid space. Their relative scaling is at our disposal and has significant influence on the performance of the iterative solver. In this paper, we follow three approaches to obtain close‐to‐optimal or even optimal scaling parameters of the subspace solvers and thus of the overall subspace correction method. We employ a Linear Program that we derive from the theory of additive subspace splittings, an algebraic transformation that produces partially negative scaling parameters that result in improved asymptotic convergence properties, and finally, we use the OptiCom method as a variable nonlinear preconditioner. Copyright © 2014 John Wiley & Sons, Ltd.     

Recent advances in linear analysis of convergence for splittings for solving ODE problems

In the nineties, Van der Houwen et al. (see, e.g., [P.J. van der Houwen, B.P. Sommeijer, J.J. de Swart, Parallel predictor–corrector methods, J. Comput. Appl. Math. 66 (1996) 53–71; P.J. van der Houwen, J.J.B. de Swart, Triangularly implicit iteration methods for ODE-IVP solvers, SIAM J. Sci. Comput. 18 (1997) 41–55; P.J. van der Houwen, J.J.B. de Swart, Parallel linear system solvers for Runge–Kutta methods, Adv. Comput. Math. 7 (1–2) (1997) 157–181]) introduced a linear analysis of convergence for studying the properties of the iterative solution of the discrete problems generated by implicit methods for ODEs. This linear convergence analysis is here recalled and completed, in order to provide a useful quantitative tool for the analysis of splittings for solving such discrete problems. Indeed, this tool, in its complete form, has been actively used when developing the computational codes BiM and BiMD [L. Brugnano, C. Magherini, The BiM code for the numerical solution of ODEs, J. Comput. Appl. Math. 164–165 (2004) 145–158. Code available at: http://www.math.unifi.it/~brugnano/BiM/index.html; L. Brugnano, C. Magherini, F. Mugnai, Blended implicit methods for the numerical solution of DAE problems, J. Comput. Appl. Math. 189 (2006) 34–50]. Moreover, the framework is extended for the case of special second order problems. Examples of application, aimed to compare different iterative procedures, are also presented.     

Coupling finite difference methods and integral formulas for elliptic problems arising in fluid mechanics

This article is devoted to the numerical analysis of two classes of iterative methods that combine integral formulas with finite‐difference Poisson solvers for the solution of elliptic problems. The first method is in the spirit of the Schwarz domain decomposition method for exterior domains. The second one is motivated by potential calculations in free boundary problems and can be viewed as a numerical analytic continuation algorithm. Numerical tests are presented that confirm the convergence properties predicted by numerical analysis. © 2003 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 20: 199–229, 2004     

Fast Fourier Transform Solvers and Preconditioners for Quadratic Spline Collocation

Quadratic Spline Collocation (QSC) methods of optimal order of convergence have been recently developed for the solution of elliptic Partial Differential Equations (PDEs). In this paper, linear solvers based on Fast Fourier Transforms (FFT)are developed for the solution of the QSC equations. The complexity of the FFT solvers is O(N ² log N), where N is the gridsize in one dimension. These direct solvers can handle PDEs with coefficients in one variable or constant, and Dirichlet, Neumann, alternating Dirichlet-Neumann or periodic boundary conditions, along at least one direction of a rectangular domain. General variable coefficient PDEs are handled by preconditioned iterative solvers. The preconditioner is the QSC matrix arising from a constant coefficient PDE. The convergence analysis of the preconditioner is presented. It is shown that, under certain conditions, the convergence rate is independent of the gridsize. The preconditioner is solved by FFT techniques, and integrated with one-step or acceleration methods, giving rise to asymptotically almost optimal linear solvers, with complexity O(N ² log N). Numerical experiments verify the effectiveness of the solvers and preconditioners, even on problems more general than the analysis assumes. The development and analysis of FFT solvers and preconditioners is extended to QSC equations corresponding to systems of elliptic PDEs.     

Nonlinear structural finite element analysis using the preconditioned Lanczos method on serial and parallel computers

The application of the Lanczos algorithm in Newton-like methods for solving non-linear systems of equations arising in nonlinear structural finite element analysis is presented. It is shown that with appropriate preconditioners iterative methods can be developed which are robust and efficient even for ill conditioned problems. Though the real advantage of iterative solvers seems to exist on distributed memory machines, even on serial machines the performance can be improved compared with direct solvers while saving memory capacity. With a specific modification of the Lanczos algorithm in combination with arc-length procedures a further speed-up of the nonlinear analysis can be achieved. For parallel implementations domain decomposition methods are used. A parallel preconditioning strategy based on an incomplete factorisation method is presented. An example is taken and the quality and efficiency of two different domain decomposition methods are discussed for a large shell structure. This work was supported by the BMBF (Bundesministerium für Bildung und Forschung) of Germany.     

Upper bound of the constant in the strengthened C.B.S. inequality for FEM 2D elasticity equations

The basic theory of the strengthened Cauchy–Buniakowskii–Schwarz (C.B.S.) inequality is the main tool in the convergence analysis of the recently proposed algebraic multilevel iterative methods. An upper bound of the constant γ in the strengthened C.B.S. inequality for the case of the finite element solution of 2D elasticity problems is obtained. It is assumed that linear triangle finite elements are used, the initial mesh consisting of right isosceles triangles and the mesh refinement procedure being uniform. For the resulting linear algebraic systems we have proved that γ²<0.75 uniformly on the mesh parameter and on Poisson's ratio ν ? (0, 1/2). Furthermore, the presented numerical tests show that the same relation holds for arbitrary initial right triangulations, even in the case of degeneracy of triangles. The theoretical results obtained are practically important for successful implementation of the finite element method to large-scale modeling of complicated structures. They allow us to construct optimal order algebraic multilevel iterative solvers for a wide class of real–life elasticity problems.     

A domain-decomposing parallel sparse linear system solver

The solution of large sparse linear systems is often the most time-consuming part of many science and engineering applications. Computational fluid dynamics, circuit simulation, power network analysis, and material science are just a few examples of the application areas in which large sparse linear systems need to be solved effectively. In this paper, we introduce a new parallel hybrid sparse linear system solver for distributed memory architectures that contains both direct and iterative components. We show that by using our solver one can alleviate the drawbacks of direct and iterative solvers, achieving better scalability than with direct solvers and more robustness than with classical preconditioned iterative solvers. Comparisons to well-known direct and iterative solvers on a parallel architecture are provided.     

Robust Multiscale Iterative Solvers for Nonlinear Flows in Highly Heterogeneous Media

In this paper, we study robust iterative solvers for finite element systems resulting in approximation of steady-state Richards' equation in porous media with highly heterogeneous conductivity fields. It is known that in such cases the contrast, ratio between the highest and lowest values of the conductivity, can adversely affect the performance of the preconditioners and, consequently, a design of robust preconditioners is important for many practical applications. The proposed iterative solvers consist of two kinds of iterations, outer and inner iterations. Outer iterations are designed to handle nonlinearities by linearizing the equation around the previous solution state. As a result of the linearization, a large-scale linear system needs to be solved. This linear system is solved iteratively (called inner iterations), and since it can have large variations in the coefficients, a robust preconditioner is needed. First, we show that under some assumptions the number of outer iterations is independent of the contrast. Second, based on the recently developed iterative methods, we construct a class of preconditioners that yields convergence rate that is independent of the contrast. Thus, the proposed iterative solvers are optimal with respect to the large variation in the physical parameters. Since the same preconditioner can be reused in every outer iteration, this provides an additional computational savings in the overall solution process. Numerical tests are presented to confirm the theoretical results.     

A domain-decomposing parallel sparse linear system solver

The solution of large sparse linear systems is often the most time-consuming part of many science and engineering applications. Computational fluid dynamics, circuit simulation, power network analysis, and material science are just a few examples of the application areas in which large sparse linear systems need to be solved effectively. In this paper, we introduce a new parallel hybrid sparse linear system solver for distributed memory architectures that contains both direct and iterative components. We show that by using our solver one can alleviate the drawbacks of direct and iterative solvers, achieving better scalability than with direct solvers and more robustness than with classical preconditioned iterative solvers. Comparisons to well-known direct and iterative solvers on a parallel architecture are provided.     

Efficient linear solvers for incompressible flow simulations using Scott‐Vogelius finite elements

Recent research has shown that in some practically relevant situations like multiphysics flows (Galvin et al., Comput Methods Appl Mech Eng, to appear) divergence‐free mixed finite elements may have a significantly smaller discretization error than standard nondivergence‐free mixed finite elements. To judge the overall performance of divergence‐free mixed finite elements, we investigate linear solvers for the saddle point linear systems arising in ((P_k)^d,P _k‐1^disc) Scott‐Vogelius finite element implementations of the incompressible Navier–Stokes equations. We investigate both direct and iterative solver methods. Due to discontinuous pressure elements in the case of Scott‐Vogelius (SV) elements, considerably more solver strategies seem to deliver promising results than in the case of standard mixed finite elements such as Taylor‐Hood elements. For direct methods, we extend recent preliminary work using sparse banded solvers on the penalty method formulation to finer meshes and discuss extensions. For iterative methods, we test augmented Lagrangian and \begin{align*}\mathcal{H}\end{align*} ‐LU preconditioners with GMRES, on both full and statically condensed systems. Several numerical experiments are provided that show these classes of solvers are well suited for use with SV elements and could deliver an interesting overall performance in several applications.© 2012 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2013     

An efficient hierarchical preconditioner for quadratic discretizations of finite element problems

Higher order finite element discretizations, although providing higher accuracy, are considered to be computationally expensive and of limited use for large‐scale problems. In this paper, we have developed an efficient iterative solver for solving large‐scale quadratic finite element problems. The proposed approach shares some common features with geometric multigrid methods but does not need structured grids to create the coarse problem. This leads to a robust method applicable to finite element problems discretized by unstructured meshes such as those from adaptive remeshing strategies. The method is based on specific properties of hierarchical quadratic bases. It can be combined with an algebraic multigrid (AMG) preconditioner or with other algebraic multilevel block factorizations. The algorithm can be accelerated by flexible Krylov subspace methods. We present some numerical results on the convection–diffusion and linear elasticity problems to illustrate the efficiency and the robustness of the presented algorithm. In these experiments, the performance of the proposed method is compared with that of an AMG preconditioner and other iterative solvers. Our approach requires less computing time and less memory storage. Copyright © 2010 John Wiley & Sons, Ltd.     

Generalized Jacobi and Gauss-Seidel Methods for Solving Linear System of Equations

The Jacobi and Gauss-Seidel algorithms are among the stationary iterative methods for solving linear system of equations. They are now mostly used as precondition-ers for the popular iterative solvers. In this paper a generalization of these methods are proposed and their convergence properties are studied. Some numerical experiments are given to show the efficiency of the new methods.     

Inexact Methods: Forcing Terms and Conditioning

In this paper, we consider inexact Newton and Newton-like methods andprovide new convergence conditions relating the forcing terms to theconditioning of the iteration matrices. These results can be exploited wheninexact methods with iterative linear solvers are used. In this framework,preconditioning techniques can be used to improve the performance ofiterative linear solvers and to avoid the need of excessively small forcingterms. Numerical experiments validating the theoretical results arediscussed.     

A comparative study of efficient iterative solvers for generalized Stokes equations

We consider a generalized Stokes equation with problem parameters ξ?0 (size of the reaction term) and ν>0 (size of the diffusion term). We apply a standard finite element method for discretization. The main topic of the paper is a study of efficient iterative solvers for the resulting discrete saddle point problem. We investigate a coupled multigrid method with Braess–Sarazin and Vanka‐type smoothers, a preconditioned MINRES method and an inexact Uzawa method. We present a comparative study of these methods. An important issue is the dependence of the rate of convergence of these methods on the mesh size parameter and on the problem parameters ξ and ν. We give an overview of the main theoretical convergence results known for these methods. For a three‐dimensional problem, discretized by the Hood–Taylor ??₂–??₁ pair, we give results of numerical experiments. Copyright © 2007 John Wiley & Sons, Ltd.     

Algebraic multigrid preconditioning for iterative eigensolvers

The objective is a comparative study of iterative solvers for eigenproblems arising from elliptic and self–adjoint partial differential operators. Typically only a few of the smallest eigenvalues of these problems are to be computed. We discuss various gradient based preconditioned eigensolvers which make use of algebraic multigrid preconditioning. We present some algorithms together with numerical results. Performance characteristics are derived by comparison with the solutions of standard problems. We show that known advantages of algebraic multigrid preconditioning (e.g. for boundary–value problems with large jumps in the coefficients) transfer to AMG–preconditioned eigensolvers. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Numerical recovery strategies for parallel resilient Krylov linear solvers

As the computational power of high‐performance computing systems continues to increase by using a huge number of cores or specialized processing units, high‐performance computing applications are increasingly prone to faults. In this paper, we present a new class of numerical fault tolerance algorithms to cope with node crashes in parallel distributed environments. This new resilient scheme is designed at application level and does not require extra resources, that is, computational unit or computing time, when no fault occurs. In the framework of iterative methods for the solution of sparse linear systems, we present numerical algorithms to extract relevant information from available data after a fault, assuming a separate mechanism ensures the fault detection. After data extraction, a well‐chosen part of missing data is regenerated through interpolation strategies to constitute meaningful inputs to restart the iterative scheme. We have developed these methods, referred to as interpolation–restart techniques, for Krylov subspace linear solvers. After a fault, lost entries of the current iterate computed by the solver are interpolated to define a new initial guess to restart the Krylov method. A well‐suited initial guess is computed by using the entries of the faulty iterate available on surviving nodes. We present two interpolation policies that preserve key numerical properties of well‐known linear solvers, namely, the monotonic decrease of the A‐norm of the error of the conjugate gradient or the residual norm decrease of generalized minimal residual algorithm for solving. The qualitative numerical behavior of the resulting scheme has been validated with sequential simulations, when the number of faults and the amount of data losses are varied. Finally, the computational costs associated with the recovery mechanism have been evaluated through parallel experiments. Copyright © 2016 John Wiley & Sons, Ltd.     

Linear and nonlinear substructured Restricted Additive Schwarz iterations and preconditioning

Iterative substructuring Domain Decomposition (DD) methods have been extensively studied, and they are usually associated with nonoverlapping decompositions. It is less known that classical overlapping DD methods can also be formulated in substructured form, i.e., as iterative methods acting on variables defined exclusively on the interfaces of the overlapping domain decomposition. We call such formulations substructured domain decomposition methods. We introduce here a substructured version of Restricted Additive Schwarz (RAS) which we call SRAS. We show that RAS and SRAS are equivalent when used as iterative solvers, as they produce the same iterates, while they are substantially different when used as preconditioners for GMRES. We link the volume and substructured Krylov spaces and show that the iterates are different by deriving the least squares problems solved at each GMRES iteration. When used as iterative solvers, SRAS presents computational advantages over RAS, as it avoids computations with matrices and vectors at the volume level. When used as preconditioners, SRAS has the further advantage of allowing GMRES to store smaller vectors and perform orthogonalization in a lower dimensional space. We then consider nonlinear problems, and we introduce SRASPEN (Substructured Restricted Additive Schwarz Preconditioned Exact Newton), where SRAS is used as a preconditioner for Newton’s method. In contrast to the linear case, we prove that Newton’s method applied to the preconditioned volume and substructured formulation produces the same iterates in the nonlinear case. Next, we introduce two-level versions of nonlinear SRAS and SRASPEN. Finally, we validate our theoretical results with numerical experiments.