Efficient parallel multigrid based solvers for large scale groundwater flow simulations

In this report we present parallel solvers for large linear systems arising from the finite-element discretization of the three-dimensional steady-state groundwater flow problem. Our solvers are based on multigrid and Krylov subspace methods. The parallel implementation is based on a domain decomposition strategy with explicit message passing using NX and MPI libraries. We have tested our parallel implementations on the Intel Paragon XP/S 150 supercomputer using up to 1024 parallel processors and on other parallel platforms such as SGI/Power Challenge Array, Cray/SGI Origin 2000, Convex Exemplar SPP-1200, and IBM SP using up to 64 processors. We show that multigrid can be a scalable algorithm on distributed memory machines. We demonstrate the effectiveness of parallel multigrid based solvers by solving problems requiring more than 70 million nodes in less than a minute. This is more than 25 times faster than the diagonal preconditioned conjugate gradient method which is one of the more popular methods for large sparse linear systems. Our results also show that multigrid as a stand alone solver works best for problems with smooth coefficients, but for rough coefficients it is best used as a preconditioner for a Krylov subspace method such as the conjugate gradient method. We show that even for extremely heterogeneous systems the multigrid pre-conditioned conjugate gradient method is at least 10 times faster than the diagonally preconditioned conjugate gradient method.

Sparsified Randomization algorithms for low rank approximations and applications to integral equations and inhomogeneous random field simulation

The Newton-based global gradient algorithm (GGA) (also known as the Todini and Pilati method) is a widely used method for computing the steady-state solution of the hydraulic variables within a water distribution system (WDS). The Newton-based computation involves solving a linear system of equations arising from the Jacobian of the WDS equations. This step is the most computationally expensive process within the GGA, particularly for large networks involving up to O(105) variables. An increasingly popular solver for large linear systems of the M-matrix class is the algebraic multigrid (AMG) method, a hierarchical-based method that uses a sequence of smaller dimensional systems to approximate the original system. This paper studies the application of AMG to the steady-state solution of WDSs through its incorporation as the linear solver within the GGA. The form of the Jacobian within the GGA is proved to be an M-matrix (under specific criteria on the pipe resistance functions), and thus able to be...

Groundwater flow simulation has become one of the top international issues in new generation of environmental applications. When managing large-scale groundwater flow problems, the intensive computational ability and large amounts of memory space required for modeling are the main bottlenecks for researchers. In order to solve three-dimensional large-scale groundwater flow problems more rapidly, the Open MP was adopted to parallelize the preconditioned conjugate gradient (PCG) algorithm in this paper. And this paper carried out a numerical experiment of the three-dimensional groundwater flow model on a computer with four cores. Based on the numerical experiment, it is found that the execution time of the original serial PCG program is about 1.74 to 2.86 times of the parallel PCG program executed with different number of threads. The experimental results also demonstrate that the PCG solver based on Open MP is an effective way for solving large-scale three-dimensional groundwater flow problem.

This article presents a parallel simulation solver for groundwater flow on CUDA. Preconditioned conjugate gradient (PCG) algorithm is used to solve the large linear systems arising from the finite-difference discretization of three-dimensional groundwater flow problems. CUDA implementing methods for the two most time-consuming operations in PCG, sparse matrix–vector multiplication and vector inner-product, are given. The experimental results show that CUDA can speed up the solving process of the groundwater simulation significantly. 1.8–3.7 speedup can be achieved with different GPUs for a transient groundwater flow problem.

A high-throughput hybrid task and data parallel Poisson solver for large-scale simulations of incompressible turbulent flows on distributed GPUs

The solution of large eigenproblems is involved in many scientific and engineering applications when for instance, stability analysis is a concern. For large simulation in material physics or thermo-acoustics, the calculation can last for many hours on large parallel platforms. On future large-scale systems, the mean time between failures (MTBF) of the system is expected to decrease so that many faults could occur during the solution of large eigenproblems. Consequently, it becomes critical to design parallel eigensolvers that can survive faults. In that framework, we investigate the relevance of approaches relying on numerical techniques, which might be combined with more classical techniques for real large-scale parallel implementations. Because we focus on numerical remedies we do not consider parallel implementations nor parallel experiments but only numerical experiments. We assume that a separate mechanism ensures the fault detection and that a system layer provides support for setting back the environment (processes,. . .) in a running state. Once the system is in a running state, after a fault, our main objective is to provide robust resilient schemes so that the eigensolver may keep converging in the presence of the fault without restarting the calculation from scratch. For this purpose, we extend the interpolation-restart (IR) strategies initially introduced for the solution of linear systems in a previous work to the solution of eigenproblems in this paper. For a given numerical scheme, the IR strategies consist of extracting relevant spectral information from available data after a fault. After data extraction, a well-selected part of the missing data is regenerated through interpolation strategies to constitute a meaningful input to restart the numerical algorithm. One of the main features of this numerical remedy is that it does not require extra resources, i.e., computational unit or computing time, when no fault occurs. In this paper, we revisit a few state-of-the-art methods for solving large sparse eigenvalue problems namely the Arnoldi methods, subspace iteration methods and the Jacobi-Davidson method, in the light of our IR strategies. For each considered eigensolver, we adapt the IR strategies to regenerate as much spectral information as possible. Through extensive numerical experiments, we study the respective robustness of the resulting resilient schemes with respect to the MTBF and to the amount of data loss via qualitative and quantitative illustrations. 1. Introduction. The computation of eigenpairs (eigenvalues and eigenvectors) of large sparse matrices is involved in many scientific and engineering applications such as when stability analysis is a concern. To name a few, it appears in structural dynamics, thermodynamics, thermo-acoustics, quantum chemistry. With the permanent increase of the computational power of high performance computing (HPC) systems by using a larger and larger number of CPU cores or specialized processing units, HPC applications are increasingly prone to faults. To guarantee fault tolerance, two classes of strategies are required. One for the fault detection and the other for fault correction. Faults such as computational node crashes are obvious to detect while silent faults may be challenging to detect. To cope with silent faults, a duplication strategy is commonly used for fault detection [18, 39] by comparing the outputs, while triple modular redundancy (TMR) is used for fault detection and correction [34, 37]. However, the additional computational resources required by such replication strategies may represent a severe penalty. Instead of replicating computational resources, studies [7, 36] propose a time redundancy model for fault detection. It consists in repeating computation twice on the same resource. The advantage of time redundancy models is the flexibility at application level; software developers can indeed select only a set of critical instructions to protect. Recomputing only some instructions instead of the whole application lowers the time redundancy overhead [25]. In some numerical simulations, data naturally satisfy well defined mathematical properties. These properties can be efficiently exploited for fault detection through a periodical check of the numerical properties during computation [10]. Checkpoint/restart is the most studied fault recovery strategy in the context of HPC systems. The common checkpoint/restart scheme consists in periodically saving data onto a reliable storage device such as a remote disk. When a fault occurs, a rollback is performed to the point of the most recent and consistent checkpoint. According to the implemented checkpoint strategy, all processes

Modern petroleum reservoir simulation serves as a primary tool for quantitatively managing reservoir production and planning new fields. It involves repeatedly solving the Jacobian of a set of strong nonlinear partial differential equations governing the mass and energy conduction and conservation. Most of the existing reservoir simulators adopt iterative solver with multiple stages of preconditioners, in which the incomplete LU (ILU) factorization is an outstanding universal smoother. However, it turns out that when the degree of freedom of each grid grows, ILU usually becomes the bottleneck of the solver. Moreover, ILU is difficult to parallelize due to its inherent data dependency. In this paper, we developed a sparse iterative solver with parallelized ILU and triangular solve using block-wise data structure. Compared with the state of art iterative solver on 14 industrial reservoir simulation matrices, the proposed ILU is 5.2x faster (on average) than the state of art iterative solver because of the block-wise data structure, which leads to 2.2x speedup on the total solver runtime. In addition, parallel ILU and triangular solve are developed to further accelerate the solver. To tackle the strong data dependency in ILU and triangular solve, we first partition the algorithm into separated tasks and construct a data flow graph to represent the data dependency. Then, tasks are scheduled in parallel according to the topological order of the data flow graph. On an 8-thread multicore architecture, we achieved another 3.6x speedup on ILU factorization, and 3.3x on triangular solve with good scalability.

The power of high-performance computing (HPC) is applied to simulate highly complex systems and processes in many scientific communities, e. g. in particle physics, weather and climate research, bio-sciences, materials science, pharmaceutics, astronomy, or finance.Current HPC systems are so complex that the design of such a system, including architecture design space exploration and performance prediction, requires HPC-like simulation capabilities. To this end, we developed an Omnest-based simulation environment that enables studying the impact of an HPC machine's communication subsystem on the overall system's performance for specific workloads.As the scale of current high-end HPC systems is in the range of hundreds of thousands of processing cores, full system simulation---at an abstraction level that still maintains a reasonably high level of detail---is infeasible without resorting to parallel simulation, the main limiting factors being simulation run time and memory footprint.We describe our experiences in adapting our simulation environment to take advantage of the parallel distributed simulation capabilities provided by Omnest. We present results obtained on a many-core SMP machine as well as a small-scale InfiniBand cluster.Furthermore, we ported our simulation environment, including Omnest itself, to the massively parallel IBM Blue Gene®/P platform. We report results from initial experiments on this platform using up to 512 cores in parallel.

We investigate and compare stable parallel algorithms for solving diagonally dominant and general narrow-banded linear systems of equations. Narrow-banded means that the bandwidth is very small compared with the matrix order and is typically between 1 and 100. The solvers compared are the banded system solvers of ScaLAPACK [12] and those investigated by Arbenz and Hegland [4, 8]. For the diagonally dominant case, the algorithms are analogs of the well-known tridiagonal cyclic reduction algorithm, while the inspiration for the general case is the lesser-known bidiagonal cyclic reduction, which allows a clean parallel implementation of partial pivoting. These divide-and-conquer type algorithms complement fine-grained algorithms which perform well only for wide-banded matrices, with each family of algorithms having a range of problem sizes for which it is superior. We present theoretical analyses as well as numerical experiments conducted on the Intel Paragon.

Graph algorithms are notorious for not getting good speedup on parallel architectures. These algorithms tend to suffer from irregular dependencies and a high synchronization cost that prevent an efficient execution on distributed memory machines. Hence such algorithms are mostly parallelized on shared memory machines. However, current commodity shared memory machines do not typically offer enough parallelism to process these problems. In this paper, we are presenting an early investigation of the scalability of such algorithms on Intel's upcoming Many Integrated Core (Intel MIC) architecture which, when it will be released in 2012, is expected to provide more than 50 physical cores with SMT capability. The Intel MIC architecture can be programmed through many programming models, here we investigate the three most popular of these models namely OpenMP, Cilk Plus and Intel's TBB. We present scalability results of a parallel graph coloring algorithm, three variations of a breadth-first search algorithm and a micro benchmark for irregular computations using these three programming models. Our results on a prototype board show that the multi-threaded architecture of Intel MIC can be effectively used for hiding latencies in irregular applications to achieve almost perfect speedup.

A large number of numerical algorithms exhibit a two-level structure with both method parallelism and system parallelism. This structure can be exploited to produce alternative parallel implementations on distributed memory machines. The compiler system TwoL, (Two Level) provides interactive and semiautomatic support for the design and realization of efficient parallel algorithms in this two-level parallel programming model. The design is structured into well-defined decision steps which are formalized in a TwoL specification language, and transformations on this language. We show how the design steps lead to a parallel algorithm, how the design is formalized in the TwoL system, how this compiler system is realized, and which algorithms are amenable to automated decision steps. Design or derivation steps are based on parameterized cost functions arising from runtime predictions for the specific parallel target machine. The design process is illustrated by the parallelization of several methods from the area of differential equations.

This paper describes an efficient implementation of a load balancing algorithm based upon a process migration mechanism, and the notion of CPU-queue length which indicates the workload at each processor. The main objective is to reduce the message overhead in conservative simulation. The algorithms were tested on four simulation workload models to study the the scalability of our algorithm on an Intel Paragon multicomputer architecture using 64 processors. A significant decrease in run time and the message overhead associated with the conservative simulation was obtained with the use of the proposed dynamic load balancing algorithm as compared to the use of a static partitioning algorithm.

Scalable domain decomposition solvers for stochastic PDEs in high performance computing

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 ((Pk)d,P k‐1disc) 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 \documentclass{article}\usepackage{mathrsfs}\usepackage{amsmath, amssymb}\pagestyle{empty}\begin{document}\begin{align*}\mathcal{H}\end{align*} \end{document} ‐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