Omar Ghattas

Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers

This paper reports on the development of a parallel numerical methodology for simulating large-scale earthquake-induced ground motion in highly heterogeneous basins. We target large sedimentary basins with contrasts in wavelengths of over an order of magnitude. Regular grid methods prove intractable for such problems. We overcome the problem of multiple physical scales by using unstructured finite elements on locally-resolved Delaunay triangulations derived from octree-based grids. The extremely large mesh sizes require special mesh generation techniques. Despite the method’s multiresolution capability, large problem sizes necessitate the use of distributed memory parallel supercomputers to solve the elastic wave propagation problem. We have developed a system that helps automate the task of writing efficient portable unstrucmred mesh solvers for distributed memory parallel supercomputers. The numerical methodology and software system have been used to simulate the seismic response of the San Fernando Valley in Southern California to an aftershock of the 1994 Northridge Earthquake. We report on parallel performance on the Cray T3D for several models of the basin ranging in size from 35 000 to 77 million tetrahedra. The results indicate that, despite the highly irregular structure of the problem, excellent performance and scalability are achieved.

p4est : Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees

We present scalable algorithms for parallel adaptive mesh refinement and coarsening (AMR), partitioning, and 2:1 balancing on computational domains composed of multiple connected two-dimensional quadtrees or three-dimensional octrees, referred to as a forest of octrees. By distributing the union of octants from all octrees in parallel, we combine the high scalability proven previously for adaptive single-octree algorithms with the geometric flexibility that can be achieved by arbitrarily connected hexahedral macromeshes, in which each macroelement is the root of an adapted octree. A key concept of our approach is an encoding scheme of the interoctree connectivity that permits arbitrary relative orientations between octrees. Based on this encoding we develop interoctree transformations of octants. These form the basis for high-level parallel octree algorithms, which are designed to interact with an application code such as a numerical solver for partial differential equations. We have implemented and tested these algorithms in the p4est software library. We demonstrate the parallel scalability of p4est on its own and in combination with two geophysics codes. Using p4est we generate and adapt multioctree meshes with up to

Parallel Lagrange--Newton--Krylov--Schur Methods for PDE-Constrained Optimization. Part I: The Krylov--Schur Solver

5.13\times10^{11}

The Dynamics of Plate Tectonics and Mantle Flow: From Local to Global Scales

octants on as many as 220,320 CPU cores and execute the 2:1 balance algorithm in less than 10 seconds per million octants per process.

A Stochastic Newton MCMC Method for Large-Scale Statistical Inverse Problems with Application to Seismic Inversion

Large-scale optimization of systems governed by partial differential equations (PDEs) is a frontier problem in scientific computation. Reduced quasi-Newton sequential quadratic programming (SQP) methods are state-of-the-art approaches for such problems. These methods take full advantage of existing PDE solver technology and parallelize well. However, their algorithmic scalability is questionable; for certain problem classes they can be very slow to converge. In this two-part article we propose a new method for steady-state PDE-constrained optimization, based on the idea of using a full space Newton solver combined with an approximate reduced space quasi-Newton SQP preconditioner. The basic components of the method are Newton solution of the first-order optimality conditions that characterize stationarity of the Lagrangian function; Krylov solution of the Karush--Kuhn--Tucker (KKT) linear systems arising at each Newton iteration using a symmetric quasi-minimum residual method; preconditioning of the KKT system using an approximate state/decision variable decomposition that replaces the forward PDE Jacobians by their own preconditioners, and the decision space Schur complement (the reduced Hessian) by a BFGS approximation initialized by a two-step stationary method. Accordingly, we term the new method {\it Lagrange--Newton--Krylov--Schur} (LNKS). It is fully parallelizable, exploits the structure of available parallel algorithms for the PDE forward problem, and is locally quadratically convergent. In part I of this two-part article, we investigate the effectiveness of the KKT linear system solver. We test our method on two optimal control problems in which the state constraints are described by the steady-state Stokes equations. The objective is to minimize dissipation or the deviation from a given velocity field; the control variables are the boundary velocities. Numerical experiments on up to 256 Cray T3E processors and on an SGI Origin 2000 include scalability and performance assessment of the LNKS algorithm and comparisons with reduced SQP for up to

High Resolution Forward And Inverse Earthquake Modeling on Terascale Computers

1,000,000

From mesh generation to scientific visualization: an end-to-end approach to parallel supercomputing

state and 50,000 decision variables. In part II of the article, we address globalization and inexactness issues, and apply LNKS to the optimal control of the steady incompressible Navier--Stokes equations.

Parameter and State Model Reduction for Large-Scale Statistical Inverse Problems

Improving Earth Models The geophysical processes responsible for shaping the planets surface and interior need largescale simulations, but to achieve high resolution at these scales is costly and tends to focus on gradual processes such as plate tectonics. By using large parallel supercomputers, Stadler et al. (p. 1033; see the Perspective by Becker; see the cover) have improved on a commonly used method—adaptive mesh refinement—to increase the resolution of global geodynamic models to the scale of a single kilometer and been able to reveal unexpected insights into localized processes, such as subduction zone mechanics, thermal anomalies in the lower mantle, and the speed of movement of oceanic plates. Computational advances enable the modeling of global geophysical processes to the scale of a kilometer. Plate tectonics is regulated by driving and resisting forces concentrated at plate boundaries, but observationally constrained high-resolution models of global mantle flow remain a computational challenge. We capitalized on advances in adaptive mesh refinement algorithms on parallel computers to simulate global mantle flow by incorporating plate motions, with individual plate margins resolved down to a scale of 1 kilometer. Back-arc extension and slab rollback are emergent consequences of slab descent in the upper mantle. Cold thermal anomalies within the lower mantle couple into oceanic plates through narrow high-viscosity slabs, altering the velocity of oceanic plates. Viscous dissipation within the bending lithosphere at trenches amounts to ~5 to 20% of the total dissipation through the entire lithosphere and mantle.

A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media

We address the solution of large-scale statistical inverse problems in the framework of Bayesian inference. The Markov chain Monte Carlo (MCMC) method is the most popular approach for sampling the posterior probability distribution that describes the solution of the statistical inverse problem. MCMC methods face two central difficulties when applied to large-scale inverse problems: first, the forward models (typically in the form of partial differential equations) that map uncertain parameters to observable quantities make the evaluation of the probability density at any point in parameter space very expensive; and second, the high-dimensional parameter spaces that arise upon discretization of infinite-dimensional parameter fields make the exploration of the probability density function prohibitive. The challenge for MCMC methods is to construct proposal functions that simultaneously provide a good approximation of the target density while being inexpensive to manipulate. Here we present a so-called Stoch...

Goal-oriented, model-constrained optimization for reduction of large-scale systems

For earthquake simulations to play an important role in the reduction of seismic risk, they must be capable of high resolution and high fidelity. We have developed algorithms and tools for earthquake simulation based on multiresolution hexahedral meshes. We have used this capability to carry out 1 Hz simulations of the 1994 Northridge earthquake in the LA Basin using 100 million grid points. Our wave propagation solver sustains 1.21 teraflop/s for 4 hours on 3000 AlphaServer processors at 80% parallel efficiency. Because of uncertainties in characterizing earthquake source and basin material properties, a critical remaining challenge is to invert for source and material parameter fields for complex 3D basins from records of past earthquakes. Towards this end, we present results for material and source inversion of high-resolution models of basins undergoing antiplane motion using parallel scalable inversion algorithms that overcome many of the difficulties particular to inverse heterogeneous wave propagation problems.