HAL Id: hal-00617249
https://hal.archives-ouvertes.fr/hal-00617249
Submitted on 18 Jun 2021
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Forward and adjoint simulations of seismic wave
propagation on fully unstructured hexahedral meshes
Daniel Peter, Dimitri Komatitsch, Yang Luo, Roland Martin, Nicolas Le Go,
Emanuele Casarotti, Pieyre Le Loher, Federica Magnoni, Qinya Liu, Céline
Blitz, et al.
To cite this version:
Daniel Peter, Dimitri Komatitsch, Yang Luo, Roland Martin, Nicolas Le Go, et al.. Forward and
adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes. Geophys-
ical Journal International, Oxford University Press (OUP), 2011, 186, pp.721-739. �10.1111/j.1365-
246X.2011.05044.x�. �hal-00617249�
Geophys. J. Int. (2011) 186, 721–739 doi: 10.1111/j.1365-246X.2011.05044.x
GJI Seismology
Forward and adjoint simulations of seismic wave propagation
on fully unstructured hexahedral meshes
Daniel Peter,
1
Dimitri Komatitsch,
2,3
Yang Luo,
1
Roland Martin,
2
Nicolas Le Goff,
2
Emanuele Casarotti,
4
Pieyre Le Loher,
2
Federica Magnoni,
4
Qinya Liu,
5
C
´
eline Blitz,
2
Tarje Nissen-Meyer,
6
Piero Basini
6
and Jeroen Tromp
1,7
1
Princeton University, Department of Geosciences, 318 Guyot Hall, Princeton, NJ 08544, USA. E-mail: dpeter@princeton.edu
2
Universit
´
e de Pau et des Pays de l’Adour, CNRS & INRIA Magique-3D, Laboratoire de Mod
´
elisation et d’Imagerie en G
´
eosciences UMR 5212, Avenue de
l’Universit
´
e, 64013 Pau Cedex, France
3
Institut universitaire de France, 103 boulevard Saint-Michel, 75005 Paris, France
4
Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143, Rome, Italy
5
Department of Physics, University of Toronto, Ontario, Canada
6
Institute of Geophysics, ETH Zurich, Sonneggstr. 5, CH-8092 Zurich, Switzerland
7
Princeton University, Program in Applied & Computational Mathematics, Princeton, NJ 08544, USA
Accepted 2011 April 12. Received 2011 April 8; in original form 2011 February 2
SUMMARY
We present forward and adjoint spectral-element simulations of coupled acoustic and
(an)elastic seismic wave propagation on fully unstructured hexahedral meshes. Simulations
benefit from recent advances in hexahedral meshing, load balancing and software optimiza-
tion. Meshing may be accomplished using a mesh generation tool kit such as CUBIT, and
load balancing is facilitated by graph partitioning based on the SCOTCH library. Coupling
between fluid and solid regions is incorporated in a straightforward fashion using domain
decomposition. Topography, bathymetry and Moho undulations may be readily included in
the mesh, and physical dispersion and attenuation associated with anelasticity are accounted
for using a series of standard linear solids. Finite-frequency Fr
´
echet derivatives are calculated
using adjoint methods in both fluid and solid domains. The software is benchmarked for a
layercake model. We present various examples of fully unstructured meshes, snapshots of
wavefields and finite-frequency kernels generated by Version 2.0 ‘Sesame’ of our widely used
open source spectral-element package SPECFEM3D.
Key words: Tomography; Interferometry; Computational seismology; Wave propagation.
1 INTRODUCTION
We present a new software package, SPECFEM3D Version 2.0
‘Sesame’, capable of simulating forward and adjoint seismic wave
propagation on fully unstructured hexahedral meshes of arbitrary
shaped model domains. In view of unrelenting growth in compu-
tational power, it has become more and more important to develop
software capable of harnessing powerful computers to address a
broad range of seismological forward and inverse problems. A well-
established numerical technique for solving such problems in a fast
and highly accurate manner is the spectral-element method (SEM).
The SEM was originally developed in computational fluid dynam-
ics (Patera 1984; Maday & Patera 1989) and has been successfully
adapted to address problems in seismic wave propagation. Early
seismic wave propagation applications of the SEM, utilizing Leg-
endre basis functions and a perfectly diagonal mass matrix, include
Cohen et al. (1993), Komatitsch (1997), Faccioli et al. (1997),
Casadei & Gabellini (1997), Komatitsch & Vilotte (1998) and
Komatitsch & Tromp (1999), whereas applications involving
Chebyshev basis functions and a non-diagonal mass matrix include
Seriani & Priolo (1994), Priolo et al. (1994) and Seriani et al. (1995).
The SEM is a continuous Galerkin technique, which may be
made discontinuous (Bernardi et al. 1994; Chaljub 2000; Kopriva
et al. 2002; Chaljub et al. 2003; Legay et al. 2005; Kopriva 2006;
Wilcox et al. 2010; Acosta Minolia & Kopriva 2011); it is then close
to a particular case of the discontinuous Galerkin technique (Reed
& Hill 1973; Ar nold 1982; Falk & Richter 1999; Hu et al. 1999;
Cockburn et al. 2000; Giraldo et al. 2002; Rivi
´
ere & Wheeler 2003;
Monk & Richter 2005; Grote et al. 2006; Ainsworth et al. 2006;
Bernacki et al. 2006; Dumbser & K
¨
aser 2006; De Basabe et al.
2008; de la Puente et al. 2009; Wilcox et al. 2010; De Basabe &
Sen 2010; Etienne et al. 2010), with optimized efficiency because
of its tensorized basis functions (Wilcox et al. 2010; Acosta Minolia
& Kopriva 2011).
An important feature of the SEM is that it can accurately han-
dle very distorted mesh elements (Oliveira & Seriani 2011), and
C
2011 The Authors 721
Geophysical Journal International
C
2011 RAS
Geophysical Journal International
Downloaded from https://academic.oup.com/gji/article/186/2/721/590417 by guest on 18 June 2021
722 D. Peter et al.
thus conforming non-structured mesh doubling bricks can effi-
ciently accommodate mesh size variations (Komatitsch & Tromp
2002a, 2004; Lee et al. 2008, 2009a,b). The method has very good
accuracy and convergence properties, such as a spectral rate of
convergence (Canuto et al. 1988; Maday & Patera 1989; Seriani
& Priolo 1994; Deville et al. 2002; Cohen 2002; De Basabe &
Sen 2007; Seriani & Oliveira 2008). In this sense the SEM is
close to the family of pseudo-spectral methods (see e.g. Canuto
et al. 1988; Carcione et al. 1988a, 1992; Carcione & Wang 1993;
Komatitsch et al. 1996), but combined with the flexibility of fi-
nite elements, in particular in terms of mesh design. For reviews
of the SEM in seismology, see for example, Komatitsch et al.
(2005), Chaljub et al. (2007), Tromp et al. (2008) and Fichtner
(2010).
The SEM is well suited to parallel implementations on very large
supercomputers (Komatitsch & Tromp 2002a; Komatitsch et al.
2003; Tsuboi et al. 2003; Komatitsch et al. 2008; Carrington et al.
2008; Komatitsch et al. 2010b) as well as on clusters of GPU accel-
erating g raphics cards (Komatitsch et al. 2009, 2010a, Komatitsch
2011). Tensor products inside each element may be optimized to
reach very high efficiency (Deville et al. 2002), and mesh point and
element numbering may be optimized to reduce processor cache
misses and improve cache reuse (Komatitsch et al. 2008). The
SEM can handle triangular (in 2-D) or tetrahedral (in 3-D) ele-
ments (Wingate & Boyd 1996; Taylor & Wingate 2000; Komatitsch
et al. 2001; Cohen 2002; Mercerat et al. 2006), as well as mixed
meshes, although with increased cost and reduced accuracy in these
non-tensorized elements, as in the discontinuous Galerkin method.
In many cases of practical seismological interest, using a con-
forming mesh and a continuous formulation are sufficient, because
in most geological models material property contrasts are not too
dramatic. When this ceases to be true, requiring a discontinuous for-
mulation, one can either turn to a discontinuous version of the SEM
(Bernardi et al. 1994; Chaljub 2000; Kopriva et al. 2002; Chaljub
et al. 2003; Legay et al. 2005; Kopriva 2006; Wilcox et al. 2010;
Acosta Minolia & Kopriva 2011) or to a discontinuous Galerkin
technique. A discontinuous formulation is particularly suitable for
dynamic rupture simulations, because high frequencies or super-
shear rupture need to be accommodated near the fault, where a
significantly denser mesh and a more sophisticated (upwind) time
scheme are required, thereby suppressing the amplification of un-
stable modes (see e.g. Benjemaa et al. 2007, 2009; de la Puente
et al. 2009; Tago et al. 2010). Another example that may require
a discontinuous formulation involves the resolution of a shallow
geotechnical layer, in which seismic shear wave speeds may be
reduced by an order of magnitude.
For seismological applications, the SEM has been success-
fully implemented for 3-D global- and regional-scale simulations
(Komatitsch & Vilotte 1998; Paolucci et al. 1999; Chaljub 2000;
Komatitsch & Tromp 2002a,b; Capdeville et al. 2003; Chaljub &
Valette 2004; Fichtner et al. 2009a), as well as local-scale simu-
lations in complex and/or densely populated regions, for example
in southern California, USA (Komatitsch et al. 2004; Tape et al.
2009, 2010), Taipei, Taiwan (Lee et
al. 2008, 2009a,b), Caracas,
Venezuela (Delavaud et al. 2006) and Grenoble, France (Chaljub
et al. 2005; Stupazzini et al. 2009; Chaljub et al. 2010). The SEM
may also be used to study elastic wave propagation on smaller
scales, for instance the propagation of ultrasonic waves in crystals
(van Wijk et al. 2004).
Two complementary SEM software packages—namely,
SPECFEM3D_GLOBE for global and regional simulations,
and SPECFEM3D for local simulations—are feature-rich, well
benchmarked and documented implementations. Data parallelism
in the SEM is efficiently exploited using the Message-Passing
Interface (MPI) standard, crucial for modern high-performance
computing. These open source packages are freely available via the
Computational Infrastructure for Geodynamics (CIG) and widely
used by the seismological community.
To extend the range of local-scale applications, easing the task of
mesh generation is paramount. The two community software pack-
ages separate a simulation into two distinct steps: first, creation
of a hexahedral mesh, and second, solution of the seismic wave
equation. This separation avoids the overhead of remeshing when
running multiple simulations for the same region, for example, re-
peated simulations at the same resolution. Focussing on local-scale
simulations, previous versions of SPECFEM3D used an internal
mesher which was explicitly tied to the specific purposes of the
package: all geological models were based on a layercake model.
Consequently, the solver was restricted by its internal mesher. It was
impossible to run spectral-element simulations on more complex 3-
D models without significant recoding, nor was it possible to run
such simulations in regions of interest for on- and off-shore explo-
ration seismology, because acoustic wave propagation in fluids was
not supported by the package.
The purpose of this paper is to present forward and ad-
joint simulations in various 3-D models using the new soft-
ware package, SPECFEM3D Version 2.0 ‘Sesame’, thereby illus-
trating its current capabilities. The original SPECFEM3D pack-
age for local simulations was extended, improved and opti-
mized in various ways. The Version 2.0 ‘Sesame’ release in-
cludes a more flexible internal mesher and accommodates more
powerful external meshers, such as CUBIT (Blacker et al.
1994; White et al. 1995; Mitchell 1996; Casarotti et al. 2008).
Adding such external meshers into the workflow greatly in-
creases flexibility for high-performance applications, as illustrated
by the GeoELSE software package (Casadei & Gabellini 1997;
Stupazzini et al. 2009; Chaljub et al. 2010). Advantages of
GeoELSE include the accommodation of viscoplastic and non-
linear rheologies, whereas benefits of SPECFEM3D include cou-
pled fluid-solid domains and adjoint capabilities; the latter enable
one to address seismological inverse problems. Load balancing par-
allel simulations in SPECFEM3D is accomplished based on the
graph partitioning software package SCOTCH (Pellegrini & Ro-
man 1996; Chevalier & Pellegrini 2008). The new package facili-
tates coupled forward and adjoint acoustic/(an)elastic simulations,
which are especially interesting for problems in exploration seismol-
ogy, ocean acoustics and medical tomography. The new software is
freely available under the GNU GPL Version 2 license via CIG.
2 GOVERNING EQUATIONS
Let us briefly summarize the equations governing seismic wave
propagation implemented in SPECFEM3D. For more technical
details, the reader is refer red to Komatitsch & Tromp (1999).
SPECFEM3D Version 2.0 ‘Sesame’ implements wave propagation
in coupled (an)elastic and acoustic materials on local scales. We
may thus safely neglect additional effects that would arise from
self-gravitation and rotation (Komatitsch & Tromp 2002b, 2005;
Chaljub et al. 2007), which are important at longer periods. In the
following, we first discuss (an)elastic wave propagation and subse-
quently consider acoustic waves.
C
2011 The Authors, GJI, 186, 721–739
Geophysical Journal International
C
2011 RAS
Downloaded from https://academic.oup.com/gji/article/186/2/721/590417 by guest on 18 June 2021
SPECFEM3D Version 2.0 ‘Sesame’ 723
2.1 Elastic domain
For elastic materials, the displacement wavefield s(x, t)isgoverned
by
ρ∂
2
t
s =∇·T + f, (1)
where ρ denotes mass density, T the stress tensor and f the seismic
source. On free surfaces, the traction vector must vanish, that is,
ˆ
n · T = 0, (2)
where
ˆ
n denotes the unit outward normal on the surface. On bound-
aries between different elastic materials, both traction
ˆ
n ·T and dis-
placement s need to be continuous. On boundaries between elastic
and acoustic domains, traction
ˆ
n · T and the normal component of
displacement
ˆ
n ·s need to be continuous. The initial conditions are
s(x, 0) = 0,∂
t
s(x, 0) = 0. (3)
We thus initiate the simulation in a medium at rest. To accommodate
simulations under pre-stressed conditions, these initial conditions
may be modified in an appropriate manner.
For elastic materials, the force f in eq. (1) represents the earth-
quake, which for a simple point source may be written as
f =−M ·∇δ(x − x
s
) S(t), (4)
where M denotes the moment tensor, x
s
the source location, δ(x −x
s
)
the Dirac delta distribution located at x
s
and S(t) the source-time
function. The software also accommodates kinematic rupture sim-
ulations, which may be captured by prescribing a moment-density
tensor field.
The stress tensor T is linearly related to the strain via the consti-
tutive relationship
T = c : ∇s , (5)
where c denotes the stiffness tensor that describes the elastic prop-
erties of the medium. The implementation is general and can handle
a fully anisotropic tensor with 21 independent parameters (Chen &
Tromp 2007; Sieminski et al. 2007a,b). Using a linear constitutive
relationship is valid under the assumption that perturbations to the
reference state are small. Note that non-linear effects are sometimes
observed, for example, non-linear soil amplification, and non-linear
constitutive relationships become important for studying such ef-
fects, for example, for risk mitigation (Xu et al. 2003; Dupros et al.
2010).
In an anelastic medium, we approximate an absorption-band solid
using a series of L standard linear solids (Liu et al. 1976), and model
the time evolution of the isotropic shear modulus μ by
μ(t) = μ
R
1 −
L
l=1
1 −
τ
l
τ
σ
l
e
−t/τ
σ
l
H(t) , (6)
where μ
R
denotes the relaxed modulus, H(t) the Heaviside function
and τ
σ
l
& τ
l
the stress and strain relaxation times of the lth standard
linear solid. Experience shows that three solids generally suffice
for simulating an absorption band (Emmerich & Korn 1987). For
further details, see Carcione et al. (1988b), Robertsson (1996),
Day & Bradley (2001), Moczo & Kristek (2005), Komatitsch et al.
(2005), Carcione (2007) and Savage et al. (2010). Simulations of
seismic wave propagation in laboratory-scale rock samples or in the
context of medical tomography involve very high frequencies (in
the kHz or even MHz range), and strong attenuation must be taken
into account.
The SEM solves the equations of motion in the weak form, which
is obtained by dotting the momentum eq. (1) with an arbitrary test
vector w and integrating by parts over the model volume . We focus
on elastic domains and consider coupling interfaces with acoustic
domains. Thus, we obtain
ρ w · ∂
2
t
s d
3
x =
∂
ˆ
n · T · w d
2
x −
∇w : T d
3
x
+ M : ∇w(x
s
) S(t) .
(7)
Note that in this formulation the traction-free surface condition is
implicitly accounted for by setting the contribution from the free
surface to zero.
When and where necessary, we use Clayton–Engquist–Stacey
absorbing conditions (Clayton & Engquist 1977; Stacey 1988;
Quarteroni et al. 1998) to absorb outgoing waves on fictitious
boundaries of the mesh, thereby representing a semi-infinite do-
main. It would be more efficient to use a Perfectly Matched Layer
(PML) (see e.g. Komatitsch & Martin 2007, 2008c; Martin &
Komatitsch 2009), but a parallel implementation with good load-
balancing properties is challenging because additional equations
need to be solved. This issue becomes impor tant when high-order
time marching is required to reduce numerical dispersion in dif-
ficult case studies that involve complex media with poroelastic or
viscoelastic rheologies (Martin et al. 2008b, 2010) or Newtonian
compressible fluids (Martin & Couder-Castaneda 2010). Conse-
quently, additional computations need to be perfor med in PML
layers, in particular in corners, where contributions along several
directions are summed (Komatitsch & Martin 2007).
At a solid–fluid boundary, the interface integral over the coupling
surface ∂ is used to exchange pressure from the fluid p
fluid
to the
solid:
ˆ
n · T =−p
fluid
ˆ
n.
2.2 Acoustic domain
We define a scalar potential φ such that the displacement s may be
written as
s = ρ
−1
∇φ. (8)
The equation of motion in terms of the potential φ becomes
κ
−1
∂
2
t
φ =∇·(ρ
−1
∇φ) + f, (9)
where κ denotes the bulk modulus. It follows that velocity v and
pressure p may be expressed as
v = ρ
−1
∇∂
t
φ, (10)
p =−κ ( ∇·s) =−∂
2
t
φ. (11)
The resulting formulation for pressure p is the reason why we choose
to define the potential φ as in eq. (8). Since pressure is continuous
across first-order discontinuities, it follows that ∂
2
t
φ and thus φ must
be continuous, a requirement which is honoured automatically by
the basis functions of the SEM. The source f may be expressed in
terms of pressure P acting at location x
s
.
f =−κ
−1
P(t ) δ(x − x
s
) . (12)
Note that the source is multiplied by a factor κ
−1
due to the formu-
lation used in eq. (9).
Using Gauss’ theorem and a scalar test function w , the weak
form becomes
κ
−1
w ∂
2
t
φ d
3
x =
∂
ρ
−1
w
ˆ
n ·∇φ d
2
x
−
ρ
−1
∇w ·∇φ d
3
x − κ
−1
P(t )w(x
s
).
(13)
C
2011 The Authors, GJI, 186, 721–739
Geophysical Journal International
C
2011 RAS
Downloaded from https://academic.oup.com/gji/article/186/2/721/590417 by guest on 18 June 2021
724 D. Peter et al.
Figure 1. Workflow for running spectral-element simulations with
SPECFEM3D Version 2.0 ‘Sesame’.
At the free surface ∂ we set the pressure p =−∂
2
t
φ = 0, thereby
enforcing φ = 0, ∂
t
φ = 0and∂
2
t
φ = 0, that is, we implement
a Dirichlet boundary condition along the surface. At a fluid–solid
boundary, the interface coupling integral may be used to exchange
the normal component of displacement between fluid and solid:
ρ
−1
ˆ
n ·∇φ =
ˆ
n · s
solid
.
3 MESHING, MESH PARTITIONING
AND LOAD BALANCING
The first step in a SEM consists of constructing a high-quality mesh
for the region of interest. In this section, we outline the key issues
based on various 3-D examples. Fig. 1 draws the schematic workflow
from meshing and partitioning to finally running spectral-element
simulations. We discuss each phase separately, focussing on the use
of an external mesher, in our case CUBIT (Blacker et al. 1994).
3.1 Hexahedral meshing
We subdivide the model volume into a set of non-overlapping,
hexahedral elements. We impose that the discretization creates a
conforming mesh, that is, elements match on a full face or edge,
and the mesh cannot be discontinuous. Using the SEM with hex-
ahedral elements leads to computational benefits over tetrahedral
finite elements (Komatitsch et al. 2001; Mercerat et al. 2006; Vos
et al. 2010). Especially for parallel implementations, taking advan-
tage of the diagonal mass matrix and optimized tensor products is
critical in terms of computational speed (Komatitsch et al. 2003;
Carrington et al. 2008; Vos et al. 2010). Hexahedral meshing is also
attractive for the SEM because it benefits from reduced errors and
generally smaller element counts compared to tetrahedral meshing
(Hesthaven & Teng 2000; Komatitsch et al. 2001; Vos et al. 2010).
Unfortunately, automatic 3-D hexahedral mesh generation
is more demanding than unstructured tetrahedral meshing
(Shepherd & Johnson 2008; Staten et al. 2010). To construct hex-
ahedral meshes, our examples make use of an external hexahedral
mesher, such as CUBIT (Blacker et al. 1994). We focus on this
particular mesh generation tool kit because it is a well-documented
and feature-rich package, on which most of our own experience is
based. One may readily use other meshing tools, such as Abaqus
(SIMULIA 2008), ANSYS (ANSYS 2011), GOCAD (Mallet 1992;
Caumon et al. 2009), GiD (Gardia-Donoro et al. 2010; Rib
´
o et al.
2011), Gmsh (Geuzaine & Remacle 2009), TrueGrid (Noble &
Nuss 2004; Rainsberger 2006) or Salome (Ribes & Caremoli 2007;
Bergeaud et al. 2010).
Fig. 2 shows several examples of fully unstr uctured hexahe-
dral meshes. In the Mount St Helens region, the mesh employs
a mesh tripling layer to increase resolution at the topographic sur-
face. Tripling is the default refinement in CUBIT for subdividing
hexahedral elements in a conforming fashion. Surface topography is
imported using Shuttle Radar Topographic Mission (SRTM) data,
converted to Universal Transverse Mercator (UTM) coordinates
with an original resolution of 90 m (Jarvis et al. 2008). Meshing is
performed automatically by CUBIT using a sweep algorithm. The
resolution of the mesh enables seismic wave simulations with fre-
quencies up to ∼1.5 Hz. The Mesh for the L’Aquila region, Italy,
consists of ∼7 M hexahedra with an element size of ∼90 m at
the top surf ace. This mesh facilitates simulations of seismic wave
propagation up to ∼5 Hz. For the exploration geophysics model,
the hexahedral mesh honours a salt dome body inside a 3-D model
capped by a water layer. The mesh for asteroid 433-Eros with a
close-bound surface has a resolution of roughly 300 m. Finally, the
filled coffee cup model discretized into hexahedra couples an elastic
domain for the cup with an acoustic domain for the coffee inside
the cup.
To ensure compatibility with previous versions of SPECFEM3D
(see e.g. Komatitsch et al. 2004; Liu et al. 2004), the in-house
mesher based on analytical linear interpolation from the top to the
bottom of the mesh has been adapted to the new code structure. It
facilitates the design of simpler, alternative meshes for layercake
models.
3.2 Partitioning and load balancing
Balancing the computational load and distributing the mesh on a
large number of cores is crucial for optimized high-performance
simulations (Martin et al.
2008a). To do so, we make use of an
e
xternal partitioner, namely SCOTCH (Pellegrini & Roman 1996;
Chevalier & Pellegrini 2008), which we use to balance spectral-
element computations on an arbitrary number of cores. An alter-
native partitioner able to fulfill these tasks is METIS (Karypis &
Kumar 1998), b ut SCOTCH is more actively maintained (Chevalier
& Pellegrini 2008) and performs better in many cases that we have
tested.
Especially for s imulations involving coupled elastic and acoustic
domains, balancing the mesh becomes paramount. Most of the com-
putation time is spent resolving the divergence of the stress tensor
in each element. The computational cost for an elastic element is
approximately four times larger than for an acoustic element, which
may be established by running simulations for one domain at a time.
During partitioning, we therefore weight each element according to
its associated domain type and computational cost to balance the
overall numerical cost rather than simply the number of elements
between partitions. The major improvement in SPECFEM3D code
performance focuses on these tensor products, using highly effi-
cient algorithms developed by Deville et al. (2002) and optimizing
cache usage. Another key aspect of mesh partitioning is minimiza-
tion of the number of edge cuts, because this reduces the amount
of MPI communications between processor cores (an edge cut oc-
curs when two contiguous elements are assigned to distinct cores).
On machines comprising a very large number of cores, it is cru-
cial to resort to non-blocking communications between compute
nodes, for instance using non-blocking MPI message passing, to
obtain good performance scaling (Danielson & Namburu 1998;
Komatitsch et al. 2008; Martin et al. 2008a).
C
2011 The Authors, GJI, 186, 721–739
Geophysical Journal International
C
2011 RAS
Downloaded from https://academic.oup.com/gji/article/186/2/721/590417 by guest on 18 June 2021