# HDG approximation of time-harmonic waves in elastic heterogeneous media with application to Full-Waveform Inversion

Hélène Barucq, Florian Faucher, Ha Pham

► To cite this version:

Hélène Barucq, Florian Faucher, Ha Pham. HDG approximation of time-harmonic waves in elastic heterogeneous media with application to Full-Waveform Inversion. Contemporary Challenges in Trefftz Methods, from Theory to Applications, CMO-BIRS workshop, May 2024, Oaxaca, Mexico. hal-04709800

HAL Id: hal-04709800

https://inria.hal.science/hal-04709800v1

Submitted on 25 Sep 2024

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
# HDG approximation of time-harmonic waves in elastic heterogeneous media with application to Full-Waveform Inversion

Hélène Barucq<fnref n="1" />, Florian Faucher<fnref n="1" /> and Ha Pham<fnref n="1" />.

## Contemporary Challenges in Trefftz Methods, from Theory to Applications

CMO-BIRS workshop, Oaxaca, 13–17 May 2024.

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

Inria

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

PROJECT
Exa-MA : Methods and
Algorithms for Exascale

<footnote n="1">¹Project-Team Makutu, Inria–TotalEnergies–UPPA, LMAP, UMR CNRS 5142.</footnote>
<figure><img src="image_4.png" /></figure>

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

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

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

# Table of Contents

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

## 1 Introduction

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

## 2 Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

* Illustration of gains

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

## 3 Inverse problem

* Iterative minimization algorithm

* Numerical experiment

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

## 4 Conclusion
<figure><img src="image_12.png" /></figure>

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

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

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

# Plan

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

## 1 Introduction
# Motivation

Seismic imaging for the **Reconstruction/Monitoring** of subsurface Earth properties, where elastic-wave problems remain computationally expensive.

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

* Accurate simulation of wave propagation in large-scale complex media,

* Efficient procedure for nonlinear reconstruction of properties.
<figure><img src="image_18.png" /></figure>

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

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

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

# Plan

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

## Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

* Illustration of gains
<figure><img src="image_23.png" /></figure>

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

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

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

# Plan

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

## Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

* Illustration of gains
# Elastic wave problem

New sensing devices such as fiber optic measure strain (ϵ). Therefore, we want to solve for (u, σ) to have maximum accuracy rather than replacing in terms of u only.

$$ \left\{ \begin{array}{l} -\omega^2 \mathbf{u} - \nabla \cdot \mathbf{\sigma} = \mathbf{f}, \\ \mathbf{\sigma} = \frac{1}{2} \mathbf{C} : \mathbf{\epsilon}, \quad \mathbf{\epsilon} = (\nabla \mathbf{u} + \nabla^T \mathbf{u}), \end{array} \right. \qquad (1) $$

The physical properties describing the medium are contained in the stiffness tensor C.

Time-harmonic wave problems:

+ Easily encode the attenuation with complex-valued parameters,

+ Direct solvers allow for multiple right-hand sides once the factorization is obtained,

- Memory cost of the matrix factorization.
<figure><img src="image_28.png" /></figure>

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

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

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

# HDG discretization

$$ \begin{cases} -\omega^2 \mathbf{u} - \nabla \cdot \boldsymbol{\sigma} = \mathbf{f}, \\ \boldsymbol{\sigma} = \frac{1}{2} \mathbf{C} : (\nabla \mathbf{u} + \nabla^T \mathbf{u}), \end{cases} $$

## Hybridizable Discontinuous Galerkin method

* **Static condensation** on first-order DG Formulations **without** increasing the number of unknown: The unknown of the global matrix is **only** the numerical trace $\hat{\mathbf{u}}$.

* **Handle complex geometry (topography) with p-adaptivity,**

* **Reduces the computational cost by removing inner dofs.**

<figure><img src="image_32.png" /><figcaption>Finite Element</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Discontinuous Galerkin</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>HDG</figcaption></figure>

<figure><img src="image_35.png" /></figure>
# HDG discretization

## Hybridizable Discontinuous Galerkin method for the discretization

* Static condensation for first-order problems **without** increasing the number of unknown: The unknown of the global matrix is **only** the numerical trace $\hat{u}$.

* Handle **complex geometry** (topography) with $p$-adaptivity,

* Reduces the computational cost by removing inner dofs.

<figure><img src="image_36.png" /><figcaption>Finite Element</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Discontinuous Galerkin</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>HDG</figcaption></figure>

HDG is more efficient with high-order polynomial (order $> 4$).
High-order $\Rightarrow$ Large cell $\Rightarrow$ variable $C$ within cells for resolution.
## HDG Variational formulation (1/2): stiffness tensor version

The local problem is written on each cell $K_e$ of the mesh, and the HDG problem is written in terms of $(u, \sigma, \hat{u})$. Using test-functions $(\psi, \phi)$,

$$ \int_{K_e} -\omega^2 \mathbf{u} \bar{\psi} - \int_{K_e} \nabla \cdot \sigma \bar{\psi} = \int_{K_e} \mathbf{f} \bar{\psi}, \quad (2a) $$

$$ \int_{K_e} \frac{1}{2} \mathbf{C} : \nabla \mathbf{u} \bar{\phi} + \int_{K_e} \frac{1}{2} \mathbf{C} : \nabla^T \mathbf{u} \bar{\phi} - \int_{K_e} \sigma \bar{\phi} = 0. \quad (2b) $$

To make appear the numerical trace $\hat{u}$, we need to integrate by parts leading to derivative of $\mathbf{C}$. When $\mathbf{C}$ is not constant per cell, one would need to provide its derivative...
<figure><img src="image_39.png" /></figure>

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

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

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

Variational formulation (2/2): compliance tensor version

$$ \text{Using } \mathbf{S} = \mathbf{C}^{-1}, \quad \begin{cases} -\omega^2 \mathbf{u} - \nabla \cdot \boldsymbol{\sigma} = \mathbf{f}, \\ \mathbf{S} : \boldsymbol{\sigma} = \frac{1}{2}(\nabla \mathbf{u} + \nabla^T \mathbf{u}), \end{cases} $$

On each cell $K_e$ of the mesh,

$$ \int_{K_e} -\omega^2 \mathbf{u} \bar{\psi} - \int_{K_e} \nabla \cdot \boldsymbol{\sigma} \bar{\psi} = \int_{K_e} \mathbf{f} \bar{\psi}, \quad (3a) $$

$$ \int_{K_e} \frac{1}{2} \nabla \mathbf{u} \bar{\phi} + \int_{K_e} \frac{1}{2} \nabla^T \mathbf{u} \bar{\phi} - \int_{K_e} \mathbf{S} : \boldsymbol{\sigma} \bar{\phi} = 0. \quad (3b) $$

⇒ using quadrature formula, $\mathbf{S}$ can easily vary within cell.

Note that under isotropy, we have an explicit formulation of $\mathbf{S}$ from Lamé parameters.
# HDG workflow in a nutshell

Unknowns are the discretized variables: $X^h = (\mathbf{u}^h, \sigma^h)$ and numerical trace $\hat{\mathbf{u}}^h$.

$$ \left\{ \begin{array}{ll} \text{local problem on each cell } K_e : & A_e X_e^h + C_e R_e \hat{\mathbf{u}}^h = \mathbf{F}, \\ \text{relation for numerical trace :} & \sum_e R_e^T (B_e X_e^h + L_e R_e \hat{\mathbf{u}}^h) = 0. \end{array} \right. $$

Reorder to write **global problem** in terms of numerical trace only

$$ \sum_e R_e^T (L_e - B_e A_e^{-1} C_e) R_e \hat{u}_h = - \sum_e R_e^T B_e A_e^{-1} F_e \quad \Leftrightarrow \quad A \hat{u}_h = B. $$

* Create local matrices on each cell (embarrassingly parallel) with *p*-adaptivity.

* Assemble global matrix $A$.

* Solve the large linear system $A\hat{\mathbf{u}}^h = B$ with MUMPS (solve multiple rhs at limited cost).

* Solve the local linear systems to obtain the volume solution $(\mathbf{u}^h, \sigma^h)$: (small matrices, embarrassingly parallel, < 2% of run time).
<figure><img src="image_43.png" /></figure>

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

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

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

Efficient stabilization for anisotropy

The efficiency of the HDG method relies onto the **stabilization matrix** $\tau$ that appears in the definition of the numerical trace:

$$ \widehat{\sigma}_h \nu = \sigma_h \nu - \tau(u_h - \widehat{u}_h). $$

A naive one (multiple of the identity) is **not robust**.
<figure><img src="image_47.png" /></figure>

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

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

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

Efficient stabilization for anisotropy

The efficiency of the HDG method relies onto the **stabilization matrix** $\tau$ that appears in the definition of the numerical trace:

$$ \widehat{\sigma}_h \nu = \sigma_h \nu - \tau(\boldsymbol{u}_h - \widehat{\boldsymbol{u}}_h). $$

A naive one (multiple of the identity) is **not robust**.

We instead construct a hybridized Godunov-upwind flux: Assuming three distinct wavespeeds in anisotropy ($c_p, c_{s1}, c_{s2}$),

$$ \tau = -i\omega \rho (c_p + c_{s1} + c_{s2}) \left(1 + \frac{\gamma}{\rho} (\Gamma + \rho \mathfrak{p})\right), $$

with $\Gamma$ the Kelvin–Christoffel matrix, $\gamma$ and $\mathfrak{p}$ auxiliary quantities depending on the wavespeed.

<figure><img src="image_51.png" /></figure>
<figure><img src="image_52.png" /></figure>

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

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

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

# Plan

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

## Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

* Illustration of gains
# Elastic-wave propagation with HDG

3D Model with topography, size 20×20×10 km³ with topography.

* Small cells required to accurately describe the topography, *p*-adaptivity

* Large cells elsewhere to benefit from HDG,

* models properties vary within the cell, here we use Lagrange basis per cell.

<figure><img src="image_57.png" /></figure>
# Elastic-wave propagation with HDG

3D Model with topography, size 20×20×10 km³ with topography.

* Small cells required to accurately describe the topography, *p*-adaptivity

* Large cells elsewhere to benefit from HDG,

* models properties vary within the cell, here we use Lagrange basis per cell.

mesh

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

polynomial orders.

<figure><img src="image_59.png" /></figure>
# Elastic-wave propagation with HDG

3D Model with topography, size 20×20×10 km³ with topography.

* Small cells required to accurately describe the topography, *p*-adaptivity

* Large cells elsewhere to benefit from HDG,

* models properties vary within the cell, here we use Lagrange basis per cell.

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

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

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

<figure><img src="image_63.png" /></figure>
# Matrix size reduction with HDG+MUMPS

3D elastic wave propagation, size $20 \times 20 \times 10$ km³ with topography.

Mesh using 120 000 cells

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

Polynomial order between 2 and 10

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

displacement field $u_z$ at 4 Hz

Global matrix size with frequency

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

elastic simulation 7 Hz with 1600 cores HDG+MUMPS:

* matrix size $N = 22.3 \times 10^6$,

* analysis time 2 min 30 sec,

* factorization time 34 min,

* factorization memory 3589 GiB,

* solve time (19 hrs) 1 min 40 s.
<figure><img src="image_67.png" /></figure>

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

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

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

Anisotropic modeling with highly-varying model parameters

2D anisotropic simulation with solar-like model background:

* Density is exponentially decreasing near surface,

* Velocity is collapsing near surface,

* Anisotropy encoded with Thomsen's parameters $\epsilon = 0.25$, $\delta = 0.15$, $\theta = 45^\circ$.

<figure><img src="image_71.png" /><figcaption>2D disk mesh used</figcaption></figure>

<figure><img src="image_72.png" /></figure>
<figure><img src="image_73.png" /></figure>

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

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

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

# Anisotropic modeling with highly-varying model parameters

2D anisotropic simulation with solar-like model background:

* Density is exponentially decreasing near surface,

* Velocity is collapsing near surface,

* Anisotropy encoded with Thomsen's parameters $\epsilon = 0.25$, $\delta = 0.15$, $\theta = 45^\circ$.

<figure><img src="image_77.png" /><figcaption>2D disk mesh used</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Radial density</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>(a) $\sqrt{\rho} \text{ Re}(u_x)$.</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>(c) $\text{Re}(\sigma_{xx}) / \sqrt{\rho}$. (d) $\text{Re}(\sigma_{zz}) / \sqrt{\rho}$.</figcaption></figure>

<caption>(b) $\sqrt{\rho} \text{ Re}(u_z)$.</caption>

<caption>(e) $\text{Re}(\sigma_{xz}) / \sqrt{\rho}$.</caption>
<figure><img src="image_81.png" /></figure>

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

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

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

# Plan

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

## 3 Inverse problem

* Iterative minimization algorithm

* Numerical experiment
<figure><img src="image_86.png" /></figure>

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

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

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

# Plan

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

## 3 Inverse problem

* Iterative minimization algorithm

* Numerical experiment
Quantitative reconstruction algorithm: FWI

Quantitative reconstruction of properties $\mathbf{m} = (\lambda, \mu, \rho)$ solving iterative minimization problem.

$\min_{\mathbf{m}} \mathcal{J}(\mathbf{m})$ with $\mathcal{J}(\mathbf{m}) = \text{dist}(\mathcal{F}(\mathbf{m}), \mathbf{d})$, $\mathcal{F}$: simulations, $\mathbf{d}$: observations.

<figure><img src="image_91.png" /></figure>
<figure><img src="image_92.png" /></figure>

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

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

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

Quantitative reconstruction algorithm: FWI

minimize $\mathcal{J}(\mathbf{m}) = \text{dist}(\mathcal{F}(\mathbf{m}), \mathbf{d})$

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

1. Repeated use of forward wave propagation solver,

$$ \mathcal{A}_e \mathbf{X}_e^h + \mathcal{C}_e \mathcal{R}_e \hat{\mathbf{u}}^h = \mathbf{F}, \quad (4a) $$

$$ \sum_e \mathcal{R}_e^T (\mathcal{B}_e \mathbf{X}_e^h + \mathcal{L}_e \mathcal{R}_e \hat{\mathbf{u}}^h) = 0. \quad (4b) $$

2. Adjoint-state method for gradient adapted to HDG Lagrangian written from $\mathcal{J}$ subject to (4a), (4b).

We prove that for HDG, the gradient is still computed from the **adjoint of the direct problem**, fundamental to avoid refactorization of the global matrix.
<figure><img src="image_97.png" /></figure>

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

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

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

Quantitative reconstruction algorithm: FWI

minimize $\mathcal{J}(\mathbf{m}) = \operatorname{dist}(\mathcal{F}(\mathbf{m}), \mathbf{d})$

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

1. Repeated use of forward wave propagation solver,

2. Adjoint-state method for gradient adapted to HDG,

3. To alleviate mesh limitations, the model parameters are represented with Lagrange basis functions.

4. Inversion is carried out with respect to the weight of the Lagrange basis functions.
<figure><img src="image_102.png" /></figure>

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

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

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

# Plan

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

## 3 Inverse problem

* Iterative minimization algorithm

* Numerical experiment
## 2D experiments: Marmousi II setup

We consider the elastic isotropic Marmousi II experiment, of size $17 \times 3.5 \text{ km}^2$. Free-surface boundary condition on top and absorbing boundary conditions elsewhere.

<figure><img src="image_107.png" /><figcaption>True P-wave speed</figcaption></figure>

<figure><img src="image_108.png" /><figcaption>True S-wave speed</figcaption></figure>

<figure><img src="image_109.png" /><figcaption>True density</figcaption></figure>
<figure><img src="image_110.png" /></figure>

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

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

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

## 2D experiments: Marmousi II setup

We consider the elastic isotropic Marmousi II experiment, of size 17×3.5 km².
Free-surface boundary condition on top and absorbing boundary conditions elsewhere.

<figure><img src="image_114.png" /><figcaption>True P-wave speed</figcaption></figure>

<figure><img src="image_115.png" /><figcaption>True S-wave speed</figcaption></figure>

<figure><img src="image_116.png" /><figcaption>True density</figcaption></figure>

Initial models are 1D-background variation.

<figure><img src="image_117.png" /><figcaption>Initial P-wave speed</figcaption></figure>

<figure><img src="image_118.png" /><figcaption>Initial S-wave speed</figcaption></figure>

<figure><img src="image_119.png" /><figcaption>Initial density</figcaption></figure>
<figure><img src="image_120.png" /></figure>

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

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

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

## 2D experiments: Marmousi II reconstructions

* Acquisition made up of 169 sources and 849 receivers near surface.

* Reconstructions using 13 frequencies from 2 to 8 Hz, 25 iterations per frequency.

* Density is not inverted and remains in its initial value.

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

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

<figure><img src="image_126.png" /></figure>
<figure><img src="image_127.png" /></figure>

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

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

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

## 2D experiments: Marmousi II reconstructions

* Acquisition made up of 169 sources and 849 receivers near surface.

* Reconstructions using 13 frequencies from 2 to 8 Hz, 25 iterations per frequency.

* Density is not inverted and remains in its initial value.

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

* Reconstruction using representation in order 1 Lagrange basis per cell (~20 × 10³ cells)

<figure><img src="image_132.png" /></figure>
<figure><img src="image_133.png" /></figure>

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

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

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

## 2D experiments: Marmousi II reconstructions

* Acquisition made up of 169 sources and 849 receivers near surface.

* Reconstructions using 13 frequencies from 2 to 8 Hz, 25 iterations per frequency.

* Density is not inverted and remains in its initial value.

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

* Computational time (18mpi × 2omp): $2 \cdot 10^4$ cells: 3h; $5 \cdot 10^4$ cells: 4h15min, i.e. -30%.

<figure><img src="image_138.png" /><figcaption>S-wave speed reconstruction using mesh with 20000 cells and Lagrange basis representation,</figcaption></figure>

<figure><img src="image_139.png" /><figcaption>S-wave speed reconstruction using mesh with 50000 cells piecewise-constant representation.</figcaption></figure>
<figure><img src="image_140.png" /></figure>

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

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

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

# Plan

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

## Conclusion
<figure><img src="image_145.png" /></figure>

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

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

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

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

# Conclusion

We propose an efficient framework for elastic inverse problem using FWI method.

## Hybridizable discontinuous Galerkin discretization

* Handle complex geometry and p-adaptivity, efficient for parallelization,

* Reduction of the size of the linear system with **high-order polynomials**,

* Model representation must allow variation within cells,

## Inversion

* Adjoint-state method for HDG re-uses forward global matrix factorization,

* Reconstruction with model representation to improve resolution.

* 3D running on toy problems, scalability analysis for regional-scale imaging is on going

* Open-source software hawen: https://ffaucher.gitlab.io/hawen-website/
<figure><img src="image_150.png" /></figure>

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

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

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

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

# Conclusion

We propose an efficient framework for elastic inverse problem using FWI method.

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

## Hybridizable discontinuous Galerkin discretization

* Handle complex geometry and p-adaptivity, efficient for parallelization,

* Reduction of the size of the linear system with **high-order polynomials**,

* Model representation must allow variation within cells,

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

## Inversion

* Adjoint-state method for HDG re-uses forward global matrix factorization,

* Reconstruction with model representation to improve resolution.

* 3D running on toy problems, scalability analysis for regional-scale imaging is on going

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

Open-source software hawen: https://ffaucher.gitlab.io/hawen-website/

THANK YOU