Developement of Algebraic Multigrid Solvers Using GPUs

We propose a new, efficient, adaptive algebraic multigrid (AMG) solver strategy for the discrete systems of partial differential equations arising from structured or unstructured grid models in reservoir simulation. The proposed strategy has been particularly tailored to linear systems of equations arising in adaptive implicit methods. The coarsening process of the AMG method designed automatically employs information on the physical structure of the models; as a smoother, an adaptive ILUT method is employed, taking care of an efficient solution of the hyperbolic parts whilst providing adequately smooth errors for the elliptic parts. To achieve a good compromise of high efficiency and robustness for a variety of problem classes – ranging from simple, small black-oil to challenging, large compositional models – an automatic, adaptive ILUT parameter and AMG solver switching strategy, α-SAMG, has been developed. Its efficiency is demonstrated for eight industrial benchmark cases by comparison against standard one-level and AMG solvers as well as the pure one- level variant of the proposed new strategy. In addition, very promising results of first parallel runs are shown.

A new high-order compact (HOC) cascadic multigrid (MG) strategy is developed for modeling large-scale complex gravity problems on nonuniform rectilinear grids. First, the governing 3D Poisson’s gravitational potential equation is discretized by a novel HOC difference scheme, which allows for Robin boundary conditions to further reduce the computational domain. Then, the resulting asymmetric linear system of equations is solved by a generalized extrapolation cascadic MG (gEXCMG) method augmented with a novel MG prolongation operator. The quartic interpolation and completed Richardson extrapolation assist this MG prolongation operator to produce a good initial guess for the biconjugate gradient stabilized smoother. Two synthetic and SEG/EAGE salt models are examined to verify the presented approach. Results indicate that our algorithm can efficiently and accurately produce the gravity signals for complex geologic models. Moreover, the comparisons with the state-of-the-art algebraic MG solvers demonstrate that our gEXCMG solver offers substantially better efficiency with significantly less memory consumption. Therefore, our newly developed method has the potential for serving as a forward engine in large-scale gravity inversions.

Commercial reservoir simulators must be very robust and fast. Moreover, current hardware requires the simulators to scale over multiple number of computing nodes and for a fixed (‘strong scalability’) as well as an increasing problem size per computing node (‘weak scalability’). In most current commercial reservoir simulators, due to the different geological structures and properties of hydrocarbon reservoirs and the use of enhanced oil recovery (EOR) techniques, the governing equations are strongly nonlinear and hard to solve. The Jacobian system is solved by FGMRES preconditioned by the two-level constrained pressure residual (CPR) preconditioner. The driving force of the CPR preconditioner is the solution of the pressure equation. The industry standard for solving the pressure equation is the algebraic multigrid (AMG) solver. AMG is well known for its ‘weak scalability’. However, in these applications, AMG has unfavorable ‘strong’ scalability properties. This degradation in scalability is due to the increased level of inter-processor communication in the algorithm. In this paper, a monotone non-Galerkin AMG (MNG-AMG) method is presented. The aim of the method is to reduce the overall communication in MNG-AMG by enforcing a predefined nonzero pattern and monotonicity property (i.e., M-matrices) on each multigrid level. This paper describes the application of the MNG-AMG method in the context of reservoir simulations. We will compare the parallel scalability of the default solver with the MNG-AMG solver and discuss the optimal values for the MNG-AMG solver for a variety of test cases based on full field reservoir simulations.

Efficient forward modeling of electromagnetic (EM) fields is the basis of interpretation and inversion of the EM data, which plays an important role in practical geophysical exploration. A novel extrapolation cascadic multigrid (EXCMG) method is developed to solve large linear systems encountered in geophysical EM modeling. The original curl-curl equation with arbitrary anisotropic conductivity is regularized by including the gradient of a scaled divergence correction term, and a linear edge element method is used to discretize the equation. We develop a novel approach to address the issue of edge unknowns in 3D edge element discretizations on nonuniform rectilinear grids. Inspired by the original EXCMG method for nodal elements, we introduce a new prolongation operator that treats edge unknowns as defined on the midpoints of edges. This operator aims to provide an accurate approximation of the finite-element solution on the refined grid. Using a good initial guess significantly reduces the number of iterations required by the preconditioned biconjugate gradient stabilized solver, which is used as a smoother for the EXCMG algorithm. Numerical experiments are carried out to validate the accuracy and efficiency of our EXCMG method, including models from magnetotellurics and controlled-source EM modeling. The testing results indicate that EXCMG is more efficient than traditional Krylov-subspace iterative solvers, the algebraic multigrid solver, and those depending on the auxiliary-space Maxwell solver, especially for large-scale problems where the number of unknowns exceeds 10 million. The EXCMG method can also be applied to solve other large-scale forward modeling problems encountered in geophysics.

Solving the Poisson equation on small aspect ratio domains using unstructured meshes

This paper describes the development of a robust multigrid, finite element‐based, Laplace solver for accurate capacitance extraction of conductors embedded in multi‐layer dielectric domains. An algebraic multigrid based on element interpolation is adopted and streamlined for the development of the proposed solver. In particular, a new, node‐based agglomeration scheme is proposed to speed up the process of agglomeration. Several attributes of this new method are investigated through the application of the Laplace solver to the calculation of the per‐unit‐length capacitance of configurations of parallel, uniform conductors embedded in multi‐layer dielectric substrates. These two‐dimensional configurations are commonly encountered as high‐speed interconnect structures for integrated electronic circuits. The proposed method is shown to be particularly robust and accurate for structures with very thin dielectric layers characterized by large variation in their electric permittivities. More specifically, it is demonstrated that for such geometries the proposed node‐based agglomeration systematically reduces the problem size and speeds up the iterative solution of the finite element matrix. Copyright © 2007 John Wiley & Sons, Ltd.

Three-precision algebraic multigrid on GPUs

Sparse Matrix Vector multiplication (SpMV) is an important kernel in both traditional high performance computing and emerging data-intensive applications. By far, SpMV libraries are optimized by either application-specific or architecture-specific approaches, making the libraries become too complicated to be used extensively in real applications. In this work we develop a Sparse Matrix-vector multiplication Auto-Tuning system (SMAT) to bridge the gap between specific optimizations and general-purpose usage. SMAT provides users with a unified programming interface in compressed sparse row (CSR) format and automatically determines the optimal format and implementation for any input sparse matrix at runtime. For this purpose, SMAT leverages a learning model, which is generated in an off-line stage by a machine learning method with a training set of more than 2000 matrices from the UF sparse matrix collection, to quickly predict the best combination of the matrix feature parameters. Our experiments show that SMAT achieves impressive performance of up to 51GFLOPS in single-precision and 37GFLOPS in double-precision on mainstream x86 multi-core processors, which are both more than 3 times faster than the Intel MKL library. We also demonstrate its adaptability in an algebraic multigrid solver from Hypre library with above 20% performance improvement reported.

As HPC system sizes grow to millions of cores and chip feature sizes continue to decrease, HPC applications become increasingly exposed to transient hardware faults. These faults can cause aborts and performance degradation. Most importantly, they can corrupt results. Thus, we must evaluate the fault vulnerability of key HPC algorithms to develop cost-effective techniques to improve application resilience.We present an approach that analyzes the vulnerability of applications to faults, systematically reduces it by protecting the most vulnerable components and predicts application vulnerability at large scales. Weinitially focus on sparse scientific applications and apply our approachin this paper to the Algebraic Multi Grid (AMG) algorithm. We empirically analyze AMG's vulnerability to hardware faults in both sequential and parallel (hybrid MPI/OpenMP) executions on up to 1,600 cores and propose and evaluate the use of targeted pointer replication to reduce it. Our techniques increase AMG's resilience to transient hardware faults by 50-80% and improve its scalability on faulty computational environments by 35%. Further, we show how to model AMG's scalability in fault-prone environments to predict execution times of large-scale runs accurately.

Our aim is to develop a physically correct simulation of blood flow through intracranial aneurysms. It shall provide means to estimate rupture risks by calculating the distribution of pressure and shear stresses in an intracranial aneurysm, in order to support the planning of clinical interventions. Due to the time-critical nature of the application, we are forced to use the most efficient state-of-the-art numerical methods and technologies together with high performance computing (HPC) infrastructures. The Navier-Stokes equations for the blood flow are discretized via the finite element method (FEM), and the resulting linear equation systems are handled by an algebraic multigrid (AMG) solver. First comparisons of our simulation results with commercial CFD (computational fluid dynamics) software already show good medical relevance for diagnostic decision support. Another challenge is the visualization of our simulation results at acceptable interaction response rates. Physicians require quick and highly interactive visualization of velocity, pressure and stress to be able to assess the rupture risk of an individual vessel morphology. To meet these demands, parallel visualization techniques and high performance computing resources are utilized. In order to provide physicians with access to remote HPC resources which are not available at every hospital, computing infrastructure of the Austrian Grid is utilized for simulation and visualization.

SUMMARYBootstrap algebraic multigrid (BAMG) is a multigrid‐based solver for matrix equations of the form Ax = b. Its aim is to automatically determine the interpolation weights used in algebraic multigrid by locally fitting a set of test vectors that have been relaxed as solutions to the corresponding homogeneous equation, Ax = 0. This paper studies an improved form of BAMG, called relaxation‐corrected bootstrap algebraic multigrid (rBAMG), that involves adding scaled residuals of the test vectors to the least‐squares equations.The basic rBAMG scheme was introduced in an earlier paper [1] and analyzed on a simple model problem. The purpose of the current paper is to further develop this algorithm by incorporating several new critical components and to systematically study its performance on an interesting model problem from quantum chromodynamics. Whereas the earlier paper introduced a new least‐squares principle involving the residuals of the test vectors, a simple extrapolation scheme is developed here to accurately estimate the convergence factors of the evolving algebraic multigrid solver. Such a capability is essential to the effective development of a fast solver, and the approach introduced here is shown numerically to be much more effective than the conventional approach of just observing successive error reduction factors. Another component of the setup process developed here is an adaptive cycling process. This component assesses the effectiveness of the V‐cycle constructed in the initial rBAMG phase by applying it to the homogeneous equation. When poor convergence is observed, the set of test vectors is enhanced with the resulting error, enabling the subsequent least‐squares fit of interpolation to produce an improved V‐cycle. A related component is the scaling and recombination Ritz process that targets the so‐called weak approximation property in an attempt to reveal the important elements of these evolving error and test vector spaces.The aim of the numerical study documented here is to provide insight into the various design choices that arise in the development of an rBAMG algorithm. With this in mind, the results for quantum chromodynamics focus on the behavior of rBAMG in terms of the number of initial test vectors used, the number of relaxation sweeps applied to them, and the size of the target matrices. Copyright © 2012 John Wiley & Sons, Ltd.

Laplacian matrices of graphs arise in large-scale computational applications such as semisupervised machine learning; spectral clustering of images, genetic data, and web pages; transportation network flows; electrical resistor circuits; and elliptic partial differential equations discretized on unstructured grids with finite elements. A lean algebraic multigrid (LAMG) solver of the symmetric linear system $Ax=b$ is presented, where $A$ is a graph Laplacian. LAMG's run time and storage are empirically demonstrated to scale linearly with the number of edges. LAMG consists of a setup phase, during which a sequence of increasingly coarser Laplacian systems is constructed, and an iterative solve phase using multigrid cycles. General graphs pose algorithmic challenges not encountered in traditional multigrid applications. LAMG combines a lean piecewise-constant interpolation, judicious node aggregation based on a new node proximity measure (the affinity), and an energy correction of coarse-level systems. This results in fast convergence and substantial setup and memory savings. A serial LAMG implementation scaled linearly for a diverse set of 3774 real-world graphs with up to 47 million edges, with no parameter tuning. LAMG was more robust than the UMFPACK direct solver and combinatorial multigrid (CMG), although CMG was faster than LAMG on average. Our methodology is extensible to eigenproblems and other graph computations.

The algebraic multigrid (AMG) method is an efficient solver for linear systems arising in finite element analyses. The AMG method is applicable at a matrix level, different from the geometric multigrid solvers. This paper proposes a combination of the parallel processing technique and the AMG method as a fast solver for electromagnetic field analyses. While the AMG method consists of a setup phase and a solution phase, parallel processing of the former phase is difficult. We present the use of long-range interpolation instead of the conventional direct interpolation for improvement of the parallel efficiency of the AMG setup phase. A magnetostatic analysis and an eddy-current analysis show the solver performance. The numerical results show that parallelized AMG is a fast solver and has sufficient scalability, as compared with the conventional solver.

Iterative linear equation solvers depend on high quality preconditioners to achieve fast convergence. For sparse symmetric systems arising from large power grid analysis problems, however, preconditioners generated by traditional incomplete Cholesky factorizations are usually of low quality, resulting in slow convergence. On the other hand, although preconditioners generated by random walks are quite effective to reduce the number of iterations, it takes considerable amount of time to compute them in a stochastic manner. In this paper, we propose a new preconditioning technique for power grid analysis, named deterministic random walk, that combines the advantages of the above two approaches. Our proposed algorithm computes the preconditioners in a deterministic manner to reduce computation time, while achieving similar quality as stochastic random walk preconditioning by modifying fill-ins to compensate dropped entries. We have proved that for such compensation scheme, our algorithm will not fail for certain matrix orderings, which otherwise cannot be guaranteed by traditional incomplete factorizations. We demonstrate that by incorporating our proposed preconditioner, a conjugate gradient solver is able to outperform a state-of-the-art algebraic multigrid preconditioned solver on public IBM power grid benchmarks for DC power grid analysis, while potentially remaining very efficient for transient analysis.

A node-based algebraic multi-grid finite element Laplace solver is demonstrated for the extraction of the multi-terminal DC resistance of the power distribution network at the board and package level. The solver can handle boards of arbitrary complexity, comprehending without numerical difficulties fine features associated with dense clusters of via holes and related voids in the metallization.