Parallel Direct Solver Research Articles

On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver

ABSTRACTWe consider the modeling of (polarized) seismic wave propagation on a rectangular domain via the discretization and solution of the inhomogeneous Helmholtz equation in 3D, by exploiting a parallel multifrontal sparse direct solver equipped with Hierarchically Semi‐Separable (HSS) structure to reduce the computational complexity and storage. In particular, we are concerned with solving this equation on a large domain, for a large number of different forcing terms in the context of seismic problems in general, and modeling in particular. We resort to a parsimonious mixed grid finite differences scheme for discretizing the Helmholtz operator and Perfect Matched Layer boundaries, resulting in a non‐Hermitian matrix. We make use of a nested dissection based domain decomposition, and introduce an approximate direct solver by developing a parallel HSS matrix compression, factorization, and solution approach. We cast our massive parallelization in the framework of the multifrontal method. The assembly tree is partitioned into local trees and a global tree. The local trees are eliminated independently in each processor, while the global tree is eliminated through massive communication. The solver for the inhomogeneous equation is a parallel hybrid between multifrontal and HSS structure. The computational complexity associated with the factorization is almost linear in the size, n say, of the matrix, viz. between O(n log n) and O(n4/3 log n), while the storage is almost linear as well, between O(n) and O(n log n). We exploit the use of a regular (Cartesian) mesh common in many seismic applications.

FANTOM: Two- and three-dimensional numerical modelling of creeping flows for the solution of geological problems

A new finite element code for the solution of the Stokes and heat transport equations is presented. It has purposely been designed to address geological flow problems in two and three dimensions at crustal and lithospheric scales. A variety of rheologies has been implemented including nonlinear thermally activated creep and brittle (or plastic) frictional models. A cloud of particles is used to track materials in the simulation domain which allows to record the integrated history of deformation; its density is variable and dynamically adapted. The code is built on the Arbitrary Lagrangian–Eulerian kinematical description: the computational grid deforms vertically and allows for a true free surface while the computational domain remains of constant width in the horizontal direction. The code can be run in sequential or parallel mode. The parallelisation is based on the MPI paradigm and the domain decomposition algorithm is presented. The solution to the large system of algebraic equations resulting from the finite element discretisation and linearisation of the set of coupled partial differential equations to be solved is obtained by means of an efficient (sequential or massively parallel) direct solver. Details of implementation concerning plasticity, cloud handing, nonlinear convergence, strain accumulation and pressure smoothing are given. The sequential 2D version of the code is used to run the numerical sandbox benchmark at normal and high resolutions. The 2D parallel version is used in the case of a thermo-mechanically coupled extension experiment in which the mantle is present. The 3D parallel version is used to run a crustal scale orogeny experiment. The overall performance of the code, as well as the respect of the incompressibility constraint, of the yield stress criterion and of the volume conservation are discussed, along with parallel scalability. Benchmark results of scalar field advection, the Rayleigh–Taylor experiment and the falling block experiment are presented in Appendix.

A parallel block LU decomposition method for distributed finite element matrices

In this work we present a new parallel direct linear solver for matrices resulting from finite element problems. The algorithm follows the nested dissection approach, where the resulting Schur complements are also distributed in parallel. The sparsity structure of the finite element matrices is used to pre-compute an efficient block structure for the LU factors. We demonstrate the performance and the parallel scaling behavior by several test examples.

A new parallel sparse direct solver: Presentation and numerical experiments in large‐scale structural mechanics parallel computing

AbstractThe main purpose of this work is to present a new parallel direct solver: Dissection solver. It is based on LU factorization of the sparse matrix of the linear system and allows to detect automatically and handle properly the zero‐energy modes, which are important when dealing with DDM. A performance evaluation and comparisons with other direct solvers (MUMPS, DSCPACK) are also given for both sequential and parallel computations. Results of numerical experiments with a two‐level parallelization of large‐scale structural analysis problems are also presented: FETI is used for the global problem parallelization and Dissection for the local multithreading. In this framework, the largest problem we have solved is of an elastic solid composed of 400 subdomains running on 400 computation nodes (3200 cores) and containing about 165 millions dof. The computation of one single iteration consumes less than 20 min of CPU time. Several comparisons to MUMPS are given for the numerical computation of large‐scale linear systems on a massively parallel cluster: performances and weaknesses of this new solver are highlighted. Copyright © 2011 John Wiley & Sons, Ltd.

Goal-Oriented Self-Adaptive hp Finite Element Simulation of 3D DC Borehole Resistivity Simulations

Abstract In this paper we present a goal-oriented self-adaptive hp Finite Element Method ( hp -FEM ) with shared data structures and a parallel multi-frontal direct solver. The algorithm automatically generates (without any user interaction) a sequence of meshes delivering exponential convergence of a prescribed quantity of interest with respect to the number of degrees of freedom. The sequence of meshes is generated from a given initial mesh, by performing h (breaking elements into smaller elements), p (adjusting polynomial orders of approximation) or hp (both) refinements on the finite elements. The new parallel implementation utilizes a computational mesh shared between multiple processors. All computational algorithms, including automatic hp goal-oriented adaptivity and the solver work fully in parallel. We describe the parallel self-adaptive hp- FEM algorithm with shared computational domain, as well as its efficiency measurements. We apply the methodology described to the three-dimensional simulation of the borehole resistivity measurement of direct current through casing in the presence of invasion.

Two-dimensional frequency-domain visco-elastic full waveform inversion: Parallel algorithms, optimization and performance

Full waveform inversion (FWI) is an appealing seismic data-fitting procedure for the derivation of high-resolution quantitative models of the subsurface at various scales. Full modelling and inversion of visco-elastic waves from multiple seismic sources allow for the recovering of different physical parameters, although they remain computationally challenging tasks. An efficient massively parallel, frequency-domain FWI algorithm is implemented here on large-scale distributed-memory platforms for imaging two-dimensional visco-elastic media. The resolution of the elastodynamic equations, as the forward problem of the inversion, is performed in the frequency domain on unstructured triangular meshes, using a low-order finite element discontinuous Galerkin method. The linear system resulting from discretization of the forward problem is solved with a parallel direct solver. The inverse problem, which is presented as a non-linear local optimization problem, is solved in parallel with a quasi-Newton method, and this allows for reliable estimation of multiple classes of visco-elastic parameters. Two levels of parallelism are implemented in the algorithm, based on message passing interfaces and multi-threading, for optimal use of computational time and the core-memory resources available on modern distributed-memory multi-core computational platforms. The algorithm allows for imaging of realistic targets at various scales, ranging from near-surface geotechnic applications to crustal-scale exploration.

Parallel Direct Solver for the Finite Integration Technique in Electrokinetic Problems

The finite integration technique allows the simulation of real-world electromagnetic field problems with complex geometries. It provides a discrete reformulation of Maxwell's equations in their integral form suitable for numerical computing. The resulting matrix equations of the discretized fields can be used for efficient numerical simulations on modern computers and can be exploited to use a parallel computing. In fact, by reordering the unknowns by the nested dissection method, it is possible to directly construct the lower triangular matrix of the Cholesky factorization with many processors without assembling the matrix system. In this paper, a parallel algorithm is proposed for the direct solution of large sparse linear systems with the finite integration technique. This direct solver has the advantage of handling singularities in the matrix of linear systems. The computational effort for these linear systems, often encountered in numerical simulation of electromagnetic phenomena by finite integration technique, is very significant in terms of run-time and memory requirements. Many numerical tests have been carried out to evaluate the performance of the parallel direct solver.

Parallel direct solver for the finite integration technique in electromagnetics

PurposeThe purpose of this paper is to develop a solving method based on the parallelization of the calculus of the finite integration technique on many processors and discuss its efficiency following the number of processors.Design/methodology/approachThe finite integration technique is used in different engineering domains as well as to compute the electromagnetic phenomena. This technique is efficient and generates a linear system with regular structure, this system can be implemented and solved in parallel with a direct solver. In fact, in reordering the unknowns by the nested dissection method, the lower triangular matrix of the Cholesky factorization can be constructed with many processors without assembling the matrix system.FindingsThis paper deals with the parallelization of the finite integration technique which is performed by data partitioning and leads to a high‐performance.Originality/valueThe paper presents a parallel implementation of the finite integration technique associated to a direct solver which is practical and efficient.

Parallel multi-frontal solver for multi-physics p adaptive problems

The paper presents a parallel direct solver for multi-physics problems. The solver is dedicated for solving problems resulting from adaptive Finite Element Method computations. The concept of finite element is actually replaced by the concept of the node. The computational mesh consists of several nodes, related to element vertices, edges, faces and interiors. The ordering of unknowns in the solver is performed on the level of nodes. The concept of the node can be efficiently utilized in order to recognize unknowns that can be eliminated at a given node of the elimination tree. The solver is tested on the exemplary three dimensional multi-physics problem involving the computations of the linear acoustics coupled with linear elasticity. The three dimensional tetrahedral mesh generation and the solver algorithm are modeled by using graph grammar formalism. The execution time and the memory usage of the solver are compared with the MUMPS solver.

Parallel multi-frontal solver for p adaptive finite element modeling of multi-physics computational problems

The paper presents a parallel direct solver for multi-physics problems. The solver is dedicated for solving problems resulting from adaptive finite element method computations. The concept of finite element is actually replaced by the concept of the node. The computational mesh consists of several nodes, related to element vertices, edges, faces and interiors. The ordering of unknowns in the solver is performed on the level of nodes. The concept of the node can be efficiently utilized in order to recognize unknowns that can be eliminated at a given node of the elimination tree. The solver is tested on the exemplary three-dimensional multi-physics problem involving the computations of the linear acoustics coupled with linear elasticity. The three-dimensional tetrahedral mesh generation and the solver algorithm are modeled by using graph grammar formalism. The execution time and the memory usage of the solver are compared with the MUMPS solver.

Graph grammar-based multi-thread multi-frontal parallel solver with trace theory-based scheduler

The paper presents the graph grammar based multi-thread multi-frontal parallel direct solver for one and two dimensional Finite Difference Method (FDM). The multi-frontal solver algorithm has been expressed by graph grammar productions. Each production represents an atomic task that internally must be executed in serial. The sequence of graph grammar productions modeling the execution of the solver has been associated with the alphabet for the trace theory analysis. The dependency relation between tasks has been introduced based on the analysis of the solver algorithm. The sequence of productions has been transformed into the Foata Normal Form (FNF). The parallel solver algorithm has been implemented and tested on NVIDIA Cuda multi-core graphic card. The tasks have been scheduled according to the classes in the FNF.

Direct Vlasov Solvers with High-Order Spectral Element Method

This paper presents the development of parallel direct Vlasov solvers using the Spectral Element Method (SEM). Instead of the standard Particle-In-Cell (PIC) approach for kinetic space plasma simulation, i.e. solving the Vlasov-Maxwell equations, the direct method has been used in this paper. There are several benefits to solve the Vlasov equation directly, such as avoiding noise associated with the finite number of particles and the capability to capture the fine structure in the plasma, etc. The most challenging part of direct Vlasov solver comes from high dimension, as the computational cost increases as N{sup 2d}, where d is the dimension of the physical space. Recently due to fast development of supercomputers, the possibility of high dimensions becomes more realistic. A significant effort has been devoted to solve the Vlasov equation in low dimensions so far, now more interests focus on higher dimensions. Different numerical methods have been tried so far, such as finite difference method, Fourier spectral method, finite volume method, etc. In this paper SEM has been successfully applied to construct these solvers. SEM has shown several advantages, such as easy interpolation due to local element structure and long time integration due to its high order accuracy. Domain decompositionmore » in high dimensions have been used for parallelization, these include scalable parallel 1D and 2D Poisson solvers. Benchmark results have been shown and simulation results have been reported for two different cases: one dimension (1P1V), and two dimensions (2P2V) in both physical and velocity spaces.« less

409 ADVENTURE SolidのBDD前処理への並列直接ソルバーの導入(OS4. 大規模並列・連成解析と関連話題(2),オーガナイズドセッション講演)

409 ADVENTURE SolidのBDD前処理への並列直接ソルバーの導入(OS4. 大規模並列・連成解析と関連話題(2),オーガナイズドセッション講演)

A Parallel Direct Solver for the Simulation of Large-Scale Power/Ground Networks

An algorithm is presented for the fast and accurate simulation of power/ground mesh structures. Our method is a direct (non-iterative) approach for simulation based upon a parallel matrix inversion algorithm. The new dimension of flexibility provided by our algorithm allows for a more accurate analysis of power/ground mesh structures using resistance, inductance, capacitance, interconnect models. Specifically, we offer a method that employs a sparse approximate inverse technique to consider more reluctance coupling terms for increased accuracy of simulation. Our algorithm shows substantial computational improvement over the best known direct and iterative numerical techniques that are applicable to these large-scale simulation problems.

Graph grammar‐driven parallel partial differential equation solver

AbstractThe paper presents an extension of the composite programmable graph grammar (CP‐graph grammar) suitable for modeling the parallel direct solver algorithm utilized by the hp finite element method (hp‐FEM). In the proposed graph grammar model, the computational mesh is represented by a CP‐graph. The presented graph grammar models the solver algorithm by a set of graph grammar productions. The graph grammar model makes it possible to examine the concurrency of the algorithm by analyzing the interdependence between the atomic tasks, tasks and super‐tasks. The atomic tasks correspond to the graph grammar productions, representing basic undividable parts of the algorithms. The level of atomic tasks models the concurrency for the shared memory architectures. On the other hand, the tasks correspond to the groups of atomic tasks with predefined inter‐task communication channels. They constitute the grain for the decomposition of the parallel algorithm for the distributed memory architecture. Finally, the super‐tasks correspond to a group of tasks resulting from the execution of load balancing algorithm. The solver algorithm is tested on distributed memory linux cluster for up to 192 processors. Copyright © 2009 John Wiley & Sons, Ltd.

A parallel direct solver for the self-adaptive hp Finite Element Method

In this paper we present a new parallel multi-frontal direct solver, dedicated for the hp Finite Element Method (hp-FEM). The self-adaptive hp-FEM generates in a fully automatic mode, a sequence of hp-meshes delivering exponential convergence of the error with respect to the number of degrees of freedom (d.o.f.) as well as the CPU time, by performing a sequence of hp refinements starting from an arbitrary initial mesh. The solver constructs an initial elimination tree for an arbitrary initial mesh, and expands the elimination tree each time the mesh is refined. This allows us to keep track of the order of elimination for the solver. The solver also minimizes the memory usage, by de-allocating partial LU factorizations computed during the elimination stage of the solver, and recomputes them for the backward substitution stage, by utilizing only about 10% of the computational time necessary for the original computations. The solver has been tested on 3D Direct Current (DC) borehole resistivity measurement simulations problems. We measure the execution time and memory usage of the solver over a large regular mesh with 1.5 million degrees of freedom as well as on the highly non-regular mesh, generated by the self-adaptive h p -FEM, with finite elements of various sizes and polynomial orders of approximation varying from p = 1 to p = 9 . From the presented experiments it follows that the parallel solver scales well up to the maximum number of utilized processors. The limit for the solver scalability is the maximum sequential part of the algorithm: the computations of the partial LU factorizations over the longest path, coming from the root of the elimination tree down to the deepest leaf.

Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic (TTI) media

A 2D finite-difference, frequency-domain method was developed for modeling viscoacoustic seismic waves in transversely isotropic media with a tilted symmetry axis. The medium is parameterized by the P-wave velocity on the symmetry axis, the density, the attenuation factor, Thomsen’s anisotropic parameters [Formula: see text] and [Formula: see text], and the tilt angle. The finite-difference discretization relies on a parsimonious mixed-grid approach that designs accurate yet spatially compact stencils. The system of linear equations resulting from discretizing the time-harmonic wave equation is solved with a parallel direct solver that computes monochromatic wavefields efficiently for many sources. Dispersion analysis shows that four grid points per P-wavelength provide sufficiently accurate solutions in homogeneous media. The absorbing boundary conditions are perfectly matched layers (PMLs). The kinematic and dynamic accuracy of the method wasassessed with several synthetic examples which illustrate the propagation of S-waves excited at the source or at seismic discontinuities when [Formula: see text]. In frequency-domain modeling with absorbing boundary conditions, the unstable S-wave mode is not excited when [Formula: see text], allowing stable simulations of the P-wave mode for such anisotropic media. Some S-wave instabilities are seen in the PMLs when the symmetry axis is tilted and [Formula: see text]. These instabilities are consistent with previous theoretical analyses of PMLs in anisotropic media; they are removed if the grid interval is matched to the P-wavelength that leads to dispersive S-waves. Comparisons between seismograms computed with the frequency-domain acoustic TTI method and a finite-difference, time-domain method for the vertical transversely isotropic elastic equation show good agreement for weak to moderate anisotropy. This suggests the method can be used as a forward problem for viscoacoustic anisotropic full-waveform inversion.

Factorization-based sparse solvers and preconditioners

Efficient solution of large-scale, ill-conditioned and highly-indefinite algebraic equations often relies on high quality preconditioners together with iterative solvers. Because of their robustness, the factorization-based algorithms could play a significant role when they are combined with iterative methods, particularly in the development of scalable solvers. We present our recent work in using the direct solver SuperLU code base to develop a new supernode-based ILU preconditioner and a domain-decomposition hybrid solver. Our ILU preconditioner is a modification of the classic ILUTP approach, incorporating a number of techniques to improve robustness and performance, which include new dropping strategies that accommodate the use of supernodal structure in the factored matrix. Our hybrid solver is based on the Schur complement method. We use parallel graph partitioning to obtain hierarchical interface/domain decomposition, and multiple parallel direct solvers to solve the subdomain problems simultaneously, and parallel preconditioned iterative solvers to solve the interface problem. We will demonstrate the effectiveness of our new techniques by applying them to two SciDAC applications, modeling next-generation particle accelerators and fusion devices.

A parallel algorithm for direct solution of large sparse linear systems, well suitable to domain decomposition methods

The goal of this paper is to develop a parallel algorithm for the direct solution of large sparse linear systems and integrate it into domain decomposition methods. The computational effort for these linear systems, often encountered in numerical simulation of structural mechanics problems by finite element codes, is very significant in terms of run-time and memory requirements.In this paper, a two-level parallelism is exploited. The exploitation of the lower level of parallelism is based on the development of a parallel direct solver with a nested dissection algorithm and to introduce it into the FETI methods. This direct solver has the advantage of handling zero-energy modes in floating structures automatically and properly. The upper level of parallelism is a coarse-grain parallelism between substructures of FETI. Some numerical tests are carried out to evaluate the performance of the direct solver.

A parallel algorithm for direct solution of large sparse linear systems, well suitable to domain decomposition methods

The goal of this paper is to develop a parallel algorithm for the direct solution of large sparse linear systems and integrate it into domain decomposition methods. The computational effort for these linear systems, often encountered in numerical simulation of structural mechanics problems by finite element codes, is very significant in terms of run-time and memory requirements.In this paper, a two-level parallelism is exploited. The exploitation of the lower level of parallelism is based on the development of a parallel direct solver with a nested dissection algorithm and to introduce it into the FETI methods. This direct solver has the advantage of handling zero-energy modes in floating structures automatically and properly. The upper level of parallelism is a coarse-grain parallelism between substructures of FETI. Some numerical tests are carried out to evaluate the performance of the direct solver.