# Méthode de décomposition de domaine basée sur une formulation de Trefftz pour l'acoustique anisotrope

Hélène Barucq, Ibrahima Djiba, Sébastien Tordeux

► To cite this version:

Hélène Barucq, Ibrahima Djiba, Sébastien Tordeux. Méthode de décomposition de domaine basée sur une formulation de Trefftz pour l'acoustique anisotrope. Smai 2023 - 11e Biennale française de Mathématiques Appliquées et Industrielles, May 2023, Le Gosier (Guadeloupe), France. hal-04112574

HAL Id: hal-04112574

https://inria.hal.science/hal-04112574

Submitted on 31 May 2023

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
# MÉTHODE DE DÉCOMPOSITION DE DOMAINE BASÉE SUR UNE FORMULATION DE TREFFTZ pour l'acoustique anisotrope

Hélène Barucq, Ibrahima Djiba, Sébastien Tordeux

EPI Makutu, INRIA, TotalEnergies, Université de Pau et des Pays de l'Adour, Bordeaux INP, CNRS

23 mai 2023 Guadeloupe
Le problème de Helmholtz

Le problème de Helmholtz s'écrit

$$
\begin{cases}
\Delta u + k^2 u = 0 & \text{dans } \Omega, \\
\frac{\partial u}{\partial n} + iku = \varphi & \text{sur } \partial\Omega.
\end{cases}
$$

Chercher $u \in H^1(\Omega)$ tel que pour tout $v \in H^1(\Omega)$

$$
\int_{\Omega} \nabla u \cdot \overline{\nabla v} - k^2 u \overline{v} + ik \int_{\partial \Omega} u \overline{v} = \int_{\partial \Omega} \varphi \overline{v}
$$

Grâce à cette forme **quasi-elliptique**, on a ainsi facilement

* le caractère bien posé (Alternative de Fredholm, prolongement unique)

• des méthodes d'éléments finis **convergentes**

Est-ce la fin de l'histoire ?

Oliver G. Ernst, and Martin J. Gander. "Why it is difficult to solve Helmholtz problems with classical iterative methods." Numerical analysis of multiscale problems (2012): 325-363.
Memory limit of a classic solver

$$Ax = F$$

- Second order system and continuous finite elements leads to less degrees of freedom.

- The Finite Element method deals correctly with the geometry and evanescent modes.

- Even efficient optimizations cannot overcome the memory issue.

- An iterative method ?

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

Memory cost (LU fact. from MUMPS$\mathcal{O}$) for different domain sizes $D_\Omega$ and approximation orders $r$. Memory = $(D_\Omega)^4$.

<table><thead><tr><th></th><th>Memory</th><th>Sockets</th><th>D<sub>Ω</sub></th><th>Power</th><th>Cost</th></tr></thead><tbody><tr><td>Laptop</td><td>32 Gb</td><td>1</td><td>10λ</td><td>200 W</td><td>0.034 €/h</td></tr><tr><td>// Computation</td><td>1 Tb</td><td>16</td><td>24λ</td><td>3.2 KW</td><td>0.544 €/h</td></tr><tr><td>HPC</td><td>320 Tb</td><td>5000</td><td>100λ</td><td>1.0 MW</td><td>170 €/h</td></tr><tr><td>HPC+</td><td>5120 Tb</td><td>80000</td><td>200λ</td><td>16 MW</td><td>2720 €/h</td></tr></tbody></table>
# Méthodes de Krylov

Théorème de Cayleigh Hamilton Si $\mathbf{A}$ est inversible, il existe un polynôme $p$ tel que

$$ \mathbf{A}^{-1} = p(\mathbf{A}). $$

## Méthode itérative de Krylov

$$ \mathbb{K}_{N_{\text{kry}}} = \text{span}(\{\mathbf{F}, \mathbf{A}\mathbf{F}, \mathbf{A}^2\mathbf{F}, \dots, \mathbf{A}^{N_{\text{kry}}-1}\mathbf{F}\}). $$

Trouver $X \in \mathbb{K}_{N_{\text{kry}}}$ qui minimise

$$ J(\mathbf{X}) = \|\mathbf{A}\mathbf{X} - \mathbf{F}\|_2^2 $$

Les méthodes **GMRES** sont **très lentement convergentes**

<figure><img src="image_3.png" /></figure>
Il existe des solutions partielles

- **Global preconditioners:**
  → A. Vion and C. Geuzaine. Double sweep preconditioner for optimized Schwarz methods applied to the Helmholtz problem. J. Comput. Phys., 2014.
  → M. Gander, L. Halpern, KS. L. Repiquet. A new coarse grid correction for RAS/AS. In Domain Decomposition Methods in Science and Engineering XXI 2014.

- **Domain decomposition methods:**
The impedance transmission condition

$$ egin{cases} Y p_T^{n+1} + \bar{v}_T^{n+1} \cdot \vec{n} = Y p_K + \bar{v}_K^n \cdot \vec{n}, \\ Y p_K^{n+1} - \bar{v}_K^{n+1} \cdot \vec{n} = Y p_T^n - \bar{v}_T^n \cdot \vec{n}, \end{cases} $$

→ A. Modave, X. Antoine et C. Geuzaine. An efficient domain decomposition method with cross-point treatment for Helmholtz problems. In CSMA 2019-14e Colloque National en Calcul des Structures. 2019.

→ F. Collino, P. Joly, and M. Lecouvez. "Exponentially convergent non overlapping domain decomposition methods for the Helmholtz equation." ESAIM: Math. Mod. and Num. An. 54.3. 775-810. 2020.

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

- **Les méthodes de Trefftz**
  → Retour à la nature hyperbolique
  → Naturellement adaptée aux méthodes GMRES
* O. Cessenat, and B. Després. "Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem." SIAM journal on numerical analysis. 1998.

* P. Monk, and D. Q. Wang. A least-squares method for the Helmholtz equation. Computer methods in applied mechanics and engineering. 1999.

* A. Moiola. Trefftz-discontinuous Galerkin methods for time-harmonic wave problems. Diss. ETH Zurich. 2011.

* R. Hiptmair, A. Moiola, and I. Perugia. "A survey of Trefftz methods for the Helmholtz equation." Building bridges: connections and challenges in modern approaches to numerical partial differential equations. Springer, 2016.

* M. Sirdey, Méthode itérative de type Trefftz pour la simulation d'ondes électromagnétiques en trois dimensions, diss. UPPA, 2022.
# Introduction

* Qu'est ce qu'une méthode de Trefftz variationnelle

* Méthode de Galerkin discontinue

* Fonctions de base solutions de l'EDP

* Ondes planes, Fonctions de Bessel, ...

* Objectifs ?

* Implémentation une **méthode de Trefftz**

* Optimiser les flux numériques (solveur de Riemann)

* Inverser la matrice à l'aide d'une **méthode GMRES** préconditionnée

* Pourquoi ?

* Pour gagner en **précision**

* Pour réduire la **taille mémoire**

* Pour réduire le **temps de calcul**

* Pour quel type de problème ?

* **Problème linéaire**

* Paramètres physiques constants par élément

* Anisotrope

* En régime fréquentile
<figure><img src="image_5.png" /></figure>

$$
\left\{
\begin{array}{ll}
A \nabla \mathbf{p}(\vec{x}, t) = -\rho_0 \frac{\partial \vec{v}}{\partial t}(\vec{x}, t) & \text{dans } \Omega, \\
\text{div } \vec{v}(\vec{x}, t) = -\xi \frac{\partial \mathbf{p}}{\partial t}(\vec{x}, t) + & \text{dans } \Omega, \\
Y \mathbf{p} - \rho_0 \vec{v} \cdot \vec{n} = \phi & \text{sur } \partial\Omega.
\end{array}
\right.
$$

Les paramètres

• $\rho_0$ la masse volumique

• $\xi$ une constante thermodynamique

$$
\xi = \frac{1}{\rho_0} \frac{d\rho}{d\mathbf{p}}(\rho_0)
$$

• c la vitesse du son

$$
c = \frac{1}{\sqrt{\xi \rho_0}}.
$$

• Y est l'admittance

Les inconnues

• $\vec{v}$ la vitesse

• p la pression

Anisotropie : A matrice sdp
# Le système normalisé

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

Le système normalisé est donné par

$$ egin{cases} A \nabla p(\vec{x}, t) = -\frac{1}{c} \frac{\partial \vec{v}}{\partial t}(\vec{x}, t) & \text{dans } \Omega, \\ \text{div } \vec{v}(\vec{x}, t) = -\frac{1}{c} \frac{\partial p}{\partial t}(\vec{x}, t) & \text{dans } \Omega, \\ Y p - \vec{v} \cdot \vec{n} = \varphi & \text{sur } \partial \Omega. \end{cases} $$

avec

$$ p(\vec{x}, t) = \bar{p}_0 p(x, t), \quad \vec{v}(\vec{x}, t) = \bar{v}_0 \vec{v}(x, t), \phi(\vec{x}, t) = \bar{\phi}_0 \varphi(x, t), \\ \text{et } \quad \frac{\bar{p}_0}{\rho_0 \bar{v}_0} = c, \quad \bar{\phi}_0 = \bar{p}_0. $$

* $\bar{p}_0$ la pression caractéristique,

* $\bar{v}_0$ la vitesse caractéristique,

* $\bar{\phi}_0$ l'amplitude caractéristique de la source,

* $p$ la pression normalisée,

* $\vec{v}$ la vitesse normalisée,

* $\varphi$ la source normalisée.
Le système en régime harmonique

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

$$
\left\{
\begin{array}{l}
A \nabla p(\vec{x}) = ik \vec{v}(\vec{x}) \text{ dans } \Omega, \\
\text{div } \vec{v}(\vec{x}) = ikp(\vec{x}) \text{ dans } \Omega, \\
Yp - \vec{v} \cdot \vec{n} = \varphi \text{ sur } \partial\Omega,
\end{array}
\right.
$$

avec $p : \Omega \to \mathbb{C}^2$, $\vec{v} : \Omega \to \mathbb{C}^2$, et $\varphi : \Omega \to \mathbb{C}^2$

$$
\left\{
\begin{array}{l}
p(\vec{x}, t) = \operatorname{Re}(p(x) \exp(-i\omega t)), \\
\vec{v}(\vec{x}, t) = \operatorname{Re}(\vec{v}(\vec{x}) \exp(-i\omega t)), \\
\varphi(\vec{x}, t) = \operatorname{Re}(\varphi(\vec{x}) \exp(-i\omega t)),
\end{array}
\right.
$$

• Le nombre d'onde $k = \frac{\omega}{c}$,

• La matrice $A$ est la matrice d'anisotropie
Maillage du domaine

• $K \in \mathcal{K}$ un élément du maillage

• $\partial K$ sa frontière

• $\mathcal{K}$ ensemble des éléments

• $\bar{\Omega} = \bigcup_{K \in \mathcal{K}} \bar{K}$.

<figure><img src="image_8.png" /><figcaption>Figure: Discrétisation du domaine $\Omega$.</figcaption></figure>
Maillage du domaine

$$
\partial K = \bigcup_{F \in \mathcal{F}_K} F
$$

• $\mathcal{F}_K$ ensemble des faces de $K$

• F une face

• $\mathcal{F}_{\text{int}}$ ensemble des **faces intérieures**

• $\mathcal{F}_{\text{ext}}$ ensemble des faces extérieures

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

• Si $F$ est une **face intérieure**, alors elle touche deux éléments $K_F^+$ et $K_F^-$.

• Si $F$ est une **face extérieure**, alors elle ne touche qu'un seul élément $K_F$.
L'espace de Trefftz

$$ \mathcal{X} = \prod_{K \in \mathcal{K}} \mathcal{X}_K $$

avec $(p_K, \vec{v}_K)$ la restriction sur $K$ de $(p, \vec{v})$ et

$$ \mathcal{X}_K = \left\{ (p_K, \vec{v}_K) \in H^1(K) \times H(\mathrm{div}, K) \mid \mathrm{div}(\vec{v}_K) = ikp_K \text{ dans } K \text{ et } A\nabla p_K = ik\vec{v}_K \text{ dans } K \text{ et } \vec{v}_K \cdot \vec{n}_K \in L^2(\partial K) \right\} $$

* solutions **locales** acoustiques

* discontinues entre deux éléments (méthode de Galerkin discontinue)
Le principe des travaux virtuels pour Trefftz

Soit $u = (p, \vec{v})$ la solution et $u' = (p', \vec{v}')$ la fonction test.

Le travail virtuel prend la forme

$$W_K = ik \int_K p_K \overline{p'_K} + \vec{v}_K \cdot A^{-1} \overline{\vec{v}'_K}.$$

La solution du système acoustique vérifie

$$ \left\{ \begin{array}{l} A \nabla p_K = ik \vec{v}_K \text{ dans } K, \\ \text{div } \vec{v}_K = ik p_K \text{ dans } K, \end{array} \right. \qquad W_K = \int_K \text{div } \vec{v}_K \overline{p'_K} + \nabla p_K \cdot \overline{\vec{v}'_K} $$
Le principe des travaux virtuels pour Trefftz

Soit $\mathcal{U} = (p, \vec{v})$ la solution et $\mathcal{U}' = (p', \vec{v}')$ la fonction test.

Le travail virtuel prend la forme

$$W_K = ik \int_K p_K \overline{p'_K} + \vec{v}_K \cdot A^{-1} \overline{\vec{v}'_K}$$

La solution du système acoustique vérifie

$$\begin{cases} A \nabla p_K = ik \vec{v}_K \text{ dans } K, \\ \operatorname{div} \vec{v}_K = ik p_K \text{ dans } K, \end{cases}$$

$$W_K = \int_K \operatorname{div} \vec{v}_K \overline{p'_K} + \nabla p_K \cdot \overline{\vec{v}'_K}$$

$$\begin{cases} A \nabla p'_K = ik \vec{v}'_K \text{ dans } K, \\ \operatorname{div} \vec{v}'_K = ik p'_K \text{ dans } K, \end{cases}$$

$$W_K = - \int_K p_K \operatorname{div} \overline{\vec{v}'_K} + \vec{v}_K \cdot \nabla \overline{p'_K}$$
Le principe des travaux virtuels pour Trefftz

Soit $\mathcal{U} = (p, \vec{v}')$ la solution et $\mathcal{U}' = (p', \vec{v}')$ la fonction test.

Le travail virtuel prend la forme

$$W_K = ik \int_K p_K \overline{p'_K} + \overrightarrow{v}_K \cdot A^{-1} \overline{\overrightarrow{v}'_K}$$

La solution du système acoustique vérifie

$$\left\{ \begin{array}{l} A \nabla p_K = ik \overrightarrow{v}_K \text{ dans } K, \\ \text{div } \overrightarrow{v}_K = ik p_K \text{ dans } K, \end{array} \right.$$

$$W_K = \int_K \operatorname{div} \overrightarrow{v}_K \overline{p'_K} + \nabla p_K \cdot \overline{\overrightarrow{v}'_K}$$

$$\left\{ \begin{array}{l} A \nabla p'_K = ik \overrightarrow{v}'_K \text{ dans } K, \\ \text{div } \overrightarrow{v}'_K = ik p'_K \text{ dans } K, \end{array} \right.$$

$$W_K = - \int_K p_K \operatorname{div} \overrightarrow{v}'_K + \overrightarrow{v}_K \cdot \nabla p'_K$$

En soustrayant les deux expressions, nous avons

$$\int_K \operatorname{div} \overrightarrow{v}_K \overline{p'_K} + \nabla p_K \cdot \overline{\overrightarrow{v}'_K} + p_K \operatorname{div} \overrightarrow{v}'_K + \overrightarrow{v}_K \cdot \nabla \overline{p}'_K = 0.$$
Le principe des travaux virtuels pour Trefftz

Soit $\mathcal{U} = (p, \vec{v})$ la solution et $\mathcal{U}' = (p', \vec{v}')$ la fonction test.

Le travail virtuel prend la forme

$$W_K = ik \int_K p_K \overline{p'_K} + \vec{v}_K \cdot A^{-1} \overline{\vec{v}'_K}$$

La solution du système acoustique vérifie

$$\left\{ \begin{array}{l} A \nabla p_K = ik \vec{v}_K \text{ dans } K, \\ \text{div } \vec{v}_K = ik p_K \text{ dans } K, \end{array} \right. \qquad W_K = \int_K \text{div } \vec{v}_K \overline{p'_K} + \nabla p_K \cdot \overline{\vec{v}'_K}$$

$$\left\{ \begin{array}{l} A \nabla p'_K = ik \vec{v}'_K \text{ dans } K, \\ \text{div } \vec{v}'_K = ik p'_K \text{ dans } K, \end{array} \right. \qquad W_K = - \int_K p_K \text{div } \overline{\vec{v}'_K} + \vec{v}_K \cdot \nabla p'_K$$

En soustrayant les deux expressions, nous avons

$$\int_K \text{div } \vec{v}_K \overline{p'_K} + \nabla p_K \cdot \overline{\vec{v}'_K} + p_K \text{div } \overline{\vec{v}'_K} + \vec{v}_K \cdot \nabla \overline{p'_K} = 0.$$

D'après le théorème de divergence, nous avons le **principe de réciprocité**

$$\int_{\partial K} \vec{v}_K \cdot \vec{n}_K \overline{p'_K} + p_K \overline{\vec{v}'_K} \cdot \vec{n}_K = 0,$$
Principe de réciproité

On somme sur tous les éléments

$$
\sum_{K \in \mathcal{K}} \int_{\partial K} \overrightarrow{v}_K \cdot \vec{n}_K \overline{p'_K} + p_K \overline{v}'_K \cdot \vec{n}_K = 0,
$$

On réécrit cette expression en termes de sauts

$$
\sum_{F \in \mathcal{F}_{\text{int}}} \int_F \overrightarrow{v}_F \overline{[p']}_F + p_F \overline{[\overrightarrow{v}']}_F + \sum_{F \in \mathcal{F}_{\text{ext}}} \int_F \overrightarrow{v}_F \overline{p'}_{K_F} \overrightarrow{n}_F + p_F \overline{\overrightarrow{v}'}_{K_F} \cdot \overrightarrow{n}_F = 0,
$$

où

• $p_F$ et $\overrightarrow{v}_F$ est la trace de $p$ et $\overrightarrow{v}$ sur $F$,

• $\vec{n}_F$ est la normale sortante à la face extérieure $F$.

Les sauts sont définis par

$$
\overline{[\overrightarrow{v}]}_F = \overrightarrow{v}_{K_F^+} \cdot \vec{n}_{K_F^+} + \overrightarrow{v}_{K_F^-} \cdot \vec{n}_{K_F^-}, \quad \overline{[p]}_F = p_{K_F^+} \vec{n}_{K_F^+} + p_{K_F^-} \vec{n}_{K_F^-}
$$

<figure><img src="image_10.png" /></figure>
## Méthodes des flux upwind

On va remplacer $\vec{v}_F$ et $p_F$

$$ \sum_{F \in \mathcal{F}_{\text{int}}} \int_F \vec{v}_F [\overrightarrow{p'}]_F + p_F [\overrightarrow{\vec{v}'}]_F + \sum_{F \in \mathcal{F}_{\text{ext}}} \int_F \vec{v}_F \overrightarrow{p'}_{K_F} \vec{n}_F + p_F \overrightarrow{\vec{v}'}_{K_F} \cdot \vec{n}_F = 0, $$

par des traces numériques qui seront construite à partir des flux sortants et entrants

entrant : $Y_K p_K - \vec{v}_K \cdot \vec{n}_K$,

sortrant : $Y_K p_K + \vec{v}_K \cdot \vec{n}_K$.

A l'interface entre deux éléments, on a la continuité de $p$ et $\vec{v} \cdot \vec{n}$

$$ p_F = \frac{1}{Y_T + Y_K} (Y_K p_K + \vec{v}_K \cdot \vec{n}_K) + \frac{1}{Y_T + Y_K} (Y_T p_T + \vec{v}_T \cdot \vec{n}_T). $$

On obtient une notion de trace définie pour les fonctions discontinues

$$ \widehat{p_F} = \frac{Y_T p_K + Y_K p_T}{Y_T + Y_K} + \frac{Y_T Y_K}{Y_T + Y_K} [\overrightarrow{v}]_F, $$

On fait de même pour $\widehat{v}_F$ et pour les faces de bord.

Pour $Y_T$, on utilise l'admitance du milieu infini :

$$ Y_T = \frac{1}{(\vec{n} \cdot \mathbf{A}^T \vec{n})^{\frac{1}{2}}} $$
# Flux de Riemann intérieur

On résout le problème acoustique posé sur tout $\mathbb{R}^2$

muni des conditions initiales

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

$$ egin{cases} p(\vec{x}, 0) = p_T & \text{pour } x_1 < 0, \\ p(\vec{x}, 0) = p_K & \text{pour } x_1 > 0, \\ \vec{v}(\vec{x}, 0) = \vec{v}_T & \text{pour } x_1 < 0, \\ \vec{v}(\vec{x}, 0) = \vec{v}_K & \text{pour } x_1 > 0. \end{cases} $$

et de condition de transmission en $x_1 = 0$.

• Pour $x_1 < 0$, on a

$$ egin{cases} A_T \nabla p(\vec{x}, t) = -\frac{1}{c} \frac{\partial \vec{v}}{\partial t}(\vec{x}, t), \\ \text{div } \vec{v}(\vec{x}, t) = -\frac{1}{c} \frac{\partial p}{\partial t}(\vec{x}, t), \end{cases} $$

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

• Pour $x_1 > 0$, on a

$$ egin{cases} A_K \nabla p(\vec{x}, t) = -\frac{1}{c} \frac{\partial \vec{v}}{\partial t}(\vec{x}, t), \\ \text{div } \vec{v}(\vec{x}, t) = -\frac{1}{c} \frac{\partial p}{\partial t}(\vec{x}, t), \end{cases} $$

On retrouve alors la même condition

$$ p(0, t) = \frac{Y_{TPK} + Y_{KPT}}{Y_T + Y_K} + \frac{Y_T Y_K}{Y_T + Y_K} [\vec{v}]_F, $$

On fait de même pour $\vec{v}$.
La formulation variationnelle

Chercher $\mathcal{U} \in \mathcal{X}$ tel que pout tout $\mathcal{U}' \in \mathcal{X}$

$$
\mathbf{a}(\mathcal{U}, \mathcal{U}') = \ell(\mathcal{U}'),
$$

avec

$$
\left\{
\begin{array}{l}
\mathbf{a}(\mathcal{U}, \mathcal{U}') = \sum_{F \in \mathcal{F}_{\text{int}}} \mathbf{a}_{F}^{\text{int}}(\mathcal{U}, \mathcal{U}') + \sum_{F \in \mathcal{F}_{\text{ext}}} \mathbf{a}_{F}^{\text{ext}}(\mathcal{U}, \mathcal{U}'), \\
\ell(\mathcal{U}') = \sum_{F \in \mathcal{F}_{\text{ext}}} \ell_{F}^{\text{ext}}(\mathcal{U}'),
\end{array}
\right.
$$

et

$$
\left\{
\begin{array}{l}
\mathbf{a}_{F}^{\text{int}}(\mathcal{U}, \mathcal{U}') \\
\mathbf{a}_{F}^{\text{ext}}(\mathcal{U}, \mathcal{U}') - \ell_{F}^{\text{ext}}(\mathcal{U}')
\end{array}
\right.
=
\begin{array}{l}
\displaystyle \int_{F} \overline{\widehat{\mathcal{V}}_{F}} \overline{[\overline{\mathcal{P}'}]_{F}} + \widehat{p_{F}} \overline{[\overline{\mathcal{V}'}]_{F}}, \\
\displaystyle \int_{F} \overline{\widehat{\mathcal{V}}_{F}} \overline{p'_{K_{F}}} \overline{\mathcal{n}}_{K_{F}} + \widehat{p_{F}} \overline{\overline{\mathcal{V}'}_{K_{F}}} \cdot \overline{\mathcal{n}}_{K_{F}}.
\end{array}
$$

La forme variationnelle obtenue est

• à partie réelle positives

$$
\mathrm{Re} \left( \mathbf{a}(\mathcal{U}, \mathcal{U}) \right) \geq 0.
$$

• définie

$$
\mathbf{a}(\mathcal{U}, \mathcal{U}) = 0 \quad \Rightarrow \quad \mathcal{U} \equiv 0.
$$

Théorie de convergence des méthodes de Trefftz :

• Cessenat et Després 98, Monk et Wank 99
Les ondes planes acoustiques anisotropes

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

L'espace de Trefftz discret est engendré par la base

$$
u_{K,n}(\vec{x}) = 0 \text{ si } \vec{x} \notin K
$$

$$
u_{K,n}(\vec{x}) = u_{(\cos(n\delta\theta), \sin(n\delta\theta))}(\vec{x}) \text{ si } \vec{x} \in K \text{ avec } \delta\theta = \frac{2\pi}{N}.
$$

On a

$$
u \in \mathbb{X}_N \iff u(\vec{x}) = \sum_{K \in \mathcal{K}} \sum_{n=1}^{N} u_{K,n} u_{K,n}(\vec{x})
$$
Quel est le nombre optimal d'ondes planes ?

* Moiola, A. (2011). Trefftz-discontinuous Galerkin methods for time-harmonic wave problems (Doctoral dissertation, ETH Zurich). (Page 156).

* Moiola, A., Hiptmair, R., Perugia, I. (2011). Vekua theory for the Helmholtz operator. Zeitschrift für angewandte Mathematik und Physik, 62(5), 779-807.

Réduction de base → amélioration du conditionnement

* T. Luostari, T. Huttunen, P. Monk, (2013). Improvements for the ultra weak variational formulation. Int. J. Numer. Meth. Engng 94, 598624.

* S. Congreve, J. Gedicke, I. Perugia, (2019). Numerical investigation of the conditioning for plane wave discontinuous Galerkin methods, Vol. 126 of Lecture Notes in Computational Science and Engineering, Springer.

* H. Barucq, A. Bendali, J. Diaz, S. Tordeux, (2021). Local strategies for improving the conditioning of the plane-wave Ultra-Weak Variational Formulation. Journal of Computational Physics, 441, 110449.

Des alternatives aux ondes planes existent

* Modes évancents

* E. Parolin, D. Huybrechs, A. MOIOLA, (2022). Stable approximation of Helmholtz solutions by evanescent plane waves. arXiv preprint arXiv:2202.05658.

* Méthodes Quasi-Trefftz

* LM. Imbert-Gerard, G. Sylvand, Three types of quasi-Trefftz functions for the 3D convected Helmholtz equation: construction and approximation properties, preprint

* H. S. Fure, S. Pernet, M. Sirdey, S. Tordeux, (2020). A discontinuous Galerkin Trefftz type method for solving the two dimensional Maxwell equations. SN Partial Differ. Equ. Appl. 1, 23.
Le squelette de la matrice

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

<table><tr><td>K<sub>13</sub></td><td>K<sub>14</sub></td><td>K<sub>15</sub></td><td>K<sub>16</sub></td></tr><tr><td>K<sub>9</sub></td><td>K<sub>10</sub></td><td>K<sub>11</sub></td><td>K<sub>12</sub></td></tr><tr><td>K<sub>7</sub></td><td>K<sub>6</sub></td><td>K<sub>7</sub></td><td>K<sub>8</sub></td></tr><tr><td>K<sub>1</sub></td><td>K<sub>2</sub></td><td>K<sub>3</sub></td><td>K<sub>4</sub></td></tr></table>

<table><tr><td>K<sub>13</sub></td><td>K<sub>14</sub></td><td>K<sub>15</sub></td><td>K<sub>16</sub></td></tr><tr><td>K<sub>9</sub></td><td>K<sub>10</sub></td><td>K<sub>11</sub></td><td>K<sub>12</sub></td></tr><tr><td>K<sub>6</sub></td><td>K<sub>7</sub></td><td>K<sub>8</sub></td><td>K<sub>8</sub></td></tr><tr><td>K<sub>1</sub></td><td>K<sub>2</sub></td><td>K<sub>3</sub></td><td>K<sub>4</sub></td></tr></table>

* Les fonctions de base de l'élt 1 interagissent avec celles des élts 1, 2 et 5

* Les fonctions de base de l'élt 6 interagissent avec celles des élts 2, 5, 6, 7 et 10.
Soit $\Omega = [0, 4]^2$

$$
\left\{
\begin{array}{l}
\nabla p(\vec{x}, t) = -\frac{1}{c} \frac{\partial \vec{v}}{\partial t}(\vec{x}, t), \quad \text{dans } \Omega, \\
\\
\text{div } \vec{v}(\vec{x}, t) = -\frac{1}{c} \frac{\partial p}{\partial t}(\vec{x}, t), \quad \text{dans } \Omega, \\
\\
p - \vec{v} \cdot \vec{n} = p_{\text{inc}} - \vec{v}_{\text{inc}} \cdot \vec{n} \quad \text{sur } \partial\Omega,
\end{array}
\right.
$$

avec comme champ incident un champ monopo-
laire

$$
p_{\text{inc}}(\vec{x}) = \frac{iH_0^{(1)}(k \| \vec{x} - \vec{x}^0 \|)}{4},
$$

$$
\vec{v}_{\text{inc}}(\vec{x}) = - \frac{H_1^{(1)}(k \| \vec{x} - \vec{x}^0 \|)}{4} \frac{\vec{x} - \vec{x}^0}{\| \vec{x} - \vec{x}^0 \|}
$$

avec

$$
\vec{x}^0 = (5, -1)^T,
$$

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

point source en $x = 5$ et $y = -1$
Surf de la solution numérique, de la solution exacte et de l'erreur

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

• La solution numérique

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

• La solution exacte

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

• L'erreur

• $\Omega = [0, 4]^2$

• $\vec{x}^0 = (5, -1)^T$

• 18 éléments en x comme en y

• 6 fonctions de base par élément
Courbes d'erreurs en norme $L^2$ de $U_h - U$ en fonction du pas de maillage $h$

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

• $\Omega = [0 \ 4]^2$

• $\frac{\partial^2}{\partial x^2} = (10, 9)^T$

• $N_b = 12, 14$ et
16 ondes planes
fixées

• $Y = 1$

• Il est inutile de raffiner trop le maillage.
Courbes d'erreurs en norme $L^2$ de $U_h - U$ en fonction du nombre de fonctions de base par élément

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

* $\Omega = [0 \ 4]^2$

* $\bar{x}^0 = (10, 9)^T$

* On a fixé le pas à $h = 0.25, 0.5$ et $1$

* $Y = 1$

* Il est inutile de trop augmenter le nombre de fonctions tests.
# Décomposition de Block Jacobi de la matrice (Cessenat-Després)

O. Cessenat and B. Despres, (1998). Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem. SIAM journal on numerical analysis, 35(1), 255-299.

$$A = M + N$$

M est bloc-diagonale

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

norme UWVF

$$||U||^2_{UWVF} = \text{Re} \begin{pmatrix} U^* & M & U \end{pmatrix}$$

N les autres blocs

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

norme DG

$$||U||^2_{DG} = \text{Re} \begin{pmatrix} U^* & A & U \end{pmatrix}$$
# Résolution itérative du système linéaire $AU = F$

Pourquoi ? Pour gagner en mémoire.

* Résolution itérative de block-Jacobi $MU_{n+1} = -NU_n + F$

* Méthode GMRES $U = GMRES(A, F)$

* Méthode GMRES préconditionnée par bloc-Jacobi $U = GMRES(M^{-1}A, M^{-1}F)$

<figure><img src="image_23.png" /><figcaption>Figure 1</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 2</figcaption></figure>

* $\Omega = [0, 4]^2$, $\vec{x}^0 = (10, 9)^T$

* 18 éléments en x et en y

* `nkry = 15` comme restart

* $Y = 1$ pour figure 1

* $Y = \frac{1}{19}$ pour la figure 2

* 6 fonctions de base par élément
Nous avons aussi fait du 3D

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

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

<figure><img src="image_27.png" /></figure>
## Passage à l'échelle

* Résultat extrait de la thèse de M. Sirdey (Doctorante ONERA-Makutu) codirigée avec S. Pernet (ONERA)

* Diffraction électromagnétique d'une frégate illuminée par le haut

$$300 \times 100 \times 100 \lambda^3.$$

* Calcul mettant en jeu plus d'un milliard d'inconnues

<table><thead><tr><td>Inconnue</td><td>Matrice</td><td>Factorisation LU</td><td>La méthode</td></tr></thead><tbody><tr><td>16 Go</td><td>3 To</td><td>1000 To</td><td>800 Go</td></tr></tbody></table>

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

* Local strategies for improving the conditioning of the plane-wave Ultra-Weak Variational Formulation Hélène Barucq, Abderrahmane Bendali, Julien Diaz, Sébastien Tordeux Journal of Computational Physics, 2021, 441, pp.110449

* Ultra-weak variational formulation for heterogeneous maxwell problem in the context of high performance computing Sébastien Pernet, Margot Sirdey, Sébastien Tordeux 2022, Esaim Proc.

* A discontinuous Galerkin Trefftz type method for solving the two dimensional Maxwell equations Hakon Sem Fure, Sébastien Pernet, Margot Sirdey, Sébastien Tordeux SN Partial Differential Equations and Applications, 2020, 1 (23), pp.19.

* A Symmetric Trefftz-DG formulation based on a local boundary element method for the solution of the Helmholtz equation Hélène Barucq, Abderrahmane Bendali, M'Barek Fares, Vanessa Mattesi, Sébastien Tordeux Journal of Computational Physics, 2017, 330, pp.1069-1092.

* GoTem3 : Un code opensource fortran Trefftz Maxwell 3D.