Parallelization of the GESIMA mesoscale atmospheric model

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 purpose of the article is to present and evaluate usefulness of a new preconditioning technique (for the Gauss-Seidel algorithm), namely incomplete WZ factorization, for iterative solving of sparse and singular linear equations systems, which arise during modeling with Markov chains. The incomplete WZ factorization proposed here will be compared with the incomplete LU factorization in respect of amount of fill-in (newly created non-zeros) and in respect of the accuracy improvement of preconditioned algorithms in relation to not preconditioned ones. In the paper, the results of some numerical experiments will be presented, which were conducted for various matrices representing Markov chains. The experiments show that the incomplete WZ factorization can be a real alternative - because it is faster than incomplete LU factorization and the fill-in generated in the process is smaller (the output matrices are sparser).

When solving large linear systems of equations arising from the discretization of elliptic boundary value problems, a combination of iterative methods and preconditioners based on incomplete LU factorizations is frequently used. Given a model problem with variable coefficients, we investigate a class of incomplete LU factorizations depending on a relaxation parameter. We show that the associated preconditioner and the factorization itself both are numerically stable. The theoretical results are complemented by numerical experiments.

A highly efficient and robust scheme is proposed for analyzing electromagnetic scattering from electrically large arbitrary shaped conductors in a half space. This scheme is based on the electric field integral equation (EFIE) with a half-space Green’s function. The precorrected fast Fourier transform (p-FFT) is first extended to a half space for general three-dimensional scattering problems. A novel enhanced dual threshold incomplete LU factorization (ILUT) is then constructed as an effective preconditioner to improve the convergence of the half-space EFIE. Inspired by the idea of the improved electric field integral operator (IEFIO), the geometrical-optics current/principle value term of the magnetic field integral equation is used as a physical perturbation to stabilize the traditional ILUT perconditioning matrix. The high accuracy of EFIE is maintained, yet good calculating efficiency comparable to the combined field integral equation (CFIE) can be achieved. Furthermore, this approach can be applied to arbitrary geometrical structures including open surfaces and requires no extra types of Sommerfeld integrals needed in the half-space CFIE. Numerical examples are presented to demonstrate the high performance of the proposed solver among several other approaches in typical half-space problems.

We present fast and highly scalable parallel computations for a number of important and fundamental matrix problems on distributed memory systems (DMS). These problems include matrix multiplication, matrix chain product, and computing the powers, the inverse, the characteristic polynomial, the determinant, the rank, the Krylov matrix, and an LU- and a QR-factorization of a matrix, and solving linear systems of equations. Our highly scalable parallel computations for these problems are based on a highly scalable implementation of the fastest sequential matrix multiplication algorithm on DMS. We show that compared with the best known parallel time complexities on parallel random access machines (PRAM), the most powerful but unrealistic shared memory model of parallel computing, our parallel matrix computations achieve the same speeds on distributed memory parallel computers (DMPC), and have an extra polylog factor in the time complexities on DMS with hypercubic networks. Furthermore, our parallel matrix computations are fully scalable on DMPC and highly scalable over a wide range of system size on DMS with hypercubic networks. Such fast (in terms of parallel time complexity) and highly scalable (in terms of our definition of scalability) parallel matrix computations were rarely seen before on any distributed memory systems.

We present fast and scalable parallel computations for a number of important and fundamental matrix problems on distributed memory systems (DMS). These problems include computing the powers, the inverse, the characteristic polynomial, the determinant, the rank, the Krylov matrix, and an LU- and a QR-factorization of a matrix, and solving linear systems of equations. These parallel computations are based on efficient implementations of the fastest sequential matrix multiplication algorithm on DMS. We show that compared with the best known time complexities on PRAM, our parallel matrix computations achieve the same speeds on distributed memory parallel computers (DMPC), and have an extra polylog factor in the time complexities on DMS with hypercubic networks. Furthermore, our parallel matrix computations are fully scalable on DMPC and highly scalable over a wide range of system size on DMS with hypercubic networks. Such fast and scalable parallel matrix computations were not seen before on any distributed memory systems.

Application demands and grand challenges in numerical simulation require for both highly capable computing platforms and efficient numerical solution schemes. Power constraints and further miniaturization of modern and future hardware give way for multi- and manycore processors with increasing fine-grained parallelism and deeply nested hierarchical memory systems -- as already exemplified by recent graphics processing units. Accordingly, numerical schemes need to be adapted and re-engineered in order to deliver scalable solutions across diverse processor configurations. Portability of parallel software solutions across emerging hardware platforms is another challenge. This work investigates multi-coloring and re-ordering schemes for block Gauß-Seidel methods and, in particular, for incomplete LU factorizations with and without fill-ins. We consider two matrix re-ordering schemes that deliver flexible and efficient parallel preconditioners. The general idea is to generate block decompositions of the system matrix such that the diagonal blocks are diagonal itself. In such a way, parallelism can be exploited on the block-level in a scalable manner. Our goal is to provide widely applicable, out-of-the-box preconditioners that can be used in the context of finite element solvers. We propose a new method for anticipating the fill-in pattern of ILU($p$) schemes which we call the power($q$)-pattern method. This method is based on an incomplete factorization of the system matrix $A$ subject to a predetermined pattern given by the matrix power $|A|^(p+1)$ and its associated multi-coloring permutation $. We prove that the obtained sparsity pattern is a superset of our modified ILU($p$) factorization applied to pi A pi^(-1). As a result, this modified ILU($p$) applied to multi-colored system matrix has no fill-ins in its diagonal blocks. This leads to an inherently parallel execution of triangular ILU($p$) sweeps. In addition, we describe the integration of the preconditioners into the HiFlow$^3$ open-source finite element package that provides a portable software solution across diverse hardware platforms. On this basis, we conduct performance analysis across a variety of test problems on multi-core CPUs and GPUs that proves efficiency, scalability and flexibility of our approach. Our preconditioners achieve a solver acceleration by a factor of up to 1.5, 8 and 85 for three different test problems. The GPU versions of the preconditioned solver are by a factor of up to 4 faster than an OpenMP parallel version on eight cores.

Program for quantum wave-packet dynamics with time-dependent potentials

Currently, there are two approaches to implementing large-scale ice sheet simulations. One approach, which is based on the use of an evenly-spaced structured grid, is to simulate the entire ice sheet at the highest resolution required. The other approach is to model the ice sheet at different levels of resolution using unstructured meshes. Structured grids lend themselves quite well to partitioning on multi-core, distributed memory parallel computers, leading to good computational efficiency, load-balancing and reduced communication costs. The disadvantage is the significant computational resources required to model a system of this size at a high enough resolution to correctly model regions of the ice sheet undergoing rapid change. Unstructured meshes can reduce such costs by allocating more computational resources to areas of the ice sheet whose dynamics are evolving quickly and fewer resources to areas evolving more slowly. However, unstructured meshes do not easily lend themselves to effective decomposition on distributed memory parallel computers, and can experience poor cache utilization, significant load imbalance, high communication costs, and increased synchronization times. In this paper, we describe our alternative approach based on embedded simulation, which provides a degree of multi-resolution capabilities to structured grids while maintaining many of the computational efficiencies that come with the structured approach. We provide experimental results focusing on the scalability characteristics and cache efficiency of each approach as a function of increasing problem size and core counts. We show that the embedded approach scales very well across all problem sizes and core counts, and that the unstructured approach scales well at smaller problem sizes and core counts but degrades significantly as the problem size becomes large.

This paper presents the results of the Connection Machine implementation of a number of preconditioners for the preconditioned conjugate gradient method. The preconditioners implemented include those based on the incomplete LU factorization, the modified incomplete LU factorization, the symmetric successive overrelaxation, and others such as several polynomial preconditioners and the hierarchical basis preconditioner. Results based on numerical experiments show that both the degree of parallelism inherent in a preconditioner and its convergence rate improvement play important roles on the overall execution time performance on parallel computers. Factors that affect the performance of the preconditioners will also be discussed. We conclude that to search for the best preconditioner on a parallel machine, we have to consider the tradeoffs between fast convergence rate and high degree of parallelism as well as the architecture of the target parallel computer.

Parallel computation of three-dimensional flows using overlapping grids with adaptive mesh refinement

Chapter 9 - Conjugate Gradient-Type Methods

The aim of this work is to build a numerical software library based on the \(\tau \)-method to solve differential problems using MATLAB. The \(\tau \)-method can be very effective in the solution of certain type of these problems, and therefore, the existence of a numerical library for its dissemination is of major importance. Furthermore, the method has been used for the solution of particular problems but has not yet been explored as a general technique. Focus will be on stability issues, namely those issued from the solution of algebraic linear systems required for the process. Additionally, preconditioners for the solution with the \(\tau \)-method will be tackled, with emphasizes on incomplete LU factorizations and (block) Jacobi preconditioners. We also propose an iterative approach, build upon an LU factorization over a moderate initial size, generating better approximations and providing a priori error estimate at each iteration. Numerical results enlightening the efficiency of the proposed methods will be presented.

A parallel preconditioned conjugate gradient package for solving sparse linear systems on a Cray Y-MP

An efficient parallel-computing method for modeling nonisothermal multiphase flow and multicomponent transport in porous and fractured media