# Fast Multipole accelerated boundary element method for elastic wave propagation in multi-region domains

Stéphanie Chaillat, Marc Bonnet, Jean-François Semblat

► To cite this version:

Stéphanie Chaillat, Marc Bonnet, Jean-François Semblat. Fast Multipole accelerated boundary element method for elastic wave propagation in multi-region domains. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422242

HAL Id: hal-01422242

https://hal.science/hal-01422242

Submitted on 24 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific 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 diffusion de documents
scientifiques 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fast Multipole accelerated boundary element method for elastic wave propagation in multi-region domains

S. Chaillat<sup>1,2,*</sup>, M. Bonnet<sup>1</sup>, J.F. Semblat<sup>2</sup>

<sup>1</sup> Laboratoire de Mécanique des Solides (UMR 7649 CNRS)
Ecole Polytechnique, 91128 Palaiseau Cedex
bonnet@lms.polytechnique.fr

<sup>2</sup> Division MSRGI, LCPC
58 bd Lefebvre, 75732 Paris Cedex 15, FRANCE
semblat@lcpc.fr

**Abstract** — Solving the 3-D elastodynamic equations using traditional boundary element methods (BEMs) is greatly hindered by the fully-populated nature of the BEM matrix. In a previous study limited to homogeneous media, we have established that the Fast Multipole (FM) method reduces the complexity of a 3-D elastodynamic BEM to $N \log N$ per GMRES iteration. Here, the methodology is extented to piecewise-homogeneous domains using a FM-accelerated multi-region BE-BE coupling. Numerical examples, and a simple preconditioning approach, are presented.

**Keywords** — Fast Multipole Method, Multi-region problems, Elastodynamics.

## 1 Introduction

Numerical methods proposed to date for computing seismic wave propagation in alluvial basins exploit series expansions, multipolar expansions of wave functions, finite elements, finite differences, spectral elements, or boundary elements, each approach having specific advantages and limitations. The main advantage of the boundary element method (BEM) is that only the domain boundaries (and possibly interfaces) are discretized, leading to a reduction of the number of degrees of freedom (DOFs), and avoiding cumulative effects of grid dispersion. However, the standard BEM leads to fully-populated matrices, which results in high computational costs in CPU time ($O(N^2)$) per iteration using an iterative solver such as GMRES) and memory requirements ($O(N^2)$), where $N$ denotes the number of DOFs of the BEM model. The appearance of accelerated boundary element (BE) methodologies, and especially the rapid development of the Fast Multipole Method (FMM) over the last 10-15 years, has dramatically improved the capabilities of BEMs for many areas of application through solution speedup, memory requirement reduction, and model size increase. To date, only few studies on the FMM in elastodynamics (including [5] for the frequency-domain case) are available. In [1], the FMM for homogeneous semi-infinite elastic propagation domains is improved by incorporating recent advances of FMM implementations for Maxwell equations [3], allowing to run BEM models of size up to $N = O(10^6)$ on a single-processor PC. The present contribution aims at extending the formulation of [1] to multi-domain situations, with emphasis on alluvial-basin configurations, by developing a FMM-based BE-BE coupling approach suitable for 3-D piecewise-homogeneous media.

* Current address: College of Computing, Georgia Institute of Technology, 266 Ferst Drive, Atlanta GA 30332-0765, USA (stephanie.chaillat@cc.gatech.edu)
## 2 Standard and fast multipole accelerated boundary element method

**Single-region boundary element method.** Let $\Omega$ denote a region of space occupied by an isotropic elastic solid characterized by $\mu$ (shear modulus), $\nu$ (Poisson's ratio) and $\rho$ (mass density). A time-harmonic motion with circular frequency $\omega$ is assumed, and the implicit factor $e^{-i\omega t}$ will be systematically omitted. Assuming the absence of body forces, the displacement $\mathbf{u}$ is governed by the well-known boundary integral equation (BIE)

$$c_{ik}(\mathbf{x})u_i(\mathbf{x}) + (\text{P.V.}) \int_{\partial\Omega} u_i(\mathbf{y}) T_i^k(\mathbf{x}, \mathbf{y}; \omega) dS_y - \int_{\partial\Omega} t_i(\mathbf{y}) U_i^k(\mathbf{x}, \mathbf{y}; \omega) dS_y = 0 \quad (1)$$

where $\mathbf{t}$ is the traction vector on the boundary $\partial\Omega$ with outward unit normal $\mathbf{n}(\mathbf{y})$, (P.V.) indicates a Cauchy principal value (CPV) singular integral, $U_i^k(\mathbf{x}, \mathbf{y}; \omega)$ and $T_i^k(\mathbf{x}, \mathbf{y}; \omega)$ denote the $i$-th components of the elastodynamic fundamental solution, i.e. of the displacement and traction, respectively, generated at $\mathbf{y} \in \mathbb{R}^3$ by a unit point force applied at $\mathbf{x} \in \mathbb{R}^3$ along the direction $k$:

$$\begin{align}
U_i^k(\mathbf{x}, \mathbf{y}; \omega) &= \frac{1}{4\pi k_S^2 \mu} \left( (\delta_{qs}\delta_{ik} - \delta_{qk}\delta_{is}) \frac{\partial}{\partial x_q} \frac{\partial}{\partial y_s} G(|\mathbf{y} - \mathbf{x}|; k_S) + \frac{\partial}{\partial x_i} \frac{\partial}{\partial y_k} G(|\mathbf{y} - \mathbf{x}|; k_P) \right), \\
T_i^k(\mathbf{x}, \mathbf{y}; \omega) &= \mu \left[ \frac{2\nu}{1-2\nu} \delta_{ij}\delta_{k\ell} + \delta_{ik}\delta_{j\ell} + \delta_{jk}\delta_{i\ell} \right] \frac{\partial}{\partial y_\ell} U_h^k(\mathbf{x}, \mathbf{y}; \omega) n_j(\mathbf{y}), \tag{2} \\
k_S^2 &= \rho\omega^2/\mu, \quad 2(1-\nu)k_P^2 = (1-2\nu)k_S^2
\end{align}$$

in which $G(r;k) = \exp(ikr)/(4\pi r)$ is the free-space Green's function for the Helmholtz equation with wavenumber $k$ corresponding to either $P$ or $S$ elastic waves, and the free-term $c_{ik}(\mathbf{x})$ is equal to $0.5\delta_{ik}$ in the usual case where $\partial\Omega$ is smooth at $\mathbf{x}$. The numerical solution of BIE (1) is based on a boundary element (BE) discretization of the surface $\partial\Omega$ and boundary traces $(\mathbf{u}, \mathbf{t})$, leading to the system:

$$[K]\{\nu\} = \{f\}, \quad (3)$$

where the $N$-vector $\{\nu\}$ collects the unknown degrees of freedom (DOFs).

**Fast Multipole Method: principle.** BEM matrix equations such as (3) are here solved iteratively using GMRES. The influence matrix $[K]$ is fully-populated. Each GMRES iteration requires one evaluation of $[K]\{\nu\}$ for given $\{\nu\}$, a task requiring a computing time of order $O(N^2)$ using traditional BEM techniques. To lower this $O(N^2)$ complexity, unacceptable for large BEM models, fast BEM solutions techniques such as the Fast Multipole Method (FMM) must be employed.

In general terms, the FMM exploits a reformulation of the fundamental solutions in terms of products of functions of $\mathbf{x}$ and of $\mathbf{y}$, so that (unlike in the traditional BEM) integrations with respect to $\mathbf{y}$ can be reused when the collocation point $\mathbf{x}$ is changed. Using the so-called multi-level FMM allows to reduce the complexity of a GMRES iteration to $O(N \log N)$; moreover, the governing BEM matrix is never explicitly formed, which leads to a storage requirement well below the $O(N^2)$ memory required for holding it. The reader is referred to [1] for details on the FMM and its implementation for single-domain elastodynamic problems.

## 3 FM-accelerated BE-BE coupling

**Continuous BEM formulations for seismic wave propagation.** This formulation, and its present implementation based on the multi-domain FM-accelerated BEM, are geared towards geometrical configurations involving a semi-infinite homogeneous reference medium with topographic irregularities and alluvial deposits (henceforth generically referred to as irregularities, Fig. 1).

In the following, $\Omega_F$ denotes the free half-space $\{\mathbf{x} = (x_1, x_2, x_3) | x_3 < 0\}$ bounded by the infinite traction-free surface $\Gamma_F = \{\mathbf{x} | x_3 = 0\}$. Configurations treated in this article are perturbations of the free half-space $\Omega_F$, where irregularities occur only in a region of finite size. For
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Propagation of seismic waves in complex geological structures (alluvial deposits, basins): various geometries and related notations.</figcaption></figure>

such configurations, the displacement vector $\mathbf{u}$ is split into $\mathbf{u} = \mathbf{u}^F + \mathbf{u}^S$, where $\mathbf{u}^F$ characterizes the free-field, a known seismic wave in the reference free half-space $\Omega_F$ composed of the incident waves and those reflected from the planar free surface $\Gamma_F$, so that $\mathbf{t}^F = 0$ on $\Gamma_F$. The scattered displacement $\mathbf{u}^S$ then arises due to the presence of irregularities. On any non-planar part of the free surface, one has $\mathbf{t}^S + \mathbf{t}^F = 0$.

The governing equation for the total field in $\Omega_1$ (using shorthand notations $U_i^k$ and $T_i^k$ instead of $U_i^k(\mathbf{x}, \mathbf{y}; \omega)$ and $T_i^k(\mathbf{x}, \mathbf{y}; \omega)$ for convenience) is:

$$c_{ik}(\mathbf{x})u_i(\mathbf{x}) + \int_{\Gamma_1 \cup \Gamma(D)} u_i^1(\mathbf{y})T_i^{k(1)}(\mathbf{y})dS_y + \sum_{m=2}^{n} \left( \int_{\Gamma_{1m}} u_i^{1m}(\mathbf{y})T_i^{k(1)}(\mathbf{y})dS_y \right) - \int_{\Gamma_1} t_i^1(\mathbf{y})U_i^{k(1)}(\mathbf{y})dS_y \\ - \sum_{m=2}^{n} \left( \int_{\Gamma_{1m}} t_i^{1m}(\mathbf{y})U_i^{k(1)}(\mathbf{y})dS_y \right) = c_{ik}^F(\mathbf{x})u_i^F(\mathbf{x}) + \int_{\Gamma_F(D)} u_i^F(\mathbf{y})T_i^{k(1)}(\mathbf{y})dS_y, \quad \forall \mathbf{x} \in \partial\Omega_1 \quad (4)$$

while the total field in subdomain $\Omega_\ell$ ($\ell > 1$) is governed by the integral equation:

$$c_{ik}(\mathbf{x})u_i(\mathbf{x}) + \int_{\Gamma_{\ell}} u_i^{\ell}(\mathbf{y})T_i^{k(\ell)}(\mathbf{y})dS_y + \sum_{m=2}^{\ell-1} \int_{\Gamma_{lm}} (u_i^{ml}(\mathbf{y})T_i^{k(\ell)}(\mathbf{y}) + t_i^{ml}(\mathbf{y})U_i^{k(\ell)}(\mathbf{y}))dS_y \\ + \sum_{m=\ell+1}^{n} \int_{\Gamma_{lm}} (u_i^{ml}(\mathbf{y})T_i^{k(\ell)}(\mathbf{y}) - t_i^{ml}(\mathbf{y})U_i^{k(\ell)}(\mathbf{y}))dS_y = 0, \forall \mathbf{x} \in \partial\Omega_\ell, (2 \le \ell \le n) \quad (5)$$

In equations (4) and (5), $U_i^{k(\ell)}$ and $T_i^{k(\ell)}$ denote the fundamental solutions defined in terms of the material parameters of $\Omega_\ell$, and $u_i(\mathbf{x})$ in the free-term stands for either $u_i^l(\mathbf{x})$ or $u_i^{lm}(\mathbf{x})$, according to whether the collocation point $\mathbf{x}$ lies on $\Gamma_\ell$ or $\Gamma_{lm}$. Moreover, use has been made of the free surface condition on $\Gamma_\ell$ and the perfect-bonding transmission conditions $\mathbf{u}^{lm} = \mathbf{u}^{m\ell}$, $\mathbf{t}^{lm} = -\mathbf{t}^{m\ell}$ to establish equations (5) [2].

The coupled BE-BE formulation to be presented next will then be based on combining discrete versions of equation (4) and equations (5) written for each subregion $\Omega_\ell$ ($\ell \ge 2$). It is similar to the one used for two subdomains in [5], but more general as (i) it is applicable to an arbitrary number of subdomains and (ii) it accommodates irregularities going above or through the free surface.

**BE-BE coupling strategy.** The present discrete coupled BE-BE formulation results from combining discrete versions of equation (4) and equations (5) written for each subregion $\Omega_\ell$ ($\ell \ge 2$). It is similar to the one used for two subdomains in [5], but more general as (i) it is applicable to an arbitrary number of subdomains and (ii) it accommodates irregularities going above or through the free surface. Its present implementation is based on three-noded triangular BEs, piecewise-linear interpolation of displacements, and piecewise-constant interpolation of tractions. Since only Neumann or transmission boundary conditions are considered here, the displacement is unknown at all mesh nodes, while the traction is unknown on each interfacial element. The chosen "element-based" traction interpolation permits traction discontinuities across edges. This is particularly convenient when the latter involve "triple points" shared by three (or more) subregions, whereas "node-based"
traction modelling would entail cumbersome adjustments due to the multiplicity of tractions from adjacent faces at such points. The proposed BE-BE coupling formulation is designed so as to invoke single-region FM-BEM computations in "black-box" fashion (here using the single-region elastodynamic FM-BEM of [1]). To this end, a boundary integral equation is formulated for each subregion $\Omega_i$ (with material properties assumed homogeneous in each $\Omega_i$), and discrete BE equations are generated by using (i) all displacement nodes and (ii) all interfacial element centers as collocation points (i) and (ii) will subsequently be referred to as "nodal collocation" and "element collocation", respectively). Each subregion is treated separately, using a separate octree for FMM computations. The matrix-vector products arising in each of these integral equations can thus be evaluated using the FM-BEM procedure for homogeneous media of [1].

The BE-BE coupling does not, however, just consist of concatenating all single-region BE equations into one global system of equations, as the latter would be overdetermined as a result. To ensure that the present BE-BE coupling yields a square global system of equations, judiciously chosen linear combinations of BE equations generated at the subregion level, arising from collocation at (a) interfacial element centers relative to either subregion adjacent to that element and (b) displacement nodes shared by more than one subregion, are defined (see [2] for details). This treatment is done externally to the FM-BEM computations, and handles easily cases of multiple displacement nodes (e.g. triple points in the case of a two-layered basin).

## 4 Propagation and amplification of seismic waves in alluvial basins

In [1], the single-domain elastodynamic FMM has been compared to the results of [6] for the scattering by an irregular homogeneous half-space of a plane vertical P-wave at normalized frequency $k_P a/\pi = 0.25$ (with $v = 0.25$), and then applied to the same configuration at a higher frequency ($k_P a/\pi = 5$). Here, the present multi-domain implementation is applied to an example where a plane vertically-incident P-wave of unit amplitude is scattered by a semi-spherical alluvial basin. This example has been run on a single-processor PC (RAM: 3GB, CPU frequency: 3.40 GHz). As in [6], we investigate the motion at the surface of the alluvial basin $\Omega_2$, for the following values of the material parameters: $\mu^{(2)} = 0.3\mu^{(1)}$, $\rho^{(2)} = 0.6\rho^{(1)}$, $v^{(1)} = 0.25$ and $v^{(2)} = 0.3$. The normalized frequency is defined by $k_P^{(1)} a/\pi$ in terms of the properties of the elastic semi-infinite medium $\Omega_1$. The radius of the discretized free surface is set to $D = 5a$.

The surface displacements computed with the present multi-domain FMM are presented, along with corresponding results from [6] (using series expansion method) and [4] (using spectral element method), for $k_P^{(1)} a/\pi = 0.5$ (Fig. 3a). All results are seen to be in good agreement. Additionally, the FMM allowed to perform computations at higher frequency $k_P^{(1)} a/\pi = 2$ (Fig. 3b), for which no published results are available for comparison purposes. The number of DOFs, the size of the leaf cells and the leaf level $\bar{l}_i$ in each subdomain $\Omega_i$ are given in the table below for the meshes used, together with the CPU time per iteration recorded. The iteration count is seen to significantly impact the computational efficiency for problem sizes for which the CPU time per iteration and the memory requirements are still moderate. An efficient preconditioning strategy is clearly needed.

<table><thead><tr><th>$k_P^{(1)} a/\pi$</th><th>$N$</th><th>$d^{\min}/\lambda_s$</th><th>$\bar{l}_1; \bar{l}_2$</th><th>CPU (s) / iter</th><th>nb iter.</th></tr></thead><tbody><tr><td>0.5</td><td>17,502</td><td>0.15</td><td>3;3</td><td>8</td><td>28</td></tr><tr><td>2</td><td>190,299</td><td>0.30</td><td>5;4</td><td>79</td><td>325</td></tr></tbody></table>

## 5 Preconditioning strategy

The main limiting factor for the size of the studied examples was the very high iteration counts reached, rather than the CPU time per iteration or the memory requirement. A preconditioning
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Propagation of an incident plane P-wave in a semi-spherical alluvial basin: notations.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Propagation of an incident plane P-wave in a semi-spherical alluvial basin: surface displacement at (a) $k_p^{(1)}a/\pi = 0.5$ and comparisons with [6] and [4] and (b) $k_p^{(1)}a/\pi = 2$.</figcaption></figure>

strategy is clearly needed to improve convergence properties for the larger models.

A feature of the FMM is that the complete system matrix is not explicitly assembled, the only explicitly available matrix being the matrix $K^{\text{near}}$ into which the near contributions are assembled. In a first step towards the development of an efficient preconditioning strategy, a nested GMRES solver has been implemented, where the preconditioning linear system based on $M = K^{\text{near}}$ used as right preconditioner is solved by the inner GMRES. Since the sparse matrix $K^{\text{near}}$ is precomputed and stored, this preconditioning strategy is efficient in terms of CPU time.

The efficiency of this preconditioning strategy has been tested on various problems involving the scattering of plane waves by canyons or alluvial basins. One such example considers the scattering of an oblique ($\theta = 30^\circ$) incident plane P-wave by a semi-spherical basin of radius $a$ (with mechanical parameters given by $v^{(1)} = 0.25$, $\mu^{(2)} = 0.3\mu^{(1)}$, $\rho^{(2)} = 0.6\rho^{(1)}$, $v^{(2)} = 0.3$). The free surface lies inside a disk of radius $D = 5a$ and the mesh features $N = 190,299$ DOFs. The non-dimensional frequency is set to $k_p^{(1)}a/\pi = 2$. The relative tolerance is set to $\varepsilon_{\text{inner}} = 10^{-1}$ for the inner solver and $\varepsilon_{\text{outer}} = 10^{-3}$ for the outer solver. The number of outer iterations is found to be greatly reduced (from 388 to 26). The preconditioning strategy additionally involves 231 inner iterations, which also influence the overall efficiency as they occur only in the preconditioned version. The cumulative CPU time was found to be reduced from 7h59'27'' to 2h30'54''.

# 6 Diffraction of a vertical incident plane P-wave by an Alpine valley

The numerical efficiency of the present FM-BEM has been shown on canonical examples. Now, the method is applied to a more realistic seismological application, namely the diffraction of a vertical incident plane P-wave by an Alpine valley (Grenoble). In the bedrock, denoted $\Omega_1$, the P- and S-velocities and mass density are set to $c_p^{(1)} = 5,600 \text{ m.s}^{-1}$, $c_s^{(1)} = 3,200 \text{ m.s}^{-1}$ and $\rho^{(1)} = 2,720 \text{ kg.m}^{-3}$. The sedimentary basin $\Omega_2$ is modelled here with just one single homogeneous
layer, with $c_p^{(2)} = 1,988 \text{ m.s}^{-1}$, $c_s^{(2)} = 526 \text{ m.s}^{-1}$ and $\rho^{(2)} = 2,206 \text{ kg.m}^{-3}$.

The diffraction of a vertical incident plane P-wave by the valley is considered for $f = 0.6$ Hz. For this example, the number of DOFs is $N = 141,288$, the leaf levels are $\bar{l}_1 = 5$ and $\bar{l}_2 = 6$. The CPU time per iteration (without preconditioning) is 77 s, the number of iterations is 747 and the cumulative CPU time (with preconditioning strategy) is 75h45'44''. In Figure 4, the modulus of the z- surface displacement component is displayed. This example shows the possibility of very high amplifications inside the aluvial basin. It also underlines the current limitation of the present FM-BEM to deal with basin problems featuring a high velocity contrast between two layers, with a loss of efficiency caused by highly non-uniform meshes.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Propagation of a vertical incident plane P-wave in the Alpine valley: modulus of the z-component of displacement for frequency $f = 0.6$ Hz.</figcaption></figure>

# 7 Conclusion

In this communication, a multi-domain fast multipole formulation has been proposed, based on a previous single-region FMM [1] and a BE-BE coupling strategy. Comparisons with previously published numerical results demonstrate the accuracy of the present implementation. The analysis of seismic wave propagation in canonical basins, for frequencies higher than in previously published results, show the numerical efficiency of the method and its suitability for dealing with realistic seismological applications. A simple and efficient preconditioning strategy has been proposed, and its efficiency shown on basin problems. Further study is expected to bring improvements on this issue. Other ongoing work deals with the formulation of multipole expansions of the half-space fundamental solutions and with the FMM for attenuating media.

# References

[1] S. Chaillat, M. Bonnet and J.F. Semblat. A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain. *Comput. Methods Appl. Mech. Engng.*, **197**:4233–4249, 2008.

[2] S. Chaillat, M. Bonnet and J.F. Semblat. A new fast multi-domain BEM to model seismic wave propagation and amplification in 3D geological structures. *Geophys. J. Int.*, 2009 (in press).

[3] E. Darve. The fast multipole method: numerical implementation. *J. Comp. Phys.*, **160**:195–240, 2000.

[4] E. Delavaud. Simulation numérique de la propagation d'ondes en milieu géologique complexe: application à l'évaluation de la réponse sismique du bassin de Caracas (Venezuela), Ph.D. thesis, Institut de Physique du Globe de Paris, 2007.

[5] H. Fujiwara, H. The fast multipole method for solving integral equations of three-dimensional topography and basin problems. *Geophys. J. Int.*, **140**:198–210, 2000.

[6] F.J. Sánchez-Sesma. Diffraction of elastic waves by 3D surface irregularities, *Bull. seism. Soc. Am.*, **73**:1621–1636, 1983.