# Understanding the invisible through numerical simulation of wave propagation

Hélène Barucq

► To cite this version:

Hélène Barucq. Understanding the invisible through numerical simulation of wave propagation. SMAI 2023 - 11e Biennale française de Mathématiques Appliquées et Industrielles, Université des Antilles, May 2023, Le Gosier (Guadeloupe), France. http://dx.doi.org/10.44013/antilles.2023.11

HAL Id: hal-04401321

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

Submitted on 17 Jan 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
Understanding the invisible through numerical simulation of wave propagation

Hélène Barucq <fnref n="1" />.

SMAI 2023 – Le Gosier, Guadeloupe – May 2023

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

MAKUTU

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

<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>

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

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

# Table of Contents

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

## 1 Introduction

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

## 2 Energy Geophysics

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

## 3 Illustration: Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

* Illustration of gains

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

## 4 Inverse problem

* Iterative minimization algorithm

* Numerical experiment

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

## 5 Helioseismology/asteroseismology

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

## 6 Illustration: Green's function computations for time-distance helioseismology

* Modal Green's kernel

* Numerical experiments

* Power spectrum, comparison with measured data

* The vector wave problem

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

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

<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>

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

# Plan

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

## 1 Introduction
<figure><img src="image_24.png" /></figure>

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

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

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

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

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

Broad scientific context

Describe a place with exactness from more or less numerous and precise memories, or guessing the content and internal structures of an object after having observed it only partially, without ever touching it because it is inaccessible or very fragile?

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

Wave propagation is helping

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

Waves are very sensitive to any change in the propagation medium.

<figure><img src="image_32.png" /><figcaption>(a) Seismogram.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(b) Reservoir model.</figcaption></figure>
<figure><img src="image_34.png" /></figure>

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

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

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

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

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

# Inverse problems

Basically, inverse wave problems are composed of:

* Emitting sources that will propagate through the medium and recording the reflected waves on a set of receivers; acquisition/forward problem

* From the acquisitions, find the propagation medium; inverse problem/backward problem

Example: Seismic imaging for the **Reconstruction/Monitoring** of subsurface Earth properties

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

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

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

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

<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>

Numerical data match real data

Compare data with simulations: computationally and possibly timely expensive

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

Comparison of observation and numerical result
Plan
# Access to energy resources

* Hydrocarbons

* Geothermal

* Hydrogen

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

Source: The Emissions Gap Report 2017. United Nations Environment Programme (UNEP)
# Facilitate safe carbon capture, utilization and storage - CCUS

* CO2 emissions: leading cause of climate change;

* Geological storage of CO2: important tool for the stabilization of atmospheric greenhouse gas concentrations;

* CO2 is injected into underground geological formations; Safe, permanent, and effective.

## How do we ensure the safety and sustainability of storage?

That's where **seismic monitoring** comes into play.

<figure><img src="image_49.png" /></figure>
# Research routine: propagation domains

* 3D large domains

* topography, heterogeneity

* a network of sources whose number is of the order of several thousands

<figure><img src="image_50.png" /><figcaption>Parameterization matters</figcaption></figure>
# Research routine: wave equations

* Acoustic wave equation,

* Elastic wave equation,

* Electromagnetic wave equation,

* Couplings

<figure><img src="image_51.png" /><figcaption>Waves in solid media</figcaption></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>

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

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

Equations 2/2

Acoustic wave equation, time and frequency domains:

$$ \begin{cases} \rho \partial_t \mathbf{v}(\mathbf{x}, t) = -\nabla p(\mathbf{x}, t), \\ \frac{1}{c^2 \rho} \partial_t p(\mathbf{x}, t) + \nabla \cdot \mathbf{v}(\mathbf{x}, t) = 0. \end{cases} $$

$$ \begin{cases} -i \omega \rho \mathbf{v}(\mathbf{x}, \omega) = -\nabla p(\mathbf{x}, \omega), \\ -\frac{1}{c^2 \rho} i \omega \rho \mathbf{v}(\mathbf{x}, \omega) + \nabla \cdot \mathbf{v}(\mathbf{x}, \omega) = 0. \end{cases} $$

Elasto-dynamic equations, time and frequency domains:

$$ \begin{cases} \rho \partial_t \mathbf{v}(\mathbf{x}, t) = \nabla \cdot \underline{\underline{\sigma}}(\mathbf{x}, t), \\ \partial_t \underline{\underline{\sigma}}(\mathbf{x}, t) = \underline{\underline{C}}(\mathbf{x})(\underline{\underline{\epsilon}}(\mathbf{v})). \end{cases} $$

$$ \begin{cases} -i\omega\rho\mathbf{v}(\mathbf{x}, \omega) = \nabla \cdot \underline{\underline{\sigma}}(\mathbf{x}, \omega), \\ -i\omega\underline{\underline{\sigma}}(\mathbf{x}, \omega) = \underline{\underline{C}}(\mathbf{x})(\underline{\underline{\epsilon}}(\mathbf{v})). \end{cases} $$
<figure><img src="image_58.png" /></figure>

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

<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>

Equations 2/2

Each approach has pros and cons

Time-domain formulation

* Using data is straightforward

* Matrix-free implementation relieves memory,

* frequency-dependent parameters,

* adjoint of the discrete problem can differ from discrete adjoint, additional developments are required,

* multi-sources

Time-harmonic formulation

* Easily handle attenuation,

* multi-sources with direct solver,

* reuse matrix factorization for adjoint-problem in inversion,

* memory cost for matrix factorization with direct solver,
<figure><img src="image_64.png" /></figure>

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

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

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

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

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

# Modelling and simulation challenges 1/2

## A multi-physics problem and...

* ✓ CO2 sequestration and monitoring are more and more recognized as a key element in the path towards energy decarbonization

* ✓ CO2 sequestration and monitoring is an excellent example gathering the leading-edge technology of many different domain to help in the prediction of the plume evolution:

* ✓ Flow and geo-mechanical simulation,

* ✓ Gravimetry,

* ✓ Seismic modelling and inverse problem, monitoring acquisition technology,

* ✓ In situ data visualization and analysis,

* ✓ Machine learning

* ✓ ...

Source: Britto

CO₂ injection site

1. Well integrity/injectivity

2. Pressure/Stress change
Fault Activation

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

3. CO₂ transport & trapping

4. Seal integrity

5. Surface deformation

6. Seismicity
<figure><img src="image_71.png" /></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>

# Modelling and simulation challenges 2/2

## HPC Issue

### Limitations:

* ✓ Multiphysics: geomechanic+Flow+Seismic

* ✓ Large Scale: 98% storage in Aquifer

* ✓ Long Term Simulation: post injection matters

### Solutions:

* ✓ Fast and scalable algorithms

* ✓ Perennial solutions: portability

## Target: Exascale Computers

Source: Blice

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

CO₂ injection site

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

1. Well integrity/injectivity

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

2. Pressure/Stress change
Fault Activation

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

3. CO₂ transport & trapping

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

4. Seal integrity

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

5. Surface deformation

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

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

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

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

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

## Finite Differences: Implementation

Easy, low cost, inaccuracy for the topography

* Finite Elements: implementation possibly tricky, expensive, accurate for the topography

* Boundary integral equations: not efficient in highly heterogeneous media

* Semi-Analytical: lack of flexibility, geometrical effects, anisotropy neglected,

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

<figure><img src="image_89.png" /></figure>
# In Makutu team 2/2

* Spectral element methods (SEM)

* Discontinuous Galerkin methods (DG)

* Non polynomial basis functions in Trefftz framework

* Explicit schemes in time

* High-order

* hp-adaptivity with DG

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

Multiphysics
Simulator on HPC

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

Hawen

https://ffaucher.gitlab.io/hawen-website/
<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>

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

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

# Plan

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

## Illustration: Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

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

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

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

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

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

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

# Elastic wave propagation with HDG

In collaboration with Florian Faucher and Ha Pham

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

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

3 Illustration: Elastic wave propagation with HDG

* Mathematical formulation and HDG algorithm

* Illustration of gains
# 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_107.png" /><figcaption>Finite Element</figcaption></figure>

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

<figure><img src="image_109.png" /><figcaption>HDG</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>

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

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

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

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

# HDG workflow in a nutshell

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

$$ 
\left\{ 
\begin{array}{ll} 
\text{local problem on each cell } \mathcal{K}_e : & \mathcal{A}_e X_e^h + \mathcal{C}_e \mathcal{R}_e \hat{\mathbf{u}}^h = \mathbf{F}, \\ 
\text{relation for numerical trace :} & \sum_e \mathcal{R}_e^T (\mathcal{B}_e X_e^h + \mathcal{L}_e \mathcal{R}_e \hat{\mathbf{u}}^h) = 0. 
\end{array} 
\right. 
$$

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

$$ \sum_{e} \mathcal{R}_{e}^{T} (\mathcal{L}_{e} - \mathcal{B}_{e} \mathcal{A}_{e}^{-1} \mathcal{C}_{e}) \mathcal{R}_{e} \hat{\mathbf{u}}_{h} = - \sum_{e} \mathcal{R}_{e}^{T} \mathcal{B}_{e} \mathcal{A}_{e}^{-1} \mathcal{F}_{e} \quad \Leftrightarrow \quad \mathcal{A} \hat{\mathbf{u}}_{h} = \mathcal{B}. $$

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

* Assemble global matrix $\mathcal{A}$.

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

* Solve the local linear systems to obtain the volume solution $(\mathbf{u}^h, \boldsymbol{\sigma}^h)$: (small matrices, embarrassingly parallel, < 2% of run time).
# Plan

3 Illustration: 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_118.png" /></figure>
<figure><img src="image_119.png" /></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>

<figure><img src="image_124.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_125.png" /></figure>

mesh

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

polynomial orders.
# 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_127.png" /></figure>

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

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

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

<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>

# Matrix size reduction with HDG+MUMPS

3D elastic wave propagation, size 20×20×10 km³ with topography.

Mesh using 120 000 cells

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

Polynomial order between 2 and 10

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

displacement field $u_z$ at 4 Hz

Global matrix size with frequency

<figure><img src="image_139.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_140.png" /></figure>

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

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

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

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

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

# Plan

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

## 4 Inverse problem

* Iterative minimization algorithm

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

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

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

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

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

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

# Plan

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

## 4 Inverse problem

* Iterative minimization algorithm

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

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

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

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

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

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

Quantitative reconstruction algorithm: FWI

Quantitative reconstruction of properties **m** = (λ, μ, ρ) solving **iterative minimization problem**.

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

<figure><img src="image_160.png" /></figure>
<figure><img src="image_161.png" /></figure>

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

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

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

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

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

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

Quantitative reconstruction algorithm: FWI

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

<figure><img src="image_168.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_169.png" /></figure>

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

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

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

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

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

Quantitative reconstruction algorithm: FWI

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

<figure><img src="image_175.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.
# Plan

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

## 4 Inverse problem

* Iterative minimization algorithm

* Numerical experiment
## 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_177.png" /><figcaption>True P-wave speed</figcaption></figure>

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

<figure><img src="image_179.png" /><figcaption>True density</figcaption></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_180.png" /></figure>

True P-wave speed

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

True S-wave speed

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

True density

Initial models are 1D-background variation.

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

Initial P-wave speed

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

Initial S-wave speed

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

Initial density
## 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_186.png" /></figure>

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

<figure><img src="image_188.png" /></figure>
<figure><img src="image_189.png" /></figure>

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

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

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

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

<figure><img src="image_194.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_195.png" /></figure>

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

<figure><img src="image_196.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_197.png" /></figure>

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

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

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

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

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

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

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

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

# Plan

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

## Helioseismology/asteroseismology
# Helioseismic studies

* Measuring solar/star oscillations,

* Processing and averaging the observations to extract the seismic data

* Interpreting the seismic data using forward and inverse methods to estimate solar internal properties.

<figure><img src="image_207.png" /><figcaption>Solar interior</figcaption></figure>
<figure><img src="image_208.png" /></figure>

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

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

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

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

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

Detecting active regions on the far side of the Sun

► Of great importance for space-weather forecasts

► Large active regions that emerge on the Sun's far side will rotate into Earth's view several days later; these may trigger coronal mass ejections, which can damage satellites and spacecraft and endanger astronauts

► It is known that far-side imaging can significantly improve models of the solar wind, which plays an important role in space-weather forecasts.

<figure><img src="image_214.png" /></figure>
<figure><img src="image_215.png" /></figure>

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

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

# Numerical simulations of acoustic waves

Acoustic waves propagate horizontally and are trapped in the vertical direction; they connect the Sun's near and far sides. As acoustic waves travel faster in magnetized regions, they can inform us about the presence of active regions along their paths of propagation.

Two helioseismic techniques:

* **Helioseismic holography:** multiply the forward and backward propagated wavefields, and subtract a reference measurement for the quiet Sun.

* **Time-distance helioseismology:** compute the cross-covariance of the wavefield

**Develop models and computational tools** that are required to recover information about Sun and eventually stellar activity from observations of oscillations.

C. Lindsey and D.C. Braun, *Helioseismic holography*. *The Astrophysical Journal*, 1997

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

L. Gizon, H. Barucq, M. Duruffé, C.S. Hanson, M. Leguèbe, A.C. Birch, J. Chabassier, D. Fournier, T. Hohage and E. Papini, *Computational helioseismology in the frequency domain: acoustic waves in axisymmetric solar models with flows*. *Astronomy & Astrophysics*, 2017

<figure><img src="image_219.png" /></figure>
<figure><img src="image_220.png" /></figure>

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

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

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

# Equations

* Galbrun's equation describes adiabatic wave motion subject to source $F$ with frequency $\omega$, on top of a time-invariant fluid background. Ignoring flows and rotation and with Cowling's approximation,

$$ -\rho_0(\omega^2+2i\omega\Gamma)\xi - \nabla[\gamma\rho_0\nabla\cdot\xi] + (\nabla\rho_0)(\nabla\cdot\xi) - \nabla(\xi\cdot\nabla\rho_0) + (\xi\cdot\nabla)\nabla\rho_0 + \rho_0(\xi\cdot\nabla)\nabla\phi_0 = F, $$
with displacement vector $\xi$ describing the solar oscillation. The density is $\rho_0$, pressure $p_0$, adiabatic index $\gamma$, gravitational potential $\phi_0$ and attenuation $\Gamma$.

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

<figure><img src="image_225.png" /></figure>
<figure><img src="image_226.png" /></figure>

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

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

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

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

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

# Equations

► **Galbrun's equation** describes adiabatic wave motion subject to source $F$ with frequency $\omega$, on top of a time-invariant fluid background. Ignoring flows and rotation and with Cowling's approximation,

$$-\rho_0(\omega^2+2i\omega\Gamma)\xi - \nabla[\gamma\rho_0\nabla\cdot\xi] + (\nabla\rho_0)(\nabla\cdot\xi) - \nabla(\xi\cdot\nabla\rho_0) + (\xi\cdot\nabla)\nabla\rho_0 + \rho_0(\xi\cdot\nabla)\nabla\phi_0 = F,$$ with displacement vector $\xi$ describing the solar oscillation. The density is $\rho_0$, pressure $p_0$, adiabatic index $\gamma$, gravitational potential $\phi_0$ and attenuation $\Gamma$.

► Further approximation ignoring the gravity, we obtain a **scalar-wave equation** for variable $u := \rho c^2 \nabla \cdot \xi$ such that,

$$-\nabla \cdot \left( \frac{1}{\rho} \nabla u \right) - \frac{\omega^2 + 2i\omega\Gamma}{\rho c^2} u = f, \quad \text{with sound-speed c.}$$

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. *SIAM J. on Appl. Math.*, 2020.

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Outgoing modal solutions for Galbrun's equation in helioseismology. *J. of Differential Equ.*, 2021.
<figure><img src="image_234.png" /></figure>

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

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

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

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

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

Spherical background solar model

Model S+Atmo: exponentially decreasing density in the solar atmosphere

$$ \rho(r) := \begin{cases} \rho_S(r), & r \le r_a \\ \rho_S(r_a) e^{-\alpha_\infty (r-r_a)}, & r > r_a \end{cases}, $$

$$ c(r) := \begin{cases} c_S(r), & r \le r_a \\ c_S(r_a), & r > r_a \end{cases} $$

with $\rho_S > \rho_0 > 0$ ; $c_S > c_0 > 0$, $R_\odot$ = the solar radius ($\sim 696 \times 10^3$km)
<figure><img src="image_240.png" /></figure>

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

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

## Spherical background solar model

### Model S+Atmo: exponentially decreasing density in the solar atmosphere

$$ \rho(r) := \begin{cases} \rho_S(r), & r \le r_a \\ \rho_S(r_a) e^{-\alpha_\infty (r-r_a)}, & r > r_a \end{cases}, $$

$$ c(r) := \begin{cases} c_S(r), & r \le r_a; \\ c_S(r_a), & r > r_a \end{cases} $$

Comparison of the S+Atmo (solid) and VAL-C (dashed) atmospheric models.

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

<table><tr><td>Height above surface (Mm) R<sub>h</sub></td><td>0.5</td><td>0.556</td><td>2.543</td><td>4</td></tr><tr><td>scaled radius = (R<sub>☉</sub> + R<sub>h</sub>)/R<sub>☉</sub></td><td>1.0007</td><td>1.0008</td><td>1.0037</td><td>1.0058</td></tr></table>
Scalar-wave problem formulation, spherical symmetry

$$ - \nabla \cdot \left( \frac{1}{\rho} \nabla u \right) - \frac{\omega^2 + 2i\omega\Gamma}{\rho c^2} u = f. $$

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. SIAM J. on Appl. Math., 2020.

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

H. Barucq, F. Faucher, H. Pham, Outgoing solutions and radiation boundary conditions for the ideal atmospheric scalar wave equation in helioseismology. ESAIM, 2020.
<figure><img src="image_246.png" /></figure>

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

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

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

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

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

Scalar-wave problem formulation, spherical symmetry

$$-\nabla \cdot \left( \frac{1}{\rho} \nabla u \right) - \frac{\omega^2 + 2i\omega\Gamma}{\rho c^2} u = f.$$

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

We use the Liouville transform to rewrite the problem in Schrödinger form,

With $u = \rho^{-1/2}u$ and $\alpha = -\rho'/\rho$, we obtain

$$\left( -\Delta - \frac{\omega^2 + 2i\omega\Gamma}{c^2} + \frac{\alpha^2}{4} + \frac{\alpha'}{2} + \frac{\alpha}{r} \right) u = f.$$

In the atmosphere, $\rho$ is exponentially decreasing $\Rightarrow \alpha$ is constant.

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. *SIAM J. on Appl. Math.*, 2020.

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

H. Barucq, F. Faucher, H. Pham, Outgoing solutions and radiation boundary conditions for the ideal atmospheric scalar wave equation in helioseismology. *ESAIM, 2020*.
<figure><img src="image_255.png" /></figure>

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

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

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

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

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

Scalar-wave problem formulation, spherical symmetry

$$-\nabla \cdot \left( \frac{1}{\rho} \nabla u \right) - \frac{\omega^2 + 2i\omega\Gamma}{\rho c^2} u = f.$$

* We use the Liouville transform to rewrite the problem in Schrödinger form,

With $u = \rho^{-1/2}u$ and $\alpha = -\rho'/\rho$, we obtain

$$\left(-\Delta - \frac{\omega^2 + 2i\omega\Gamma}{c^2} + \frac{\alpha^2}{4} + \frac{\alpha'}{2} + \frac{\alpha}{r}\right) u = f,$$

In the atmosphere, $\rho$ is exponentially decreasing $\Rightarrow \alpha$ is constant.

* Allows us to build the exact Dirichlet-to-Neumann map and new classes of Radiation Boundary conditions for outgoing solutions.

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. SIAM J. on Appl. Math., 2020.

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

H. Barucq, F. Faucher, H. Pham, Outgoing solutions and radiation boundary conditions for the ideal atmospheric scalar wave equation in helioseismology. ESAIM, 2020.
# Plan

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

## Illustration: Green's function computations for time-distance helioseismology

* Modal Green's kernel

* Numerical experiments

* Power spectrum, comparison with measured data

* The vector wave problem
<figure><img src="image_264.png" /></figure>

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

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

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

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

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

Numerical method for computing Green's functions

In collaboration with Florian Faucher, Ha Pham, and Damien Fournier, Laurent Gizon from Max Planck Institute for Sun (MPS, Göttingen, Germany)

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

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

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

<figure><img src="image_273.png" /></figure>
# Plan

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

## Illustration: Green's function computations for time-distance helioseismology

* Modal Green's kernel

* Numerical experiments

* Power spectrum, comparison with measured data

* The vector wave problem
<figure><img src="image_275.png" /></figure>

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

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

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

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

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

Modal Green's kernel

Decomposition on harmonic mode $\ell$, the modal Green's kernel $G_\ell$ solves,

$$ \mathcal{L} G_\ell(r, s) = \delta(r - s), \quad \text{with} \quad \mathcal{L} := \left( -\partial_r^2 - \frac{\omega^2 + 2i\omega\Gamma}{c^2} + \frac{\alpha^2}{4} + \frac{\alpha'}{2} + \frac{\alpha}{r} + \frac{\ell(\ell+1)}{r^2} \right). $$

with Neumann boundary condition at origin:

$$ \partial_n G_\ell(r = 0) = 0, $$

and radiation condition at $r = r_{\max}$:

$$ (\partial_n G_\ell - \mathcal{Z} G_\ell)_{r=r_{\max}} = 0. $$

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

H. Barucq, F. Faucher, H. Pham, Outgoing solutions and radiation boundary conditions for the ideal atmospheric scalar wave equation in helioseismology.
ESAIM, 2020.
<figure><img src="image_282.png" /></figure>

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

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

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

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

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

Modal Green's kernel

Decomposition on harmonic mode $\ell$, the modal Green's kernel $G_\ell$ solves,

$$ \mathcal{L} G_\ell(r, s) = \delta(r-s), \quad \text{with} \quad \mathcal{L} := \left( -\partial_r^2 - \frac{\omega^2 + 2i\omega\Gamma}{c^2} + \frac{\alpha^2}{4} + \frac{\alpha'}{2} + \frac{\alpha}{r} + \frac{\ell(\ell+1)}{r^2} \right). $$

with Neumann boundary condition at origin:
and radiation condition at $r = r_{\max}$:

$$ \begin{aligned} \partial_n G_\ell(r=0) &= 0, \\ (\partial_n G_\ell - \mathcal{Z} G_\ell)_{r=r_{\max}} &= 0. \end{aligned} $$

Computation of the modal Green's kernel:

* **Approach 1 (naive):** for each source position $s$, solving equation $\mathcal{L} G_\ell(r, s) = \delta(r-s)$ gives $G_\ell(\cdot, s)$.

* **Approach 2:** Use an assembling formula, the entire kernel is obtained from the solutions of only two boundary-value problems.

<figure><img src="image_288.png" /></figure>
<figure><img src="image_289.png" /></figure>

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

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

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

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

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

Computation of the Green's kernel: Approach 2

► **Approach 1 (naive):** for each source position $s$, solve equation $\mathcal{L} G_\ell(r, s) = \delta(r - s)$.

► **Approach 2:** Use an assembling formula, the entire kernel is obtained from the solutions of only two boundary-value problems.

Find $\psi_1$ that solves
$$
\begin{cases}
\mathcal{L} \psi_1 = 0, & \text{on } (0, r_{\max}), \\
(\partial_n \psi_1)_{r=0} = 0, & (\psi_1)_{r=r_b} = 1.
\end{cases}
$$

Find $\psi_2$ that solves
$$
\begin{cases}
\mathcal{L} \psi_2 = 0, & \text{on } (r_a, r_{\max}), \\
(\psi_2)_{r=r_a} = 1, & (\partial_n \psi_2 - \mathcal{Z} \psi_2)_{r=r_{\max}} = 0.
\end{cases}
$$

<figure><img src="image_295.png" /></figure>
Computation of the Green's kernel: Approach 2

► **Approach 1 (naive):** for each source position $s$, solve equation $\mathcal{L} G_\ell(r, s) = \delta(r - s)$.

► **Approach 2:** Use an assembling formula, the entire kernel is obtained from the solutions of only two boundary-value problems.

* Find $\psi_1$ that solves
$$ \begin{cases} \mathcal{L} \psi_1 = 0, & \text{on } (0, r_{\max}), \\ (\partial_n \psi_1)_{r=0} = 0, & (\psi_1)_{r=r_b} = 1. \end{cases} $$

* Find $\psi_2$ that solves
$$ \begin{cases} \mathcal{L} \psi_2 = 0, & \text{on } (r_a, r_{\max}), \\ (\psi_2)_{r=r_a} = 1, & (\partial_n \psi_2 - Z \psi_2)_{r=r_{\max}} = 0. \end{cases} $$

* Assemble the Green's kernel, with $H$ the Heaviside and $W$ the Wronskian
  $W(s) := W\{\psi_1(s), \psi_2(s)\}$,

$$ G_{\ell}(r, s) = \frac{-H(s-r)\psi_{1}(r)\psi_{2}(s) - H(r-s)\psi_{1}(s)\psi_{2}(r)}{W(s)}, \quad \forall (r, s) \in (r_a, r_{\max}), $$

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. SIAM J. on Appl. Math., 2020.
<figure><img src="image_297.png" /></figure>

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

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

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

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

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

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

Computation of the Green's kernel: Approach 2

► **Approach 1 (naive):** for each source position $s$, solve equation $\mathcal{L} G_\ell(r, s) = \delta(r-s)$.

► **Approach 2:** Use an assembling formula, the entire kernel is obtained from the solutions of only two boundary-value problems.

**Approach 1:** the solution of one problem for a Dirac in $s_0$ only gives $G_\ell(r, s = s_0)$ and $G_\ell(r = s_0, s)$.

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

**Approach 2:** from the solutions of two boundary value problems, $G_\ell(r, s)$ is obtained for any position between $r_a$ and $r_b$.

It also avoids the singularity of the Dirac source.

<figure><img src="image_305.png" /></figure>
# Plan

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

## Illustration: Green's function computations for time-distance helioseismology

* Modal Green's kernel

* Numerical experiments

* Power spectrum, comparison with measured data

* The vector wave problem
## Solar Green's kernels

* The background sound-speed and density are given by model S+Atmo or model Val-C.

<figure><img src="image_307.png" /></figure>
<figure><img src="image_308.png" /></figure>

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

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

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

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

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

# Solar Green's kernels

* The background sound-speed and density are given by model S+Atmo or model Val-C.

* The Green's kernel is computed with Approach 2 which only need the solutions of two problems.

Imaginary part of the Solar modal Green's functions at 7 mHz for mode $ℓ = 100$ using different background.

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

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. SIAM J. on Appl. Math., 2020

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

F. Faucher

hawen: time-harmonic wave modeling and inversion using hybridizable discontinuous Galerkin discretization

Journal of Open Source Software, 6 (57), 2021
# Solar Green's kernels

* The background sound-speed and density are given by model S+Atmo or model Val-C.

* The Green's kernel is computed with Approach 2 which only needs the solutions of two problems.

Approach 2 is more **efficient** as it only needs two boundary value problems to assemble the entire kernel, and more **accurate** as it avoids the Dirac singularity.

* To obtain $G_\ell(s, s)$ with Approach 1, we need 4000 simulations,

* Approach 1 needs a refined discretization mesh to capture correctly the singularity.

<figure><img src="image_317.png" /></figure>
# Plan

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

## Illustration: Green's function computations for time-distance helioseismology

* Modal Green's kernel

* Numerical experiments

* Power spectrum, comparison with measured data

* The vector wave problem
<figure><img src="image_319.png" /></figure>

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

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

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

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

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

# Power spectrum

Observables are an average over height of $G_\ell$, in first approximation (depending on the assumption on the source), one can consider the Power Spectrum to be directly related to the Imaginary part of $G_\ell$,

$$ P_{\ell}(\omega) \propto |\mathrm{Im}(G_{\ell}(\omega)(r_0, s_0))|^2. \qquad (9) $$

We can compare the numerical power spectrum, i.e., using Green's function computed for harmonic degrees (modes) $\ell$ and frequencies $\omega/(2\pi)$, with observations.

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology. SIAM J. on Appl. Math., 2020.
# Power spectrum

We evaluate the imaginary part of the Green's function $G_ℓ(r = 1, s = 1)$ with frequencies and modes, using model S+Atmo.

$$ℓ \in (0, 1000), \quad ω/(2π) \in (0, 10)mHz.$$

<figure><img src="image_326.png" /></figure>
# Power spectrum

We evaluate the imaginary part of the Green's function $G_\ell(r = 1, s = 1)$ with frequencies and modes, using model S+Atmo.

Comparison with HMI data (white crosses), for $\ell \in (0, 100)$, $\omega/(2\pi) \in (0.5, 4)$mHz.

<figure><img src="image_327.png" /></figure>
# Plan

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

## Illustration: Green's function computations for time-distance helioseismology

* Modal Green's kernel

* Numerical experiments

* Power spectrum, comparison with measured data

* The vector wave problem
<figure><img src="image_329.png" /></figure>

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

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

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

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

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

Vector-wave equation with gravity

The same method can be applied to the vector-wave problem (Vector Spherical Harmonics),

$$ -\rho_0(\omega^2+2i\omega\Gamma)\xi - \nabla[\gamma\rho_0\nabla\cdot\xi] + (\nabla\rho_0)(\nabla\cdot\xi) - \nabla(\xi\cdot\nabla\rho_0) + (\xi\cdot\nabla)\nabla\rho_0 + \rho_0(\xi\cdot\nabla)\nabla\phi_0 = F, $$

That is, the vector Green's kernels (depending on the direction) can be obtained from the computation of only two boundary-value problems.

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

H. Barucq, F. Faucher, D. Fournier, L. Gizon and H. Pham, Efficient computation of modal Green's kernels for vectorial equations in helioseismology under spherical symmetry. Research Report, 2021.
<figure><img src="image_336.png" /></figure>

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

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

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

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

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

Towards 3D simulations

<figure><img src="image_342.png" /><figcaption>(a) spherical solar background model</figcaption></figure>

<figure><img src="image_343.png" /><figcaption>(b) Active region as velocity perturbation</figcaption></figure>

<figure><img src="image_344.png" /><figcaption>(c) Simulations with spherical background (top) and perturbations (bottom).</figcaption></figure>

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

<figure><img src="image_346.png" /></figure>
# Plan

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

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

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

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

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

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

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

Recap

* Simulation of wave propagation provides a non invasive tool for probing the invisible

* Advanced numerical methods and HPC are mandatory

* Inversion remains very sensitive

* Recovering real data is also very sensitive and mathematical modeling has its part to play
<figure><img src="image_354.png" /></figure>

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

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

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

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

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

Ongoing works

► Coupling SEM and DG for time-dependent problems in geophysics

► Conducting porous media

► Imaging 3D anisotropic elastic media with new acquisition methods (DAS)

► Solve the 3D Galbrun's equation with and without Cowling's approximation

► Construction of Butterfly diagrams for stars

► and so many very exciting problems

<figure><img src="image_360.png" /></figure>
<figure><img src="image_361.png" /></figure>

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

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

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

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

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

Makutu, September 2022

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

THANK YOU