Some observations on weighted GMRES

We investigate the convergence of the weighted GMRES method for solving linear systems. Two different weighting variants are compared with unweighted GMRES for three model problems, giving a phenomenological explanation of cases where weighting improves convergence, and a case where weighting has no effect on the convergence. We also present new alternative implementations of the weighted Arnoldi algorithm which may be favorable in terms of computational complexity, and examine stability issues connected with these implementations. Two implementations of weighted GMRES are compared for a large number of examples. We find that weighted GMRES may outperform unweighted GMRES for some problems, but more often this method is not competitive with other Krylov subspace methods like GMRES with deflated restarting or BICGSTAB, in particular when a preconditioner is used

Iterative Solution Methods for Reduced-Order Models of Parameterized Partial Differential Equations

This dissertation considers efficient computational algorithms for solving parameterized discrete partial differential equations (PDEs) using techniques of reduced-order modeling. Parameterized equations of this type arise in numerous mathematical models. In some settings, e.g. sensitivity analysis, design optimization, and uncertainty quantification, it is necessary to compute discrete solutions of the PDEs at many parameter values. Accuracy considerations often lead to algebraic systems with many unknowns whose solution via traditional methods can be expensive. Reduced-order models use a reduced space to approximate the parameterized PDE, where the reduced space is of a significantly smaller dimension than that of the discrete PDE. Solving an approximation of the problem on the reduced space leads to reduction in cost, often with little loss of accuracy. In the reduced basis method, an offline step finds an approximation of the solution space and an online step utilizes this approximation to solve a smaller reduced problem, which provides an accurate estimate of the solution. Traditionally, the reduced problem is solved using direct methods. However, the size of the reduced system needed to produce solutions of a given accuracy depends on the characteristics of the problem, and it may happen that the size is significantly smaller than that of the original discrete problem but large enough to make direct solution costly. In this scenario, it is more effective to use iterative methods to solve the reduced problem. To demonstrate this we construct preconditioners for the reduced-order models or construct well-conditioned reduced-order models. We demonstrate that by using iterative methods, reduced-order models of larger dimension can be effective. There are several reasons that iterative methods are well suited to reduced- order modeling. In particular, we take advantage of the similarity of the realizations of parameterized systems, either by reusing preconditioners or by recycling Krylov vectors. These two approaches are shown to be effective when the underlying PDE is linear. For nonlinear problems, we utilize the discrete empirical interpolation method (DEIM) to cheaply evaluate the nonlinear components of the reduced model. The method identifies points in the PDE discretization necessary for representing the nonlinear component of the reduced model accurately. This approach incurs online computational costs that are independent of the spatial dimension of the discretized PDE. When this method is used to assemble the reduced model cheaply, iterative methods are shown to further improve efficiency in the online step. Finally, when the traditional offline/online approach is ineffective for a given problem, reduced-order models can be used to accelerate the solution of the full model. We follow the solution model of Krylov subspace recycling methods for sequences of linear systems where the coefficient matrices vary. A Krylov subspace recycling method contains a reduced-order model and an iterative method that searches the space orthogonal to the reduced space. We once again use iterative solution techniques for the solution of the reduced models that arise in this context. In this case, the iterative methods converge quickly when the reduced basis is constructed to be naturally well conditioned

Global convection in Earth's mantle : advanced numerical methods and extreme-scale simulations

The thermal convection of rock in Earth's mantle and associated plate tectonics are modeled by nonlinear incompressible Stokes and energy equations. This dissertation focuses on the development of advanced, scalable linear and nonlinear solvers for numerical simulations of realistic instantaneous mantle flow, where we must overcome several computational challenges. The most notable challenges are the severe nonlinearity, heterogeneity, and anisotropy due to the mantle's rheology as well as a wide range of spatial scales and highly localized features. Resolving the crucial small scale features efficiently necessitates adaptive methods, while computational results greatly benefit from a high accuracy per degree of freedom and local mass conservation. Consequently, the discretization of Earth's mantle is carried out by high-order finite elements on aggressively adaptively refined hexahedral meshes with a continuous, nodal velocity approximation and a discontinuous, modal pressure approximation. These velocity--pressure pairings yield optimal asymptotic convergence rates of the finite element approximation to the infinite-dimensional solution with decreasing mesh element size, are inf-sup stable on general, non-conforming hexahedral meshes with "hanging nodes,'' and have the advantage of preserving mass locally at the element level due to the discontinuous pressure. However, because of the difficulties cited above and the desired accuracy, the large implicit systems to be solved are extremely poorly conditioned and sophisticated linear and nonlinear solvers including powerful preconditioning techniques are required. The nonlinear Stokes system is solved using a grid continuation, inexact Newton--Krylov method. We measure the residual of the momentum equation in the H⁻¹-norm for backtracking line search to avoid overly conservative update steps that are significantly reduced from one. The Newton linearization is augmented by a perturbation of a highly nonlinear term in mantle's rheology, resulting in dramatically improved nonlinear convergence. We present a new Schur complement-based Stokes preconditioner, weighted BFBT, that exhibits robust fast convergence for Stokes problems with smooth but highly varying (up to 10 orders of magnitude) viscosities, optimal algorithmic scalability with respect to mesh refinement, and only a mild dependence on the polynomial order of high-order finite element discretizations. In addition, we derive theoretical eigenvalue bounds to prove spectral equivalence of our inverse Schur complement approximation. Finally, we present a parallel hybrid spectral--geometric--algebraic multigrid (HMG) to approximate the inverses of the Stokes system's viscous block and variable-coefficient pressure Poisson operators within weighted BFBT. Building on the parallel scalability of HMG, our Stokes solver demonstrates excellent parallel scalability to 1.6 million CPU cores without sacrificing algorithmic optimality.Computational Science, Engineering, and Mathematic

An embedded boundary integral solver for the stokes equations

We present a new method for the solution of the Stokes equations. Our goal is to develop a robust and scalable methodology for two and three dimensional, moving-boundary, flow simulations. Our method is based on Anita Mayo\u27s method for the Poisson\u27s equation: “The Fast Solution of Poisson\u27s and the Biharmonic Equations on Irregular Regions”, SIAM J. Num. Anal., 21 (1984), pp. 285– 299. We embed the domain in a rectangular domain, for which fast solvers are available, and we impose the boundary conditions as interface (jump) conditions on the velocities and tractions. We use an indirect boundary integral formulation for the homogeneous Stokes equations to compute the jumps. The resulting integral equations are discretized by Nystrom\u27s method. The rectangular domain problem is discretized by finite elements for a velocity-pressure formulation with equal order interpolation bilinear elements (Q1-Q1). Stabilization is used to circumvent the inf-sup condition for the pressure space. For the integral equations, fast matrix vector multiplications are achieved via a NlogN algorithm based on a block representation of the discrete integral operator, combined with (kernel independent) singular value decomposition to sparsify low-rank blocks. Our code is built on top of PETSc, an MPI based parallel linear algebra library. The regular grid solver is a Krylov method (Conjugate Residuals) combined with an optimal two-level Schwartz-preconditioner. For the integral equation we use GMRES. We have tested our algorithm on several numerical examples and we have observed optimal convergence rates

Next-Generation Computational Fluid Dynamics Capability for Aircraft Aeroelasticity and Loads

This paper presents some of the first results obtained from the recently implemented linear frequency domain solver in CODA, the next generation flow solver framework. It uses automatic differentiation capability to compute the exact product of the Jacobian matrix with an arbitrary vector. Perturbation results for subsonic and transonic LANN wing cases show good agreement with experiments and results computed using DLR-TAU code. Computations have been performed using both one- and two-equation turbulence models for the NASA Common Research Model. Scalability assessment of the frequency-domain solver demonstrates the advantage of hybrid MPI/OpenMP partitioning over pure MPI partitioning