An Efficient Solver for Two‐way Coupling Rigid Bodies with Incompressible Flow

We present an efficient solver for monolithic two‐way coupled simulation of rigid bodies with incompressible fluids that is robust to poor conditioning of the coupled system in the presence of large density ratios between the solid and the fluid. Our method leverages ideas from the theory of Domain Decomposition, and uses a hybrid combination of direct and iterative solvers that exploits the low‐dimensional nature of the solid equations. We observe that a single Multigrid V‐cycle for the fluid equations serves as a very effective preconditioner for solving the Schur‐complement system using Conjugate Gradients, which is the main computational bottleneck in our pipeline. We use spectral analysis to give some theoretical insights behind this observation. Our method is simple to implement, is entirely assembly‐free besides the solid equations, allows for the use of large time steps because of the monolithic formulation, and remains stable even when the iterative solver is terminated early. We demonstrate the efficacy of our method on several challenging examples of two‐way coupled simulation of smoke and water with rigid bodies. To illustrate that our method is applicable to other problems, we also show an example of underwater bubble simulation.     

Additive polynomial preconditioners for parallel computers

We present a polynomial preconditioner that can be used with the conjugate gradient method to solve symmetric and positive definite systems of linear equations. Each step of the preconditioning is achieved by simultaneously taking an iteration of the SOR method and an iteration of the reverse SOR method (equations taken in reverse order) and averaging the results. This yields a symmetric preconditioner that can be implemented on parallel computers by performing the forward and reverse SOR iterations simultaneously. We give necessary and sufficient conditions for additive preconditioners to be positive definite.

We find an optimal parameter, ω, for the SOR-Additive linear stationary iterative method applied to 2-cyclic matrices. We show this method is asymptotically twice as fast as SSOR when the optimal ω is used.

We compare our preconditioner to the SSOR polynomial preconditioner for a model problem. With the optimal ω, our preconditioner was found to be as effective as the SSOR polynomial preconditioner in reducing the number of conjugate gradient iterations. Parallel implementations of both methods are discussed for vector and multiple processors. Results show that if the same number of processors are used for both preconditioners, the SSOR preconditioner is more effective. If twice as many processors are used for the SOR-Additive preconditioner, it becomes more efficient than the SSOR preconditioner when the number of equations assigned to a processor is small. These results are confirmed by the Blue Chip emulator at the University of Washington.     

Discrete Fundamental Solution Preconditioning for Hyperbolic Systems of PDE

We present a new preconditioner for the iterative solution of systems of equations arising from discretizations of systems of first order partial differential equations (PDEs) on structured grids. Such systems occur in many important applications, including compressible fluid flow and electromagnetic wave propagation. The preconditioner is a truncated convolution operator, with a kernel that is a fundamental solution of a difference operator closely related to the original discretization. Analysis of a relevant scalar model problem in two spatial dimensions shows that grid independent convergence is obtained using a simple one-stage iterative method. As an example of a more involved problem, we consider the steady state solution of the non-linear Euler equations in a two-dimensional, non-axisymmetric duct. We present results from numerical experiments, verifying that the preconditioning technique again achieves grid independent convergence, both for an upwind discretization and for a centered second order discretization with fourth order artificial viscosity.     

Improving the preconditioning of linear systems from interior point methods

This paper deals with preconditioners for solving linear systems arising from interior point methods, using iterative methods. The main focus is the development of a set of results that allows a more efficient computation of the splitting preconditioner. During the interior point methods iterations, the linear system matrix becomes ill conditioned, leading to numerical difficulties to find a solution, even with iterative methods. Therefore, the choice of an effective preconditioner is essential for the success of the approach. The paper proposes a new ordering for a splitting preconditioner, taking advantage of the sparse structure of the original matrix. A formal demonstration shows that performing this new ordering the preconditioned matrix condition number is limited; numerical experiments reinforce the theoretical results. Case studies show that the proposed idea has better sparsity features than the original version of the splitting preconditioner and that it is competitive regarding the computational time.     

An efficient solver for the RANS equations and a one-equation turbulence model

A three-stage Runge-Kutta (RK) scheme with multigrid and an implicit preconditioner has been shown to be an effective solver for the fluid dynamic equations. Using the algebraic turbulence model of Baldwin and Lomax, this scheme has been used to solve the compressible Reynolds-averaged Navier–Stokes (RANS) equations for transonic and low-speed flows. In this paper we focus on the convergence of the RK/Implicit scheme when the effects of turbulence are represented by the one-equation model of Spalart and Allmaras. With the present scheme the RANS equations and the partial differential equation of the turbulence model are solved in a loosely coupled manner. This approach allows the convergence behavior of each system to be examined. Point symmetric Gauss-Seidel supplemented with local line relaxation is used to approximate the inverse of the implicit operator of the RANS solver. To solve the turbulence equation we consider three alternative methods: diagonally dominant alternating direction implicit (DDADI), symmetric line Gauss-Seidel (SLGS), and a two-stage RK scheme with implicit preconditioning. Computational results are presented for airfoil flows, and comparisons are made with experimental data. We demonstrate that the two-dimensional RANS equations and a transport-type equation for turbulence modeling can be efficiently solved with an indirectly coupled algorithm that uses RK/Implicit schemes.     

A Matrix-free Two-grid Preconditioner for Solving Boundary Integral Equations in Electromagnetism

In this paper, we describe a matrix-free iterative algorithm based on the GMRES method for solving electromagnetic scattering problems expressed in an integral formulation. Integral methods are an interesting alternative to differential equation solvers for this problem class since they do not require absorbing boundary conditions and they mesh only the surface of the radiating object giving rise to dense and smaller linear systems of equations. However, in realistic applications the discretized systems can be very large and for some integral formulations, like the popular Electric Field Integral Equation, they become ill-conditioned when the frequency increases. This means that iterative Krylov solvers have to be combined with fast methods for the matrix-vector products and robust preconditioning to be affordable in terms of CPU time. In this work we describe a matrix-free two-grid preconditioner for the GMRES solver combined with the Fast Multipole Method. The preconditioner is an algebraic two-grid cycle built on top of a sparse approximate inverse that is used as smoother, while the grid transfer operators are defined using spectral information of the preconditioned matrix. Experiments on a set of linear systems arising from real radar cross section calculation in industry illustrate the potential of the proposed approach for solving large-scale problems in electromagnetism.     

自顶向下聚集型代数多重网格预条件的边权选择

针对基于图划分的自顶向下聚集型代数多重网格预条件,考察了利用METIS软件包进行多重网格构建的方法,并就该软件包只能处理整型权重,不能处理实型权重的问题,提出了一种将实型边权转化为整型边权的有效方法。之后将这种转化方法应用到METIS图划分软件中的边权选择,并用其给出了对自顶向下聚集型代数多重网格预条件的一种改进算法。通过对二维与三维模型偏微分方程离散所得稀疏线性方程组的数值实验表明,带边权的改进型算法大大提高了多重网格预条件共轭斜量法的迭代效率,特别是对各向异性问题,改进效果更加显著。     

Multi-dimensional gyrokinetic parameter studies based on eigenvalue computations

Plasma microinstabilities, which can be described in the framework of the linear gyrokinetic equations, are routinely computed in the context of stability analyses and transport predictions for magnetic confinement fusion experiments. The GENE code, which solves the gyrokinetic equations, has been coupled to the SLEPc package for an efficient iterative, matrix-free, and parallel computation of rightmost eigenvalues. This setup is presented, including the preconditioner which is necessary for the newly implemented Jacobi–Davidson solver. The fast computation of instabilities at a single parameter set is exploited to make parameter scans viable, that is to compute the solution at many points in the parameter space. Several issues related to parameter scans are discussed, such as an efficient parallelization over parameter sets and subspace recycling.     

A reduced unconstrained system for the cloth dynamics solver

Modern direct solvers have been more and more widely used by computer graphics community for solving sparse linear systems, such as those that arise in cloth simulation. However, external constraints usually prevent a direct method from being used for cloth simulation due to the singularity of the constrained system. This paper makes two major contributions towards the re-introduction of direct methods for cloth dynamics solvers. The first one is an approach which eliminates all the constrained variables from the system so that we obtain a reduced, nonsingular and unconstrained system. As alternatives to the well-known MPCG algorithm, not only the original, unmodified PCG method, but also any direct method can be used to solve the reduced system at a lower cost. Our second contribution is a novel direct-iterative scheme applied for the reduced system, which is basically the conjugate gradient method using a special preconditioner. Specifically, we use the stiff part of the coefficient matrix, which we call the matrix core, as the preconditioner for the PCG. The inverse of this preconditioner is computed by any eligible direct solver. The direct-iterative method has proved to be more efficient than both direct and iterative methods. Our experiments show a factor of two speedup over direct methods when stiff springs are used, even greater improvements over the MPCG iterative method.     

Preconditioned conjugate residual methods for the solution of spectral equations

Conjugate residual methods for the solution of spectral equations are described. An inexact finite-difference operator is introduced as a preconditioner in the iterative procedures. Application of these techniques is limited to problems for which the symmetric part of the coefficient matrix is positive definite. (The symmetric part of the coefficient matrix A is defined by (A + A^T)/2.) Although the spectral equation is a very ill-conditioned and full matrix problem, the computational effort of the present iterative methods for solving such a system is comparable to that for the sparse matrix equations obtained from the application of either finite-difference or finite-element methods to the same problems. Numerical experiments are shown for a self-adjoint elliptic partial differential equation with Dirichlet boundary conditions, and comparison with other solution procedures for spectral equations is presented.     

Preconditioning for a Class of Spectral Differentiation Matrices

We propose an efficient preconditioning technique for the numerical solution of first-order partial differential equations (PDEs). This study has been motivated by the computation of an invariant torus of a system of ordinary differential equations. We find the torus by discretizing a nonlinear first-order PDE with a full two-dimensional Fourier spectral method and by applying Newton’s method. This leads to large nonsymmetric linear algebraic systems. The sparsity pattern of these systems makes the use of direct solvers prohibitively expensive. Commonly used iterative methods, e.g., GMRes, BiCGStab and CGNR (Conjugate Gradient applied to the normal equations), are quite slow to converge. Our preconditioner is derived from the solution of a PDE with constant coefficients; it has a fast implementation based on the Fast Fourier Transform (FFT). It effectively increases the clustering of the spectrum, and speeds up convergence significantly. We demonstrate the performance of the preconditioner in a number of linear PDEs and the nonlinear PDE arising from the Van der Pol oscillator     

Shape design optimization of stationary fluid-structure interaction problems with large displacements and turbulence

This paper presents general and efficient methods for analysis and gradient based shape optimization of systems characterized as strongly coupled stationary fluid-structure interaction (FSI) problems. The incompressible fluid flow can be laminar or turbulent and is described using the Reynolds-averaged Navier-Stokes equations (RANS) together with the algebraic Baldwin–Lomax turbulence model. The structure may exhibit large displacements due to the interaction with the fluid domain, resulting in geometrically nonlinear structural behaviour and nonlinear interface coupling conditions. The problem is discretized using Galerkin and Streamline-Upwind/Petrov–Galerkin finite element methods, and the resulting nonlinear equations are solved using Newtons method. Due to the large displacements of the structure, an efficient update algorithm for the fluid mesh must be applied, leading to the use of an approximate Jacobian matrix in the solution routine. Expressions for Design Sensitivity Analysis (DSA) are derived using the direct differentiation approach, and the use of an inexact Jacobian matrix in the analysis leads to an iterative but very efficient scheme for DSA. The potential of gradient based shape optimization of fluid flow and FSI problems is illustrated by several examples.     

Efficient solution of fluid flow using the generalised conjugate grandient algorithm on a transputer-based machine

The discretisation of the equations governing fluid flow gives rise to coupled, quasi-linear and non-symmetric systems. The solution is usually obtained by iteration using a guess-and-correct procedure where each iteration aims to improve the solution of the previous step. Each step or outer iteration of the process involves the solution of nominally linear algebraic systems. These systems are normally solved using methods based on the Gauss-Seidel iteration—such as the TDMA. However, these methods generally converge very slowly and can be very time consuming for realistic applications. In this paper, these equations are solved using the Generalised Conjugate Gradient (GCG) algorithm with a simple-to-implement Gauss-Seidel-based preconditioner on a distributed memory message-passing machine. We take advantage of the fact that only tentative improvements to the flow-field are sought during each iteration and study the convergence behaviour of the parallel implementation on a multi-processor environment.     

On the performance of SPAI and ADI-like preconditioners for core collapse supernova simulations in one spatial dimension

The simulation of core collapse supernovæ calls for the time accurate solution of the (Euler) equations for inviscid hydrodynamics coupled with the equations for neutrino transport. The time evolution is carried out by evolving the Euler equations explicitly and the neutrino transport equations implicitly. Neutrino transport is modeled by the multi-group Boltzmann transport (MGBT) and the multi-group flux limited diffusion (MGFLD) equations. An implicit time stepping scheme for the MGBT and MGFLD equations yields Jacobian systems that necessitate scaling and preconditioning. Two types of preconditioners, namely, a sparse approximate inverse (SPAI) preconditioner and a preconditioner based on the alternating direction implicit iteration (ADI-like) have been found to be effective for the MGFLD and MGBT formulations. This paper compares these two preconditioners. The ADI-like preconditioner performs well with both MGBT and MGFLD systems. For the MGBT system tested, the SPAI preconditioner did not give competitive results. However, since the MGBT system in our experiments had a high condition number before scaling and since we used a sequential platform, care must be taken in evaluating these results.     

Preconditioned conjugate gradients,radial basis functions,and Toeplitz matrices

Radial basis functions provide highly useful and flexible interpolants to multivariate functions. Further, they are beginning to be used in the numerical solution of partial differential equations. Unfortunately, their construction requires the solution of a dense linear system. Therefore, much attention has been given to iterative methods. In this paper, we present a highly efficient preconditioner for the conjugate gradient solution of the interpolation equations generated by gridded data. Thus, our method applies to the corresponding Toeplitz matrices. The number of iterations required to achieve a given tolerance is independent of the number of variables.     

Structure Preserving Preconditioners for Image Deblurring

Regularizing preconditioners for accelerating the convergence of iterative regularization methods without spoiling the quality of the approximated solution have been extensively investigated in the last twenty years. Several strategies have been proposed for defining proper preconditioners. Usually, in methods for image restoration, the structure of the preconditioner is chosen Block Circulant with Circulant Blocks (BCCB) because it can be efficiently exploited by Fast Fourier Transform (FFT). Nevertheless, for ill-conditioned problems, it is well-known that BCCB preconditioners cannot provide a strong clustering of the eigenvalues. Moreover, in order to get an effective preconditioner, it is crucial to preserve the structure of the coefficient matrix. The structure of such a matrix, in case of image deblurring problem, depends on the boundary conditions imposed on the imaging model. Therefore, we propose a technique to construct a preconditioner which has the same structure of the blurring matrix related to the restoration problem at hand. The construction of our preconditioner requires two FFTs like the BCCB preconditioner. The presented preconditioning strategy represents a generalization and an improvement with respect to both circulant and structured preconditioning available in the literature. The technique is further extended to provide a non-stationary preconditioning in the same spirit of a recent proposal for BCCB matrices. Some numerical results show the importance of preserving the matrix structure from the point of view of both restoration quality and robustness of the regularization parameter.     

Stable discretization of poroelasticity problems and efficient preconditioners for arising saddle point type matrices

Poroelastic models arise in reservoir modeling and many other important applications. Under certain assumptions, they involve a time-dependent coupled system consisting of Navier–Lamé equations for the displacements, Darcy’s flow equation for the fluid velocity and a divergence constraint equation. Stability for infinite time of the continuous problem and, second and third order accurate, time discretized equations are shown. Methods to handle the lack of regularity at initial times are discussed and illustrated numerically. After discretization, at each time step this leads to a block matrix system in saddle point form. Mixed space discretization methods and a regularization method to stabilize the system and avoid locking in the pressure variable are presented. A certain block matrix preconditioner is shown to cluster the eigenvalues of the preconditioned matrix about the unit value but needs inner iterations for certain matrix blocks. The strong clustering leads to very few outer iterations. Various approaches to construct preconditioners are presented and compared. The sensitivity of the number of outer iterations to the stopping accuracy of the inner iterations is illustrated numerically.     

Supercloseness of Continuous Interior Penalty Methods on Shishkin Triangular Meshes and Hybrid Meshes for Singularly Perturbed Problems with Characteristic Layers

A fast preconditioned penalty method is developed for a system of parabolic linear complementarity problems (LCPs) involving tempered fractional order partial derivatives governing the price of American options whose underlying asset follows a geometry Lévy process with multi-state regime switching. By means of the penalty method, the system of LCPs is approximated with a penalty term by a system of nonlinear tempered fractional partial differential equations (TFPDEs) coupled by a finite-state Markov chain. The system of nonlinear TFPDEs is discretized with the shifted Grünwald approximation by an upwind finite difference scheme which is shown to be unconditionally stable. Semi-smooth Newton’s method is utilized to solve the finite difference scheme as an outer iterative method in which the Jacobi matrix is found to possess Toeplitz-plus-diagonal structure. Consequently, the resulting linear system can be fast solved by the Krylov subspace method as an inner iterative method via fast Fourier transform (FFT). Furthermore, a novel preconditioner is proposed to speed up the convergence rate of the inner Krylov subspace iteration with theoretical analysis. With the above-mentioned preconditioning technique via FFT, under some mild conditions, the operation cost in each Newton’s step can be expected to be \(\mathcal{O}(N\mathrm{log}N)\), where N is the size of the coefficient matrix. Numerical examples are given to demonstrate the accuracy and efficiency of our proposed fast preconditioned penalty method.     

SOLVING SPARSE LEAST SQUARES PROBLEMS WITH PRECONDITIONED CGLS METHOD ON PARALLEL DISTRIBUTED MEMORY COMPUTERS

Abstract

In this paper we study the parallel aspects of PCGLS, a basic iterative method based on the conjugate gradient method with preconditioner applied to normal equations and Incomplete Modified Gram-Schmidt (IMGS) preconditioner, for solving sparse least squares problems on massively parallel distributed memory computers. The performance of these methods on this kind of architecture is usually limited because of the global communication required for the inner products. We will describe the parallelization of PCGLS and IMGS preconditioner by two ways of improvement. One is to accumulate the results of a number of inner products collectively and the other is to create situations where communication can be overlapped with computation. A theoretical model of computation and communication phases is presented which allows us to determine the optimal number of processors that minimizes the runtime. Several numerical experiments on the Parsytec GC/PowerPlus are presented.     

Convergence analysis of a subdomain iterative method for the finite element approximation of the coupling of Stokes and Darcy equations

We consider a Galerkin Finite Element approximation of the Stokes-Darcy problem which models the coupling between surface and groundwater flows. Then we propose an iterative subdomain method for its solution, inspired to the domain decomposition theory. The convergence analysis that we develop is based on the properties of the discrete Steklov-Poincaré operators associated to the given coupled problem. An optimal preconditioner for Krylov methods is proposed and analyzed.