With increasingly finescale parameterization of global tomographic models comes growth of the linear systems that underlie most iterative inversion schemes (e.g. SEMum: Lekic and Romanowicz, 2010). As assembly and solution of these systems becomes intractable on modern workstations, distributedmemory parallel solvers provide an attractive alternative. We summarize the development of a parallel solver tailored to the regularized GaussNewton scheme often used in global tomography.
In a separate effort, we explore alternative datapartitioning schemes for waveform inversion using a synthetic dataset. The SPICE tomographic benchmark was developed to evaluate the effects of inversion methodology on the resolving power of tomographic images (Qin et al., 2008). We discuss our use of the simpler SPICE ``preliminary'' waveform dataset to examine datapartitioning schemes different from those employed in Berkeley global tomographic efforts to date.
In waveform tomography, we seek to invert a dataset of longperiod waveforms
for a model of earth structure . The nonlinear dependence of the
wavefield upon requires structural sensitivity to be
linearized and iteratively estimated.
The generalized leastsquares formalism is popular choice for such
problems (e.g. Tarantola, 2005), and gives rise to the following
regularized GaussNewton scheme:
While possesses a sizable zero component, it is generally nonsparse. Thus, the by matrix , where is the dimension of , is likely fully dense. In the case of the SEMum model of Lekic and Romanowicz (2010), where , wallclock times of approximately one day were observed during assembly and sequential solution of the above system using MATLAB on a modern workstation. Further, repeated solution is required at each linearized iteration ( ) in order to evaluate model sensitivity to regularization. Thus, parallel solution for the regularized GaussNewton update described above is a desirable capability.
A parallel solver mimicking the functionality of the sequential code was
developed and validated, utilizing the Scalable Linear Algebra PACKage
(ScaLAPACK
)  named for its attractive weakscaling behavior on large dense problems
such as these (Choi et al., 1996). A robust parallel implementation of LU
decomposition with partial pivoting was used. Benchmark tests were
performed for realizations of the above system arising in the development of
SEMum, and wallclock solve times for a range of processor grids are shown in
Table 2.1.

Perhaps most importantly, this savings in computetime allows for faster and more detailed enumeration of model sensitivity to regularization  hopefully leading to greater flexibility in scheduling of compute resources for far heavier workloads, such as wavefield modeling with the Spectral Element Method (SEM: e.g. Komatitsch and Vilotte, 1998).
The simple isotropic model underlying the SPICE ``preliminary'' benchmark contains three large ( km) and five small ( km) columnar perturbations of % and % respectively, from which and are scaled (Qin et al., 2008). These structures have roughly Gaussian profile and extend from the surface to the CMB. The associated SEM synthetic dataset is valid to 50 s period, which we further bandpass filtered with corners at 80 and 250 s and cutoff at 60 and 400 s. Following the waveform partitioning and weighting scheme of Li and Romanowicz (1996), first and second orbit Rayleigh wave fundamental mode and overtone wavepackets were selected from vertical component waveforms for all 27 synthetic events recorded at the 256 stations in the dataset.
While maintaining a 1D crust, we inverted for mantle structure using the pathaverage approximation (PAVA: Woodhouse and Dziewonski, 1984) and nonlinear asymptotic coupling theory (NACT: Li and Romanowicz, 1995) for fundamental mode and overtone data, respectively. We partitioned the inversion at km, below which the resolving power of fundamental mode data rapidly decreases, and adjusted the relative weight assigned to the fundamental mode and overtone datasets in each partition. We found that we could more effectively extend model sensitivity into the transition zone over simply merging the datasets with a uniform weighting. The cubic spline model parametrization was kept continuous across the partition so that the model would remain smooth. This represents a methodology broadly consistent with Berkeley tomographic studies to date.
To explore alternative approaches to optimizing depth sensitivity, we proposed a simple multifrequency scheme. Overlapping filters with cutoff periods of 60/140 s, 120/200 s, and 180/260 s were applied to SPICE vertical component waveforms, which were further weighted according to path redundancy as in Li and Romanowicz (1996), as well as maximum waveform amplitude. As the filtered waveforms are dominated by fundamental mode surface waves, we chose to invert full waveforms using PAVA only. During the inversion, sensitivities of each of the filtered datasets were independently weighted to maximize smoothness of total sensitivity with depth. Unsurprisingly, resolution of structure through the transition zone, best constrained by overtone surface waves, was notably degraded relative to the wavepacketpartitioned inversion. A far more complex multifrequency scheme (e.g. AMI: Lebedev et al., 2005) would likely be required in order to compete with the wavepacketbased scheme.
As approximate forward modeling schemes are supplanted by SEM in waveform tomography (e.g. SEMum: Lekic and Romanowicz, 2010), mitigation of computational cost takes on new importance. One proposed solution is the sourcestacking technique (Capdeville et al., 2005). Linearity of the elastic wave equation with respect to source implies that SEM waveforms stacked over individual sources are equivalent to those resulting from simultaneous triggering of sources in a single simulation. Thus, inversion of sourcestacked traces may provide a novel avenue for waveform tomography at vastly reduced cost. As a naïve test, we sought to invert the SPICE preliminary benchmark using PAVA for both the sensitivity kernels and as a standin for SEM.
Vertical component SPICE synthetic waveforms were bandpassed with corners at 120 and 250 s (cutoff at 100 and 400 s), aligned on event origin time, weighted by scalar moment magnitude, and stacked over all 27 events. PAVA sensitivity kernels and synthetics were similarly weighted and stacked. Using these approximate methods alone, the structure of the SPICE model became recognizable as early as the second inversion iteration (see Figure 2.30). Of course, while these results are highly preliminary and do not reflect realistic data reliability and noise conditions, they reiterate the promise of the source stacking method. Overcoming the shortcomings present in realistic datasets, as well as developing methods for optimized eventstation clustering, are of primary concern moving forward.
This research was supported by National Science Foundation grant NSF/EAR0738284.
Capdeville, Y., Y. Gung, and B. Romanowicz, Towards global earth tomography using the spectral element method: a technique based on source stacking, Geophys. J. Int, 162, 541554, 2005.
Choi, J., J. Demmel, I. Dhillon, J. Dongarra, S. Ostrouchov, A. Petitet, K. Stanley, D. Walker, and R.C. Whaley, ScaLAPACK: a portable linear algebra library for distributed memory computers  design issues and performance, Comp. Phys. Comm., 97, 115, 1996.
Komatitsch, D. and J. Vilotte, The spectral element method: an efficient tool to simulate the seismic responds of 2D and 3D geological structures, Bull. Seism. Soc. Am., 88, 368392, 1998.
Lebedev, S., G. Nolet, T. Meier, and R.D. van der Hilst, Automated multimode inversion of surface and S waveforms, Geophys. J. Int., 162, 951964, 2005.
Lekic, V. and B. Romanowicz, Inferring upper mantle structure by full waveform tomography with the Spectral Element Method, Submitted to Geophys. J. Int., 2010.
Li, X.D. and B. Romanowicz, Comparison of global waveform inversions with and without considering crossbranch modal coupling, Geophys. J. Int., 121, 695709, 1995.
Li, X.D. and B. Romanowicz, Global mantle shearvelocity model developed using nonlinear asymptotic coupling theory, Geophys. J. R. Astr. Soc., 101, 22,24522,272, 1996.
Qin, Y., Y. Capdeville, V. Maupin, J.P. Montagner, S. Lebedev, and E. Beucler, SPICE Benchmark for global tomographic methods, Geophys. J. Int., 175, 598616, 2008.
Tarantola, A., Inverse problem theory and methods for model parameter estimation, SIAM, 2005.
Woodhouse, J.H. and A.M. Dziewonski, Mapping the upper mantle: Three dimensional modeling of earth structure by inversion of seismic waveforms, J. Geophys. Res, 89, 59535986, 1984.
Berkeley Seismological Laboratory
215 McCone Hall, UC Berkeley, Berkeley, CA 947204760
Questions or comments? Send email: www@seismo.berkeley.edu
© 2007, The Regents of the University of California