Jacobo Bielak

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.

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory

This article reports on the development of a modular two-step, finite- element methodology for modeling earthquake ground motion in highly heteroge- neous localized regions with large contrasts in wavelengths. We target complex geo- logical structures such as sedimentary basins and ridges that are some distance away from the earthquake source. We overcome the problem of multiple physical scales by subdividing the original problem into two simpler ones. The first is an auxiliary problem that simulates the earthquake source and propagation path effects with a model that encompasses the source and a background structure from which the lo- calized feature has been removed. The second problem models local site effects. Its input is a set of equivalent localized forces derived from the first step. These forces act only within a single layer of elements adjacent to the interface between the exterior region and the geological feature of interest. This enables us to reduce the domain size in the second step. If the background subsurface structure is simple, one can replace the finite-element method in the first step with an alternative efficient method. The methodology is illustrated in a companion paper (Yoshimura et al., 2003) for several 3D problems of increasing physical and computational complexity. We consider first a flat-layered, stratigraphic system. For this simple case, the first step can be carried out by means of 3D Greens function evaluations. The extension to more general problems is illustrated by two examples: a basin and a hill, with the same background stratigraphy. To verify the two-step procedure with a problem for which the finite-element method is used throughout, we model ground motion in a small region of the Los Angeles Basin, using both the two-step domain-reduction method and the traditional approach in which the computational domain contains both the source and the geological region of interest.

High Resolution Forward And Inverse Earthquake Modeling on Terascale Computers

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.

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

Parallel supercomputing has traditionally focused on the inner kernel of scientific simulations: the solver. The front and back ends of the simulation pipeline - problem description and interpretation of the output - have taken a back seat to the solver when it comes to attention paid to scalability and performance, and are often relegated to offline, sequential computation. As the largest simulations move beyond the realm of the terascale and into the petascale, this decomposition in tasks and platforms becomes increasingly untenable. We propose an end-to-end approach in which all simulation components - meshing, partitioning, solver, and visualization - are tightly coupled and execute in parallel with shared data structures and no intermediate I/O. We present our implementation of this new approach in the context of octree-based finite element simulation of earthquake ground motion. Performance evaluation on up to 2048 processors demonstrates the ability of the end-to-end approach to overcome the scalability bottlenecks of the traditional approach

A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion

We present a nonlinear optimization method for large-scale 3D elastic full-waveform seismic inversion. The method combines outer Gauss–Newton nonlinear iterations with inner conjugate gradient linear iterations, globalized by an Armijo backtracking line search, solved on a sequence of finer grids and higher frequencies to remain in the vicinity of the global optimum, inexactly terminated to prevent oversolving, preconditioned by L-BFGS/Frankel, regularized by a total variation operator to capture sharp interfaces, finely discretized by finite elements in the Lame parameter space to provide flexibility and avoid bias, implemented in matrix-free fashion with adjoint-based computation of reduced gradient and reduced Hessian-vector products, checkpointed to avoid full spacetime waveform storage, and partitioned spatially across processors to parallelize the solutions of the forward and adjoint wave equations and the evaluation of gradient-like information. Several numerical examples demonstrate the grid independence of linear and nonlinear iterations, the effectiveness of the preconditioner, the ability to solve inverse problems with up to 17 million inversion parameters on up to 2048 processors, the effectiveness of multiscale continuation in keeping iterates in the basin of attraction of the global minimum, and the ability to fit the observational data while reconstructing the model with reasonable resolution and capturing sharp interfaces.

Model for Basin Effects on Long-Period Response Spectra in Southern California

We propose a model for the effect of sedimentary basin depth on long-period response spectra. The model is based on the analysis of 3-D numerical simulations (finite element and finite difference) of long-period (2–10 s) ground motions for a suite of sixty scenario earthquakes (Mw 6.3 to Mw 7.1) within the Los Angeles basin region. We find depth to the 1.5 km/s S-wave velocity isosurface to be a suitable predictor variable, and also present alternative versions of the model based on depths to the 1.0 and 2.5 km/s isosurfaces. The resulting mean basin-depth effect is period dependent, and both smoother (as a function of period and depth) and higher in amplitude than predictions from local 1-D models. The main requirement for the use of the results in construction of attenuation relationships is determining the extent to which the basin effect, as defined and quantified in this study, is already accounted for implicitly in existing attenuation relationships, through (1) departures of the average “rock” site from our idealized reference model, and (2) correlation of basin depth with other predictor variables (such as Vs30).

Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part II: Verification and Applications

Several examples are used to verify the domain reduction method (DRM), a two-step finite-element methodology described in a companion article for modeling earthquake ground motion in highly heterogeneous three-dimensional localized regions. The first set involves a simple, flat-layered system. Verification of the DRM for this problem is carried out by comparing the results to those calculated directly by the theoretical Green9s function method. Its applicability to more general problems is illustrated by two examples: a basin and a hill with and without a weathered surface layer and with the same stratigraphy. We use a Green9s function approach for the first step, which for the examples under consideration needs to be performed only once. For the second step, the domain of computation is reduced in each case to a small neighborhood of the geological feature at hand. The second application considers the ground motion due to a strike-slip double couple buried 14 km below the free surface in an 80-km × 80-km × 30-km region that encloses entirely the Los Angeles basin. This problem is solved first by the finite-element method using the single-step traditional approach, in which the ground motion is calculated simultaneously near the seismic source, along the propagation path, and within the region of interest with a single model that encompasses the entire geological structure, from the source to the region of interest. The DRM is then used to determine anew the ground motion over a much smaller (6-km × 6-km × 0.6-km) region contained within the original domain, and the results of the two methods within this region of interest are compared. These examples serve to demonstrate that in many applications the DRM can be significantly more efficient than the traditional approach. The DRM can be particularly advantageous (1) if the source is far from the local structure and the local structure is much softer than that of the exterior region, (2) if the localized feature exhibits nonlinear behavior, or (3) if for a prescribed source, one wishes to consider a sequence of simulations in which the properties of the local feature, which might include man-made structures, are varied from one simulation to the next.

Earthquake ground motion modeling on parallel computers

We describe the design and discuss the performance of a parallel elastic wave propagation simulator that is being used to model and study earthquake-induced ground motion in large sedimentary basins. The components of the system include mesh generators, a mesh partitioner, a parceler, and a parallel code generator, as well as parallel numerical methods for applying seismic forces, incorporating absorbing boundaries, and solving the discretized wave propagation problem. We discuss performance on the Cray T3D for unstructured mesh wave propagation problems of up to 14 million tetrahedra. By paying careful attention to each step of the process, we obtain excellent performance despite the highly irregular structure of the coefficient matrices of the problem. The mesh generator, partitioner, parceler, and code generator have been incorporated into an integrated toolset/compiler. This system, called Archimedes, automates the solution of unstructured mesh PDE problems on parallel computers, and is being used for other unstructured mesh applications beyond ground motion modeling.

Two-Dimensional Nonlinear Earthquake Response Analysis of a Bridge-Foundation-Ground System

This paper presents a two-dimensional advanced nonlinear FE model of an actual bridge, the Humboldt Bay Middle Channel (HBMC) Bridge, and its response to seismic input motions. This computational model is developed in the new structural analysis software framework OpenSees. The foundation soil is included to incorporate soil-foundation-structure interaction effects. Realistic nonlinear constitutive models for cyclic loading are used for the structural (concrete and reinforcing steel) and soil materials. The materials in the various soil layers are modeled using multi-yield-surface plasticity models incorporating liquefaction effects. Lysmer-type absorbing/transmitting boundaries are employed to avoid spurious wave reflections along the boundaries of the computational soil domain. Both procedures and results of earthquake response analysis are presented. The simulation results indicate that the earthquake response of the bridge is significantly affected by inelastic deformations of the supporting soil medium due to lateral spreading induced by soil liquefaction.

Full Waveform Inversion for Seismic Velocity and Anelastic Losses in Heterogeneous Structures

We present a least-squares optimization method for solving the nonlinear full waveform inverse problem of determining the crustal velocity and intrinsic at- tenuation properties of sedimentary valleys in earthquake-prone regions. Given a known earthquake source and a set of seismograms generated by the source, the in- verse problem is to reconstruct the anelastic properties of a heterogeneous medium with possibly discontinuous wave velocities. The inverse problem is formulated as a constrained optimization problem, where the constraints are the partial and ordinary differential equations governing the anelastic wave propagation from the source to the receivers in the time domain. This leads to a variational formulation in terms of the material model plus the state variables and their adjoints. We employ a wave propaga- tion model in which the intrinsic energy-dissipating nature of the soil medium is mod- eled by a set of standard linear solids. The least-squares optimization approach to inverse wave propagation presents the well-known difficulties of ill posedness and multiple minima. To overcome ill posedness, we include a total variation regulariza- tion functional in the objective function, which annihilates highly oscillatory material property components while preserving discontinuities in the medium. To treat multi- ple minima, we use a multilevel algorithm that solves a sequence of subproblems on increasingly finer grids with increasingly higher frequency source components to re- main within the basin of attraction of the global minimum. We illustrate the metho- dology with high-resolution inversions for two-dimensional sedimentary models of the San Fernando Valley, under SH-wave excitation. We perform inversions for both the seismic velocity and the intrinsic attenuation using synthetic waveforms at the observer locations as pseudoobserved data.