Fictitious domain method for an inverse problem in volcanoes

An inverse problem applied to volcanology is studied. It consists in the determination of the pressure applied to a crack in order to ﬁt 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 ﬁctitious domain method, using a ﬁnite element discretization of XFEM type. The gradient of the cost function is derived from a Lagrangian. A numerical test is presented


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 1 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 Ω be a bounded open set in R d , d = 2, 3 with smooth boundary ∂Ω := Γ D ∪ Γ N where Γ D and Γ N are of nonzero measure and Γ D ∩ Γ N = ∅.We assume that Ω is occupied by an elastic solid and we denote by u the displacement field of the solid and the density of external forces by f ∈ L 2 (Ω).The Cauchy stress σ(u) and strain ε(u) are given by where (λ, µ) are the Lamé coefficients, I R d denotes the identity tensor, and Tr(•) represents the matrix trace.Consider a crack Γ C ⊂ Ω represented by a line (d = 2) or a surface (d = 3) parametrized by an injective mapping.Around the crack, Ω is split into Ω − and Ω + .The deformation field of the solid is supposed to satisfy the following elastostatic system: where n is the outward unit normal to its boundaries.Typically in such a situation, Γ N is the ground surface and free to move.Practically, the displacement field can be observed on Γ N , whereas the pressure p exerted on the crack is unknown most of the time.
Consider the following function defined on L 2 (Γ C ): where u d ∈ L 2 (Γ 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 α > 0 is a regularization parameter.The aim of this work is to study the following problem, of optimal control type: 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.

Domain decomposition : the direct solver
To solve the direct problem (1), we use a domain decomposition method.More precisely, following [1], the domain Ω 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 Γ C (e.g.Γ 0 in Figure 1).Therefore, instead of the crack Fig. 1 Splitting the volcanic cracked domain problem (1), we have to solve two Neumann-type boundary problems such that for each problem we impose a pressure on Γ C , which is more convenient from both theoretical and numerical points of view.More precisely we solve the following system: where u + = u |Ω + and u − = u |Ω , and v denotes the jump of v across Γ 0 .The two last conditions in (4) enforce the continuity of displacement and stress across Γ 0 .Notice that the boundary conditions on Γ 0 ensure the construction of a global displacement field in H 1 (Ω \ Γ C ) solving the original problem (1).
Let us define the following Hilbert spaces and their dual spaces V ± and W , endowed with their usual norms.Prescribing the continuity of displacement across the Γ 0 via a Lagrangian formulation, the mixed weak formulation of Problem (4) reads as follows: with a(u ± , v ± ) = bilinear, symmetric, coercive and 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 n = pn + = −pn − and the discretized form of system (5) has the linear algebraic formulation The system (6) can be solved by a Uzawa Conjugate gradient/domain decomposition method [1].The method can be classically stabilized and the convergence of the numerical scheme can be proved as h → 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.

The crack inverse Problem
First, we have the following result.
Proof The proof is classical: applying the same method as in [4], one easily shows that J is strictly convex and coercive on L 2 (Γ C ).
The objective function J being strictly convex, first order optimality conditions can be computed to implement an suitable optimization method (in our case the conjugate gradient).Let us introduce the adjoint system where u is a solution of system (1).Is is easy to prove that this adjoint system admits a unique solution φ ∈ H 1 (Ω).We have the following Proposition 2 Let p * be the solution of problem (3) and (u * , φ * ) be the associated solutions of (1) and (7).Then, the following optimality condition holds: 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 ∈ L 2 (Γ C ), compute (u, φ) which solve (1) and (7).Then, the Gâteaux derivative of J is given in L 2 (Γ C ) by For a given pressure p ∈ L 2 (Γ 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 Notice that the observed data u d can be interpolated on two sub-domains Ω ± to obtain u ± d corresponding to u ± .In view of (6), denoting R the reduction matrix R : X → U and the discrete cost function is defined as where X is the solution of (6), M N and M F are the mass matrices on Γ N and Γ C , respectively.This finite dimensional problem then boils down to finding the saddle point of the following Lagrangian Computing the KKT conditions for this problem allows to compute the gradient of J d : for a given vector p, let X be the solution of (6) and Φ be the solution of the adjoint problem Then, we have 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 computational cost.

Numerical experiments
Aiming at practical applications, we applied the technique to a realistic volcano, 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 (Figure 2, right).The convergence curves in Figure 3, highlight the efficiency of adapted optimization algorithm.The conjugate gradient minimization performs efficiently, even for fine meshes.

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)

Fig. 2
Fig. 2 Triangular surface mesh [2] representing the crack (left), amplitude of the displacement of a realistic volcano (right).