# Fictitious domain method for an inverse problem in volcanoes

Olivier Bodart, Valérie Cayol, Farshid Dabaghi, Jonas Koko

► To cite this version:

Olivier Bodart, Valérie Cayol, Farshid Dabaghi, Jonas Koko. Fictitious domain method for an inverse problem in volcanoes. Domain Decomposition Methods in Science and Engineering XXV Series: Lecture Notes in Computational Science and Engineering, 138, Springer International Publishing, pp.498, 2020. [hal-02325032](https://hal HAL Id:)

HAL Id: hal-02325032

https://uca.hal.science/hal-02325032

Submitted on 13 Dec 2020

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
Fictitious domain method for an
inverse problem in volcanos

Oliver Bodart, Valérie Cayol, Farshid Dabaghi and Jonas Koko

**Abstract** An inverse problem applied to volcanology is studied. It consists in the determination of the pressure applied to a crack in order to fit observed ground deformations. The deformation of the volcano is assumed to be governed by a linear elasticity PDE. The direct problem is solved via a fictitious domain method, using a finite element discretization of XFEM type. The gradient of the cost function is derived from a Lagrangian. A numerical test is presented.

# 1 General framework and problem setting

Problems in volcanology often involve elasticity models in presence of cracks (see e.g. [3]). Most of the time the force exerted on the crack is unknown, and the position and shape of the crack are also frequently unknown or partially known (see e.g. [2]). Most of the time the model is approximated *via* boundary element methods. These methods are quite convenient to take into account the crack since the problem is then reformulated into an external problem where the crack is the only object to be meshed. However these methods do not allow to take the heterogeneity and/or the anisotropy of the medium into account. Another drawback is that, when it comes to identifying the shape

Oliver Bodart
Université Jean Monnet Saint-Étienne, CNRS UMR 5208, Institut Camille Jordan, F-42023 Saint-Etienne, France, e-mail: olivier.bodart@univ-st-etienne.fr

Valérie Cayol, Farshid Dabaghi
Laboratoire Magmas et Volcans, Université Jean Monnet-CNRS-IRD, Saint-Etienne F-42023, France e-mail: v.cayol@opgc.fr, farshid.dabaghi@univ-st-etienne.fr

Jonas Koko
LIMOS, UMR 6158, Université Clermont Auvergne, BP 10448, F-63173 Aubière Cedex, France e-mail: jonas.koko@uca.fr
and/or location of the crack, the variation of the latter implies a remeshing and assembling of all the matrices of the problem.

Using a domain decomposition technique then appears as the natural solution to these problems. In [1], a first step was made with the development of a direct solver implementing a domain decomposition method. The present work represents a step further with the use of such a solver, which has been improved since the publication of [1], to solve inverse problems in the field of earth sciences. To our knowledge, this is the first work using these kind of techniques in this field of application. The next step of our project will be the shape optimization problem to identify the shape and location of the crack.

Let $\Omega$ be a bounded open set in $\mathbb{R}^d$, $d = 2,3$ with smooth boundary $\partial\Omega := \Gamma_D \cup \bar{\Gamma}_N$ where $\Gamma_D$ and $\Gamma_N$ are of nonzero measure and $\Gamma_D \cap \Gamma_N = \emptyset$. We assume that $\Omega$ is occupied by an elastic solid and we denote by $\mathbf{u}$ the displacement field of the solid and the density of external forces by $\mathbf{f} \in \mathbf{L}^2(\Omega)$. The Cauchy stress $\sigma(\mathbf{u})$ and strain $\varepsilon(\mathbf{u})$ are given by

$$ \sigma(\mathbf{u}) = \lambda(\mathrm{Tr}\varepsilon(\mathbf{u}))\mathbf{I}_{\mathbb{R}^d} + 2\mu\varepsilon(\mathbf{u}) \quad \text{and} \quad \varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T), $$

where $(\lambda, \mu)$ are the Lamé coefficients, $\mathbf{I}_{\mathbb{R}^d}$ denotes the identity tensor, and $\mathrm{Tr}(\cdot)$ represents the matrix trace. Consider a crack $\Gamma_C \subset \Omega$ represented by a line ($d = 2$) or a surface ($d = 3$) parametrized by an injective mapping. Around the crack, $\Omega$ is split into $\Omega^-$ and $\Omega^+$. The deformation field of the solid is supposed to satisfy the following elastostatic system:

$$ \left\{ \begin{array}{ll} \text{Find } \mathbf{u} \in \mathbf{H}^1(\Omega \setminus \Gamma_C) \text{ such that :} & \\ -\mathrm{div} \, \sigma(\mathbf{u}) = \mathbf{f} & \text{in } \Omega \setminus \Gamma_C, \\ \mathbf{u} = 0 & \text{in } \Gamma_D, \\ \sigma(\mathbf{u}) \cdot \mathbf{n} = 0 & \text{on } \Gamma_N, \\ \sigma(\mathbf{u}) \cdot \mathbf{n}^{\pm} = p\mathbf{n}^{\pm} & \text{on } \Gamma_C. \end{array} \right. \qquad (1) $$

where $\mathbf{n}$ is the outward unit normal to its boundaries. Typically in such a situation, $\Gamma_N$ is the ground surface and free to move. Practically, the displacement field can be observed on $\Gamma_N$, whereas the pressure $p$ exerted on the crack is unknown most of the time.

Consider the following function defined on $\mathbf{L}^2(\Gamma_C)$:

$$ J(p) := \frac{1}{2} \int_{\Gamma_N} (\mathbf{u} - \mathbf{u}_d)^{-1} (\mathbf{u} - \mathbf{u}_d)^T d\Gamma_N + \frac{\alpha}{2} \|p\|^2_{\mathbf{L}^2(\Gamma_C)}, \qquad (2) $$

where $u_d \in \mathbf{L}^2(\Gamma_N)$ is the measured displacement field and $u$ is the solution of (1) associated with $p$. Moreover, the matrix $C$ is the covariance operator of the measurements uncertainties, and is assumed to be positive definite (see e.g. [5]), and finally $\alpha > 0$ is a regularization parameter. The aim of this work is to study the following problem, of optimal control type:
$$ \min_{p \in L^2(\Gamma_C)} J(p). \qquad (3) $$

The paper is organized as follows : the next section will be devoted to the presentation of the domain decomposition method and its discretization. Section 3 gives the optimality conditions for the problem (3) and establishes their discrete version. A special focus will be made on the adaptation of the problem to a domain decomposition formulation. Finally we present a relevant numerical test in section 4 and discuss the next steps of our project.

## 2 Domain decomposition : the direct solver

To solve the direct problem (1), we use a domain decomposition method. More precisely, following [1], the domain $\Omega$ is split into two subdomains such that each point of the domain lies on one side of the crack or on the crack. Moreover, the global unknown solution $u$ is decoupled in two sub-solutions for each side of the crack. For this purpose, we are using an artificial extension of the considered crack $\Gamma_C$ (e.g. $\Gamma_0$ in Figure 1). Therefore, instead of the crack

<figure><img src="image_2.png" /><figcaption>Fig. 1 Splitting the volcanic cracked domain</figcaption></figure>

problem (1), we have to solve two Neumann-type boundary problems such that for each problem we impose a pressure on $\Gamma_C$, which is more convenient from both theoretical and numerical points of view. More precisely we solve the following system:
$$
\left\{
\begin{array}{ll}
\text{Find } \mathbf{u} \in \mathbf{H}^1(\Omega) \text{ such that :} & \\
-\operatorname{div} \sigma(\mathbf{u}^{\pm}) = \mathbf{f}^{\pm} & \text{in } \Omega^{\pm}, \\
\mathbf{u}^{\pm} = 0 & \text{on } \Gamma_{\mathrm{D}} \cap \partial\Omega^{\pm}, \\
(\sigma(\mathbf{u}) \cdot \mathbf{n})^{\pm} = 0 & \text{on } \Gamma_{\mathrm{N}} \cap \partial\Omega^{\pm}, \\
(\sigma(\mathbf{u}) \cdot \mathbf{n})^{\pm} = p\mathbf{n}^{\pm} & \text{on } \Gamma_{\mathrm{C}}, \\
[\mathbf{u}] = 0 & \text{on } \Gamma_0, \\
[\sigma(\mathbf{u})] \cdot \mathbf{n}^{+} = 0 & \text{on } \Gamma_0,
\end{array}
\right.
\qquad (4)
$$

where $\mathbf{u}^+ = \mathbf{u}|_{\Omega^+}$ and $\mathbf{u}^- = \mathbf{u}|_{\Omega^-}$, and $[\mathbf{v}]$ denotes the jump of $v$ across $\Gamma_0$. The two last conditions in (4) enforce the continuity of displacement and stress across $\Gamma_0$. Notice that the boundary conditions on $\Gamma_0$ ensure the construction of a global displacement field in $\mathbf{H}^1(\Omega \setminus \Gamma_C)$ solving the original problem (1).

Let us define the following Hilbert spaces

$$
\mathbf{V}^{\pm} = \{v \in \mathbf{H}^{1}(\Omega^{\pm}) \mid v = 0 \text{ on } \Gamma_{\mathrm{D}} \cap \partial\Omega^{\pm}\}, \quad \mathbf{W} = (\mathbf{H}^{\frac{1}{2}}(\Gamma_{0})).
$$

and their dual spaces $\mathbf{V}'^{\pm}$ and $\mathbf{W}'$, endowed with their usual norms. Pre-
scribing the continuity of displacement across the $\Gamma_0$ via a Lagrangian for-
mulation, the mixed weak formulation of Problem (4) reads as follows:

$$
\left\{
\begin{array}{ll}
\text{Find } \mathbf{u}^{\pm} \in \mathbf{V}^{\pm} \text{ and } \boldsymbol{\lambda} \in \mathbf{W}' \text{ such that :} & \\
\quad a(\mathbf{u}^{\pm}, \mathbf{v}^{\pm}) \pm b(\boldsymbol{\lambda}, \mathbf{v}^{\pm}) = l^{\pm}(\mathbf{v}^{\pm}) & \quad \forall \mathbf{v}^{\pm} \in \mathbf{V}^{\pm}, \\
\quad b(\boldsymbol{\mu}, [\mathbf{u}]) = 0 & \quad \forall \boldsymbol{\mu} \in \mathbf{W}',
\end{array}
\right.
\qquad (5)
$$

with

$$
a(\mathbf{u}^{\pm}, \mathbf{v}^{\pm}) = \int_{\Omega^{\pm}} \sigma(\mathbf{u}^{\pm}) : \varepsilon(\mathbf{v}^{\pm}) \, d\Omega^{\pm}
$$

bilinear, symmetric, coercive and

$$
l^{\pm}(\mathbf{v}^{\pm}) = \int_{\Omega^{\pm}} f \cdot \mathbf{v}^{\pm} \, d\Omega^{\pm} + \int_{\Gamma_C} (p\mathbf{n})^{\pm} \cdot \mathbf{v}^{\pm} \, d\Gamma_C
$$

linear and continuous. Moreover, *b* is defined as the duality pairing between **W**' and **W**: *b*(*λ*, **v**^±) = ⟨*λ*, **v**^±⟩_{**W**',**W**}. Therefore, it is straightforward to prove the existence and uniqueness of a solution to Problem (5) (see e.g. [1] ans references within).

Denoting then **f**± and **p** the approximations of *f* and *p* in **V**± and **W**̂h,
setting *p*<sub>*n*</sub> = **p****n**<sup>+</sup> = −**p****n**<sup>−</sup> and
$$
\begin{align*}
K &= \begin{pmatrix} A^+ & 0 & B^{+T} \\ 0 & A^- & -B^{-T} \\ B^+ & -B^- & 0 \end{pmatrix}, &
\mathbf{X} &= \begin{pmatrix} \mathbf{u}^+ \\ \mathbf{u}^- \\ \boldsymbol{\lambda} \end{pmatrix}, \\
\mathbf{F} &= \begin{pmatrix} \mathbf{F}^+ \\ \mathbf{F}^- \\ 0 \end{pmatrix} = \begin{pmatrix} M_{\Omega}^+ \cdot \mathbf{f}^+ \\ M_{\Omega}^- \cdot \mathbf{f}^- \\ 0 \end{pmatrix} + \begin{pmatrix} +M_c^+ \cdot p_n \\ -M_c^- \cdot p_n \\ 0 \end{pmatrix} := L_{\Omega} \cdot \mathbf{f} + L_c \cdot \mathbf{p},
\end{align*}
$$

the discretized form of system (5) has the linear algebraic formulation

$$
KX = F. \tag{6}
$$

The system (6) can be solved by a Uzawa Conjugate gradient/domain de-
composition method [1]. The method can be classically stabilized and the
convergence of the numerical scheme can be proved as $h \to 0$.

In what follows, we will focus on the adaptation of a crack inverse problem
to this domain decomposition formulation and its application to a realistic
problem.

3 The crack inverse Problem

First, we have the following result.

**Proposition 1** For any $\alpha > 0$, the problem (3) admits a unique solution $p^*$ in $L^2(\Gamma_C)$.

*Proof* The proof is classical: applying the same method as in [4], one easily shows that *J* is strictly convex and coercive on *L*²(*Γ*<sub>C</sub>). □

The objective function *J* being strictly convex, first order optimality condi-
tions can be computed to implement an suitable optimization method (in our
case the conjugate gradient). Let us introduce the adjoint system

$$
\left\{
\begin{array}{ll}
-\operatorname{div} \sigma(\phi) = 0 & \text{in } \Omega, \\
\phi = 0 & \text{on } \Gamma_{\mathrm{D}}, \\
\sigma(\phi) \cdot \mathbf{n} = C^{-1}(\mathbf{u} - \mathbf{u}_d) & \text{on } \Gamma_{\mathrm{N}},
\end{array}
\right.
\qquad (7)
$$

where **u** is a solution of system (1). Is is easy to prove that this adjoint system
admits a unique solution **φ** ∈ H¹(Ω). We have the following

**Proposition 2** Let $p^*$ be the solution of problem (3) and $(\mathbf{u}^*, \phi^*)$ be the associated solutions of (1) and (7). Then, the following optimality condition holds:

$$
\alpha p^* + (\phi^* \cdot \mathbf{n}^\pm) = 0. \tag{8}
$$

This result can be proved using a classical sensitivity analysis technique. The
important point here is that it gives a way to compute the gradient of the
function $J$: for a given $p \in L^2(\Gamma_C)$, compute $(\mathbf{u}, \phi)$ which solve (1) and (7). Then, the Gâteaux derivative of $J$ is given in $L^2(\Gamma_C)$ by

$$J'(p) = \alpha p + (\phi \cdot \mathbf{n}^{\pm}), \qquad (9)$$

For a given pressure $p \in L^2(\Gamma_C)$, the computation of the gradient $J'(p)$ then requires to solve two systems.

Since we transformed our direct problem into system (4), we now need to adapt the inverse problem to this formulation. The cost function $J$ defined by (2) then rewrites into

$$J(p) := \frac{1}{2} \int_{\Gamma_N^{\pm}} (\mathbf{u}^{\pm} - \mathbf{u}_d) C^{-1} (\mathbf{u}^{\pm} - \mathbf{u}_d)^{\top} d\Gamma_N^{\pm} + \frac{\alpha}{2} \|p\|_{L^2(\Gamma_C)}^2. \quad (10)$$

Notice that the observed data $\mathbf{u}_d$ can be interpolated on two sub-domains $\Omega^{\pm}$ to obtain $\mathbf{u}_d^{\pm}$ corresponding to $\mathbf{u}^{\pm}$.

In view of (6), denoting $R$ the reduction matrix $R : \mathbf{X} \to \mathbf{U}$ and

$$\mathbf{U} = \begin{pmatrix} \mathbf{u}^{+} \\ \mathbf{u}^{-} \end{pmatrix}, \quad \mathbf{U}_d = \begin{pmatrix} \mathbf{u}_d^{+} \\ \mathbf{u}_d^{-} \end{pmatrix},$$

the discrete cost function is defined as

$$J_d(\mathbf{p}) = \frac{1}{2}(R\mathbf{X} - \mathbf{U}_d)^{\top}C^{-1}M_N(R\mathbf{X} - \mathbf{U}_d) + \frac{\alpha}{2}(\mathbf{p}^{\top}M_F\mathbf{p}), \quad (11)$$

where $\mathbf{X}$ is the solution of (6), $M_N$ and $M_F$ are the mass matrices on $\Gamma_N$ and $\Gamma_C$, respectively. This finite dimensional problem then boils down to finding the saddle point of the following Lagrangian

$$\mathcal{L}(\mathbf{X}, \mathbf{p}, \mathbf{\Phi}) = J_d(\mathbf{p}) - \langle K\mathbf{X} - (L_{\Omega}\mathbf{f} + L_c\mathbf{p}), \mathbf{\Phi} \rangle.$$

Computing the KKT conditions for this problem allows to compute the gradient of $J_d$: for a given vector $\mathbf{p}$, let $\mathbf{X}$ be the solution of (6) and $\mathbf{\Phi}$ be the solution of the adjoint problem

$$K^{\top}\mathbf{\Phi} = C^{-1}M_N(R\mathbf{X} - \mathbf{U}_d). \quad (12)$$

Then, we have

$$\nabla J_d(\mathbf{p}) = \alpha M_F \mathbf{p} + L_c^{\top} \mathbf{\Phi}, \quad (13)$$

The system (12) and the gradient (13) are the discrete counterparts of (7) and (9). As (6), the adjoint system (12) is solved by a Uzawa conjugate gradient/domain decomposition method.

**Computational aspects:** the problem studied here is actually of quadratic type. Hence it is natural to use a suitable minimization technique, namely a conjugate gradient algorithm. It is important to notice that, using the underlying quadratic form, one can determine the optimal step size. Therefore no
line search algorithm is necessary, which consequently reduces the computa-
tional cost.

4 Numerical experiments

Aiming at practical applications, we applied the technique to a realistic vol-
cano, the Piton de la Fournaise, Île de la Réunion, France. The mesh was
built from a digital elevation model (DEM), provided by the french institute
IGN (Institut Géographique National, French National Geographic Institute).
Both the boundary and volume mesh for the whole domain were generated
by Gmsh software (Figure 2, left). The crack geometry is assumed to be
quadrangular and intersecting the surface. It is constructed following [2] (see
Figure 2, left). The crack mesh does not match the volume mesh. Moreover,
it can be easily extended in order to split the domain. We assume that the
crack is submitted to an initial pressure p^0. The inverse problem will consist
in determining the unknown pressure from the surface displacements (Fig-
ure 2, right). The convergence curves in Figure 3, highlight the efficiency of
adapted optimization algorithm. The conjugate gradient minimization per-
forms efficiently, even for fine meshes.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Triangular surface mesh [2] representing the crack (left), amplitude of the displacement of a realistic volcano (right).</figcaption></figure>

5 Conclusion

We have studied a conjugate gradient type method for an interface pressure
inverse problem using a Uzawa conjugate gradient domain decomposition
method (from [1]) as inner solver. Further study is underway to derive a
single-loop conjugate gradient domain decomposition method by (directly)
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Decay of the norm of gradient (left), the error of displacement on the ground in each iteration (right) and number of iteration to the converged **p** after each refinement of the mesh (center)</figcaption></figure>

considering the constrained minimization problem (1)-(2) and using sensitiv-
ity and adjoint systems techniques.

**References**

1. Olivier Bodart, Valérie Cayol, Sébastien Court, and Jonas Koko. XFEM-based fictitious domain method for linear elasticity model with crack. *SIAM Journal on Scientific Computing*, 38(2):B219–B246, 2016.

2. Y. Fukushima, V. Cayol, and P. Durand. Finding realistic dike models from interferometric synthetic aperture radar data: The february 2000 eruption at piton de la fournaise. *Journal of Geophysical Research: Solid Earth*, 110(B3).

3. Cayol V. Kervyn F., Wauthier C. and dOreye N. Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an insar analysis. *Journal of Geophysical Research*, 117, 2012.

4. Jacques Louis Lions. *Optimal Control of Systems Governed by Partial Differential Equations*. Grundlehren Der Mathematishen Wissenschaften 170. Springer, 1971.

5. A. Tarantola. *Inverse Problem Theory and Methods for Model Parameter Estimation*. Society for Industrial and Applied Mathematics, 2005.