# Solveurs itératifs pour des méthodes de Trefftz par ondes planes évanescentes pour l'équation de Helmholtz

Vincent Robert

► To cite this version:

Vincent Robert. Solveurs itératifs pour des méthodes de Trefftz par ondes planes évanescentes pour l'équation de Helmholtz. Numerical Analysis [cs.NA]. 2024. hal-05007030

HAL Id: hal-05007030

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

Submitted on 26 Mar 2025

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
Solveurs itératifs pour des méthodes de Trefftz
par ondes planes évanescentes
pour l'équation de Helmholtz

Rapport de stage

Vincent Robert

sous l'encadrement d'Emile Parolin
Centre INRIA Paris - Equipe Alpines

27 septembre 2024
Je remercie très chaleureusement Emile Parolin pour la confiance qu'il m'a accordée,
sa grande disponibilité et ses patientes explications.

Merci aussi à Caroline et mes parents pour leur soutien indéfectible.
Table des matières

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1 Introduction</td><td>1</td></tr><tr><td>1.1 Contexte</td><td>1</td></tr><tr><td>1.2 Objectifs du stage</td><td>1</td></tr><tr><td>1.3 Contenu</td><td>1</td></tr><tr><td>1.4 Travail d'implémentation</td><td>1</td></tr><tr><td>2 Méthode numérique</td><td>2</td></tr><tr><td>2.1 Equation de Helmholtz</td><td>2</td></tr><tr><td>2.2 Méthodes de Trefftz</td><td>3</td></tr><tr><td>2.3 Formulation ultra faible</td><td>4</td></tr><tr><td>2.4 Espaces d'approximation</td><td>7</td></tr><tr><td>2.5 Opérateur T</td><td>10</td></tr><tr><td>2.6 Dérivation de la formulation avec Λ.</td><td>11</td></tr><tr><td>3 Implémentation du système linéaire</td><td>13</td></tr><tr><td>3.1 Implémentation des matrices</td><td>13</td></tr><tr><td>3.2 Vérification de la structure des matrices</td><td>16</td></tr><tr><td>4 Etude spectrale</td><td>17</td></tr><tr><td>4.1 Propriétés de D</td><td>17</td></tr><tr><td>4.2 Valeurs propres de D<sup>-1</sup>C</td><td>18</td></tr><tr><td>4.3 SVD régularisée</td><td>20</td></tr><tr><td>4.4 Valeurs propres de D<sub>ε</sub><sup>-1</sup>C</td><td>21</td></tr><tr><td>5 Algorithmes séquentiels de résolution</td><td>24</td></tr><tr><td>5.1 Résolution directe</td><td>24</td></tr><tr><td>5.2 Algorithme de Richardson</td><td>24</td></tr><tr><td>5.3 GMRES</td><td>25</td></tr><tr><td>6 Expériences numériques, algorithmes séquentiels</td><td>28</td></tr><tr><td>6.1 Tests de validité</td><td>28</td></tr><tr><td>6.2 Tests de convergence en N.</td><td>31</td></tr><tr><td>6.3 Tests de convergence en h.</td><td>39</td></tr><tr><td>6.4 Vers le régime haute fréquence : Gmres, nombre d'ondes planes fixé, k croissant</td><td>42</td></tr><tr><td>7 Algorithmes parallèles de résolution</td><td>46</td></tr><tr><td>7.1 Présentation</td><td>46</td></tr><tr><td>7.2 Construction des matrices et vecteurs</td><td>47</td></tr><tr><td>7.3 Assemblage du système linéaire</td><td>48</td></tr><tr><td>8 Expériences numériques, algorithmes parallèles</td><td>50</td></tr><tr><td>8.1 Test de validation</td><td>50</td></tr><tr><td>8.2 Weak scaling</td><td>50</td></tr><tr><td>8.3 strong scaling</td><td>52</td></tr><tr><td>8.4 Une recherche de précision maximale</td><td>53</td></tr><tr><td>8.5 Somme aléatoire de Fourier-Bessels</td><td>54</td></tr></tbody></table>
<table><tr><td>9</td><td>Problème de diffraction sur un écran</td><td>56</td></tr><tr><td>9.1</td><td>Présentation - problème aux limites et formulation variationnelle</td><td>56</td></tr><tr><td>9.2</td><td>Maillages et espaces locaux</td><td>57</td></tr><tr><td>9.3</td><td>Implémentation et premiers tests de validation</td><td>57</td></tr><tr><td>9.4</td><td>Tests de convergence</td><td>59</td></tr><tr><td>9.5</td><td>Test 1</td><td>61</td></tr><tr><td>9.6</td><td>Test 2</td><td>65</td></tr><tr><td>9.7</td><td>Approximation d'une singularité de coin</td><td>67</td></tr><tr><td>10</td><td>Déflation</td><td>71</td></tr><tr><td>10.1</td><td>Présentation</td><td>71</td></tr><tr><td>10.2</td><td>Expériences numériques</td><td>72</td></tr><tr><td>11</td><td>Conclusion</td><td>77</td></tr></table>
# 1 Introduction

## 1.1 Contexte

Ce rapport documente le travail réalisé lors de mon stage de fin de Master 2 de Mathématiques de la modélisation, qui s'est déroulé au centre INRIA de Paris entre mai et septembre 2024, au sein de l'équipe Alpines. Les recherches de cette équipe se concentrent sur le calcul scientifique haute-performance, et comprennent le développement de méthodes numériques pour la résolution d'équations aux dérivées partielles comme d'algorithmes pour l'algèbre linéaire creuse et dense (voir [1]). Au sein de l'équipe, les recherches d'Emile Parolin se concentrent sur l'étude mathématique de phénomènes de propagation d'ondes acoustiques et électromagnétiques à l'aide de différentes méthodes numériques (voir [2]).

## 1.2 Objectifs du stage

L'objet du stage est d'explorer une nouvelle méthode de Trefftz pour l'approximation de solution de l'équation de Helmholtz sur un domaine borné en deux dimensions, dans l'objectif de parvenir à effectuer, aussi efficacement que possible, à l'aide de méthodes numériques directes comme itératives, des simulations en régime haute fréquence, problème notoirement difficile.

## 1.3 Contenu

Les nouveaux espaces d'approximation proposés par Emile Parolin, Daan Huybrechs et Andrea Moiola dans [13], n'ont pas encore été utilisés dans des formulations de Trefftz habituelles (UWVF, TDG, etc...). Ce travail explore pour la première fois ces espaces avec la formulation variationnelle ultra-faible présentée par Bruno Desprès et Olivier Cessenat dans [5].

Comme habituellement pour les méthodes de Trefftz, on observe des instabilités numériques qui limitent la méthode. On propose une technique de régularisation similaire à [4] basée sur une SVD tronquée. Une analyse spectrale empirique valide l'approche (chapitre 4).

On s'intéresse ensuite aux méthodes de résolution des systèmes linéaires associés, notamment GMRES. On propose et implémente des versions séquentielle et parallèle (chapitres 5-6-7-8), et l'on conduit une analyse empirique paramétrique par rapport aux principaux paramètres de la méthode : raffinement du maillage, augmentation des espaces locaux, nombre de threads (pour la version parallèle). On montre ainsi la capacité de cette méthode à approcher efficacement des problèmes haute-fréquence, et atteindre un degré de précision élevé dépendant de $\kappa$, de $N$ et du maillage.

On expérimente la méthode sur un problème de diffraction (chapitre 9).

On présente enfin certains des axes de recherches menés qui n'ont pas abouti : la généralisation des conditions de transmission aux interfaces dans la formulation variationnelle, ainsi qu'une méthode de déflation (chapitre 10).

## 1.4 Travail d'implémentation

Une part importante du travail produit consiste en du code, écrit en langage *Julia*. Emile Parolin a implémenté ces nouveaux espaces d'approximation dans ce langage préalablement à notre stage, nous avons donc contribué au développement de ce code, en y ajoutant plusieurs algorithmes itératifs, plusieurs modes de construction des éléments des systèmes linéaires ainsi que plusieurs types de conditions aux limites, et enfin de nombreuses fonctions utilitaires.

Il n'est malheureusement pas encore public, nous ne pouvons fournir d'adresse où le consulter en ligne. Notons ici que ce stage a été pour nous l'occasion d'apprendre ce langage, et que cet apprentissage a constitué une part non négligeable de notre travail.
# 2 Méthode numérique

## 2.1 Equation de Helmholtz

Nous considérons ici un domaine borné $\Omega \in \mathbb{R}^2$, de bord $\Gamma$ (de normale $\nu$ sortante par rapport à $\Omega$). Nous partons de l'équation des ondes scalaires décrivant la propagation du son en milieu homogène dans $\mathbb{R}^2$, prenant la forme

$$ \frac{\partial^2 p}{\partial t^2} - c^2 \Delta p = 0, \quad x \in \mathbb{R}^2, $$

où $c > 0$ est la célérité de l'onde et $p = p(x, t)$ sa pression inconnue. L'équation est complétée de conditions aux limites classiques, aboutissant par exemple au problème de Dirichlet homogène si la condition de bord est $p|_{\Gamma} = 0$ ou au problème de Neumann homogène si la condition au bord est $\frac{\partial p}{\partial \nu} = 0$.

On considère des solutions harmoniques de cette équation, c'est-à-dire de la forme

$$ U(x, t) = \mathfrak{Re}\{u(x)e^{-i\omega t}\}, $$

de pulsation $\omega > 0$, correspondant à une période temporelle $T := \frac{2\pi}{\omega}$ et à une fréquence $f := \frac{1}{T} = \frac{\omega}{2\pi}$.

La partie complexe de $U$ dépendante en espace $u$ satisfait l'équation de Helmholtz :

$$ -\Delta u - \kappa^2 u = 0, $$

où le nombre d'onde $\kappa$ est donné par $\kappa := \frac{\omega}{c}$. Cette équation est également complétée de conditions aux limites sur le domaine considéré, de type Dirichlet, Neumann ou Robin.

### 2.1.1 Problème modèle

Soit $\Omega \subset \mathbb{R}^2$ borné, de bord $\Gamma$ suffisamment régulier. Le problème modèle que nous considérons est le suivant : pour $f \in L^2(\Omega)$ et $g \in L^2(\Gamma)$,

$$ \left\{ \begin{array}{ll} \text{Trouver } u \in \mathbb{H}^1(\Omega) \text{ tel que} & \\ -\Delta u - \kappa^2 u = f \quad \text{dans } \Omega, & (1) \\ (\partial_\nu + i\kappa)u = g \quad \text{sur } \Gamma, & \end{array} \right. $$

où l'inconnue $u$ est à valeurs complexes. Dans notre travail, nous prendrons toujours le terme source volumique $f = 0$.

Rappelons le résultat d'existence et d'unicité suivant ([6, Theorem 1 p. 4]):

#### Théorème 1

$Si (f,g) \in L^2(\Omega) \times L^2(\Gamma)$, alors (1) admet une unique solution dans $\mathbb{H}^1(\Omega)$, qui satisfait $\Delta u \in L^2(\Omega)$ et $\partial_\nu u \in L^2(\Gamma)$. De plus il existe une constante $C_\Omega > 0$ telle que

$$ ||u||_{\mathbb{H}^1(\Omega)} + ||\partial_{\nu} u||_{L^2(\Gamma)} \leq C_{\Omega} (||f||_{L^2(\Omega)} + ||g||_{L^2(\Gamma)}). $$

### 2.1.2 Domaine $\Omega$

Notre domaine de référence sera un simple carré :

$$ \Omega = [0; 1] \times [-0.5; 0.5], \qquad (2) $$
que nous matérialiserons sous la forme d'un maillage conforme, constitué de huit triangles dans sa forme la plus simple. Nous présentons en figure 1 la répartition des éléments qui le composent. C'est ce maillage que nous raffinerons progressivement dans nos expériences numériques.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Domaine Ω, maillage de référence.</figcaption></figure>

## 2.2 Méthodes de Trefftz

### 2.2.1 Introduction

On se base ici sur la présentation [14]. Les méthodes de Trefftz sont des schémas numériques discrets dans lesquelles les fonctions de test et d'approximation sont solutions locales (i.e. dans chaque cellule du maillage considéré) de l'EDP étudiée. Elles sont particulièrement utilisées pour les problèmes d'ondes harmoniques, car leurs espaces d'approximation contiennent des fonctions de bases oscillantes et peuvent donc avoir de meilleures propriétés d'approximation que des espaces polynomiaux par morceaux.

Comme les méthodes d'éléments finis, les méthodes de Trefftz réalisent une discrétisation du volume, mais contrairement à celles-ci elles peuvent nécessiter bien moins de degrés de libertés pour atteindre une précision équivalente. Elles incorporent de plus les propriétés de la solution du problème dans les espaces d'approximation. Elles partagent avec les méthodes de volumes finis certaines caractéristiques comme l'intégration sur des variétés de dimension inférieure, mais à leur différence elles ont l'avantage de ne pas nécessiter l'évaluation d'intégrales singulières.

### 2.2.2 Méthodes de Trefftz et équation de Helmholtz

Les solutions de l'équation de Helmholtz oscillant avec une longueur d'onde proportionnelle à l'inverse de $κ$, l'utilisation de polynômes par morceaux n'est pas efficace pour en obtenir une bonne approximation pour des fréquences élevées (i.e. des valeurs de $κ$ de $10^2$ à $10^4$). C'est ici au contraire que les méthodes de Trefftz montrent leur plus grande qualité, du fait des propriétés de leurs fonctions de bases.

Un autre intérêt de ces méthodes est que tout espace de Trefftz discret peut être employé dans toute formulation variationnelle de Trefftz. Il existe de nombreuses méthodes de Trefftz pour l'équation de Helmholtz, et il est difficile de conclure de leur étude comparative qu'une formulation est clairement préférable à une autre. Une conclusion d'ordre général est que le choix de l'espace d'approximation compte plus que la formulation variationnelle.
Les méthodes de trefftz souffrent de deux problèmes principaux : le mauvais conditionnement du système linéaire, dû à une faible indépendance linéaire des fonctions de base, et le besoin de méthodes de quadrature numérique pour les intégrandes oscillantes. Heureusement pour ce dernier, dans le cas de l'espace d'approximation à base d'ondes planes que nous utiliserons, les intégrales sur tout maillage polygonal ou polyhédral peuvent être calculées analytiquement facilement (voir [13] section 4). Nous n'échapperons malheureusement pas aux difficultés liées au premier.

## 2.3 Formulation ultra faible

La formulation variationnelle que nous utiliserons pour notre méthode de Trefftz tout au long de notre stage est la formulation variationnelle ultra-faible (UWVF) développée par Bruno Desprès et Olivier Cessenat, et exposée dans leur article de 1998 [5]. Nous renvoyons à cet article pour une présentation détaillée, et nous contentons ici de rappeler les principaux résultats qui nous seront utiles.

### 2.3.1 Cadre continu

Pour un maillage dont les $K$ éléments sont $\Omega_k$ de bords $\partial\Omega_k$ pour $1 \le k \le K$, le cadre fonctionnel de la formulation ultra faible est le suivant : $V := \prod_{k=1}^K \mathbf{L}^2(\partial\Omega_k)$, avec le produit scalaire associé définissant une norme et une norme induite : $(x, y)_{V} = \sum_k \int_{\partial\Omega_k} x_{|_{\partial\Omega_k}} \bar{y}_{|_{\partial\Omega_k}}$.

Soit $\partial_{\nu_k}$ la dérivée normale sortante à $\partial\Omega_k$. La valeur de l'inconnue $x$ de la formulation ultra-faible est définie de $u$ solution de (1), en supposant l'hypothèse de régularité $x \in V$, comme $x_{|_{\partial\Omega_k}} = ((-\partial_{\nu_k} + i\kappa)u_{|_{\partial\Omega_k}})_{|_{\partial\Omega_k}}$.

Nous reproduisons ici le théorème principal de l'article[5, théorème 4.3 p. 6], définissant la formulation ultra faible.

**Théorème 2**

Soit $u \in \mathbb{H}^1(\Omega)$ satisfaisant l'hypothèse de régularité $\partial_{\nu_k} u \in \mathbb{L}^2(\partial\Omega_k)$ pour tout $k$, et solution du problème de Helmholtz (1).

Alors $x$ définie par $x_{|_{\partial\Omega_k}} = x_k$, avec $x_k = ((-\partial_{\nu_k} + i\kappa)u_{|_{\partial\Omega_k}})_{|_{\partial\Omega_k}}$ satisfait

$$ 
\begin{aligned}
& \sum_k \int_{\partial\Omega_k} x_k \overline{(-\partial\Omega_k + i\kappa)e_k} - \sum_{k,j} \int_{\sum_{k,j}} x_j \overline{(\partial\Omega_k + i\kappa)e_k} \\
& = -2i\kappa \sum_k \int_{\Omega_k} f \overline{e_k} + \sum_k \int_{\Gamma_k} g \overline{(\partial\Omega_k + i\kappa)e_k},
\end{aligned}
\quad (3) $$

pour tout

$$ e \in H, \quad e = (e_k)_{k=1,...,K}, \quad H = \prod_{k=1}^{K} H_k, \quad (4) $$

avec

$$ H_k = \left\{ v_k \in \mathbb{H}^1(\Omega_k) \mid \begin{cases} (-\Delta - \kappa^2)v_k = 0 & \text{in } \Omega_k \\ (-\partial_{\nu_k} + i\kappa)(v_k)_{|_{\partial\Omega_k}} \in \mathbb{L}^2(\partial\Omega_k) \end{cases} \right\} \quad (5) $$

Réciproquement, si $x$ est solution de (3), alors la fonction $u$ définie par

$$ 
\left\{ 
\begin{array}{l}
u|_{\Omega_k} = u_k \\
(-\Delta - \kappa^2)u_k = f|_{\Omega_k} \\
(-\partial_{\nu_k} + i\kappa)u_k = x_k
\end{array} 
\right. 
\quad (6) $$

où $u_k \in \mathbb{H}^1(\Omega_k)$, est l'unique solution de (1).
Afin de pouvoir réécrire ce nouveau problème sous forme d'opérateurs, on introduit les définitions suivantes.

**Définition 3**

On définit les applications d'extension E et $E_f$ par $E_f := \begin{cases} V \to H = \prod_{k=1}^{K} H_k , \\ z \to e = (e_k), e_k = e_{|\Omega_k} \end{cases}$ où $e_k$ est l'unique solution de (par définition de $H$) $\begin{cases} (-\partial_\nu + i\kappa)e_k = z_{|\partial\Omega_k} \text{ sur } \partial\Omega_k \\ (-\Delta - \kappa^2)e_k = f_{|\Omega_k} \text{ dans } \Omega_k \end{cases}$ et $E := E_0$.

$E$ est un opérateur linéaire, tandis que $E_f$, ($f \neq 0$) est un opérateur affine.

**Définition 4**

On définit l'opérateur $F \in \mathcal{L}(V)$, qui à la trace sortante $(-\partial_{\nu_k} + i\kappa)e_k$ associe la trace entrante $(\partial_{\nu_k} + i\omega)e_k$

$$F := \begin{cases} V \to V \\ z \to Fz = ((\partial_{\nu_k} + i\kappa)(E(z)_{\Omega_k})\partial_{\Omega_k})_k. \end{cases} \qquad (7)$$

F est analogue à un opérateur "Neumann - Dirichlet", mais appliqué à des conditions de Robin sur les bords des éléments du maillage. On a $Fz \in \prod_k (\mathbb{L}^2(\partial\Omega_k) = V)$.

**Définition 5**

On définit l'opérateur linéaire $\Pi$ par $\Pi \in \mathcal{L}(V) \begin{cases} \Pi z_{|\Sigma_{kj}} = z_{|\Sigma_{jk}} \\ \Pi z_{|\Gamma_k} = 0 \end{cases}$

L'opérateur $\Pi$ satisfait $||\Pi||_V \leq 1$.

**Définition 6**

Soit $F^* \in \mathcal{L}(V)$ ($V$ identifié à son dual) l'adjoint de $F$. On définit $A \in \mathcal{L}(V)$ par $A := F^*\Pi$.

Nous pouvons à présent reformuler le théorème (3) sous forme d'opérateur.

**Théorème 7**

1. Le problème (3) est équivalent à

$$\begin{cases} \text{Trouver } x \in V \text{ tel que } \forall y \in V, \\ (x, y)_V - (\Pi x, Fy)_V = (b, y)_V \end{cases} \qquad (8)$$

où le terme de droite $b \in V$ est défini, via le théorème de représentation de Riesz, par

$$\forall y \in V, \quad (b, y)_V = -2i\omega \sum_k \int_{\Omega_k} f \overline{E(y)_{\Omega_k}} + \sum_k \int_{\Gamma_k} g \overline{F(y)_{\Gamma_k}} \qquad (9)$$

2. Si $u$ est solution de (1), alors $x = (\partial_{\nu_k} + i\omega)u$ est solution de (9).

3. Réciproquement, si $x$ est solution de (9), alors $u = E_f(x)$ est l'unique solution de (1).

Le problème (9) est équivalent à

$$\begin{cases} \text{Etant donné } b \in V, \text{ trouver } x \in V, \\ (I - A)x = b. \end{cases} \qquad (10)$$

L'opérateur $A = F^*\Pi$ satisfait notamment les propriétés suivantes, particulièrement importantes puisqu'elles garantissent la convergence des méthodes itératives :

— $||A|| \leq 1$,

— L'opérateur $(I - A)$ est injectif.
### 2.3.2 Discrétisation de la formulation variationnelle.

On introduit un espace discret $V_h \subset V$ de dimension finie pour obtenir

$$ \begin{cases} \forall y_h \in V_h, \\ (x_h, y_h)_V - (\Pi s_h, Fy_h)_V = (b, y_h)_V \end{cases} \qquad (11) $$

#### Théorème 8

Le problème (11) a une unique solution.

### 2.3.3 Fonctions de base

Dans chaque élément $\Omega_k$, on introduit un nombre fini de fonctions $e_{kl}$ à support compact dans $\Omega_k$, solutions indépendantes de l'équation de Helmholtz homogène dans $\Omega_k$ :

$$ \begin{cases} (e_{kl})_{|\Omega_l} = 0 & \text{si } k \neq j \\ (-\Delta - \kappa^2)(e_{kl})_{|\Omega_k} = 0 & \text{et } e_{kl} \neq 0. \end{cases} \qquad (12) $$

L'utilisation d'une base à support compact a pour conséquence que la matrice du système linéaire est creuse. L'index $k$ renvoie à l'élément du maillage dans lequel $e_{kl}$ n'est pas nul. L'index $l$, compris entre 1 et $p$, correspond à l'index local de fonctions de la base dans $\Omega_k$. On considère pour simplifier un nombre $p$ constant de fonctions de base par élément.

Les fonctions $z_{kl}$ qui forment la base de $V_h$ sont définies par

$$ z_{kl} = (-\partial_{\nu_k} + i\kappa)e_{kl} $$

$e_{kl}$ étant solution de (12). On a donc

$$ \begin{cases} (z_{kl})_{|\Omega_j} = 0 & \text{si } k \neq j \\ (z_{kl})_{|\Omega_k} = (-\partial_{\nu_k} + i\kappa)e_{kl}. \end{cases} \qquad (13) $$

#### Lemme 9

Les fonctions $\{z_{kl}\}_{1\le l \le p}$ définies par (13) sont linéairement indépendantes dans $V$ si et seulement si $\{e_{kl}\}_{1\le l \le p}$ sont linéairement indépendantes dans $H_k$.

Sous les hypothèses de ce lemme, $\{z_{kl}\}_{k,l}$ forme une base de $V_h$, que l'on définit donc de la façon suivante :

$$ V_h = \text{Vect} \left\{ (z_{kl})_{\substack{1 \le l \le p \\ 1 \le k \le K}} \right\}. $$

Une solution approchée $x_h$ est alors définie par ses $p \times K$ coefficients complexes $x_{kl}$ dans la base $z_{kl}$ comme suit :

$$ (x_h)|_{\partial \Omega_k} = \sum_l x_{kl} z_{kl} = \sum_l x_{kl} (-\partial_{\nu_k} + i\kappa) e_{kl}. $$

### 2.3.4 Construction des opérateurs discrets

Il nous faut calculer les termes correspondant à la formulation (10). La formulation discrétisée (11) de (8) aboutit à (14)

$$ \begin{cases} \text{Find } X \in \mathbb{C}^{p^K}, \\ (D - C)X = b. \end{cases} \qquad (14) $$

La matrice $D$ correspond au produit scalaire dans $V_h$, la matrice $C$ est la matrice correspondant à la forme bilinéaire $(\Pi x_h, Fy_h)$.
— Les coefficients de $D$, définis par $D_{k,j}^{l,m} = (z_{j,m}, z_{k,l})_V$ sont

$$D_{k,j}^{l,m} = \delta_{kj} \int_{\partial\Omega_k} (-\partial_{\nu_k} + i\kappa)e_{jm} \overline{(-\partial_{\nu_k} + i\kappa)e_{kl}}, \quad (15)$$

— Les coefficients de $C$, définis par $C_{k,j}^{l,m} = (\Pi z_{j,m}, F z_{k,l})_V$ sont

$$C_{k,j}^{l,m} = \int_{\Sigma_{kj}} (\partial_{\nu_k} + i\kappa)e_{jm} \overline{(\partial_{\nu_k} + i\kappa)e_{kl}}, \quad (16)$$

— Les coefficients du second membre $b$, défini par $b_{k,l} = (b, z_{kl})_V$ sont

$$b_{k,l} = -2i\kappa \int_{\Omega_k} f \overline{e_{kl}} + \int_{\Gamma_k} g \overline{(\partial_{\nu_k} + i\kappa)e_{kl}}. \quad (17)$$

## 2.4 Espaces d'approximation

Nous décrivons ici dans les grandes lignes les nouveaux espaces d'approximation que nous allons utiliser, qui sont des sous-espaces de dimension finie de $H$ (voir (4)). Nous renvoyons à l'article [13] (ou à [9] pour la version à trois dimensions) pour la présentation détaillée. Leur construction est particulièrement technique, et puisque cela n'est pas l'objet de notre stage, nous ne rentrerons pas dans ces détails, et tâcherons simplement de donner les éléments suffisants à leur manipulation.

L'idée générale est d'enrichir une base d'ondes planes propagatives par un choix judicieux d'ondes planes évanescentes générées selon un procédé précis, ce qui permet d'atteindre dans les simulations un haut degré de précision, de façon numériquement stable malgré le mauvais conditionnement des systèmes à résoudre.

Commençons par définir les différents types d'ondes en jeu.

### 2.4.1 Ondes circulaires

On suppose ici sans perte de généralité que le domaine est le disque de rayon $\rho > 0$ noté $B_\rho := \{\mathbf{x} \in \mathbb{R}^2 \mid \|\mathbf{x}\| < \rho\}$, ce qui permet une analyse modale via séparation de variables. Les ondes circulaires sont les solutions bornées de l'équation de Helmholtz dans ce domaine qui sont séparables en coordonnées polaires. On s'y réfèrera dans ce rapport sous le nom de *Fonctions de Fourier-Bessel*.

Si l'on définit les produits scalaires et norme suivants : pour tous $u, v \in \mathbb{H}^1(B_1)$,

$$(u, v)_{\mathcal{B}} := (u, v)_{L^2(B_1)} + \kappa^{-2} (\nabla u, \nabla v)_{L^2(B_1)}, \quad \|u\|_{\mathcal{B}}^2 := (u, u)_{\mathcal{B}},$$

nous pouvons alors définir les ondes circulaires de la façon suivante : pour tout $p \in \mathbb{Z}$,

$$ \begin{cases} \tilde{b}_p(\mathbf{x} := J_p(\kappa\rho r)e^{ip\theta}, \quad \forall \mathbf{x} = (r, \theta) \in B_1, \\ b_p := \beta_p \tilde{b}_p, \quad \text{où } \beta_p := \|\tilde{b}_p\|_{\mathcal{B}}^{-1}, \end{cases} \quad \text{et } \mathcal{B} := \overline{\operatorname{span}\{b_p\}_{p \in \mathbb{Z}}} \|\cdot\|_{\mathcal{B}} \subset \mathbb{H}^1(B_1) $$

présentées ici sous une forme normalisée. $J_p$ est la fonction de Bessel de première espèce définie, pour $\alpha \in \mathbb{R}$, par

$$J_{\alpha}(z) = \sum_{p=0}^{\infty} \frac{(-1)^p}{p!\Gamma(p + \alpha + 1)} \left(\frac{z}{2}\right)^{2p+\alpha}.$$

Cette solution de l'équation de Bessel est analytique pour $z \in \mathbb{C}$ sauf en $z = 0$ lorsque $\alpha$ n'est pas entier.

Deux résultats fondamentaux sont le [13, Lemme 2.2] :
**Lemme 10**

L'espace $(\mathcal{B}, \| \cdot \|_{\mathcal{B}})$ est un espace de Hilbert, et la famille $\{b_p\}_{p \in \mathbb{Z}}$ est une base de Hilbert,

ainsi que le [13, Lemme 2.3] :

**Lemme 11**

$u \in \mathbb{H}^1(B_\rho)$ satisfait l'équation de Helmholtz si et seulement si $u \in \mathcal{B}$.

Les ondes circulaires sont utilisées dans de nombreuses méthodes de Trefftz. On peut les répartir en plusieurs catégories : en particulier, les ondes circulaires de mode $|p| < \kappa$ sont nommées modes propagatifs, leur énergie étant située dans l'intérieur du domaine, tandis que celles pour lesquelles $|p| \gg \kappa$ sont nommées évanescentes, et leur énergie se concentre sur le bord du domaine.

Lors de nos expériences numériques, nous évaluerons la qualité de nos approximations en mesurant l'erreur commise (le plus souvent l'erreur $\mathbb{H}^1$ relative) comparativement à une solution de référence du problème d'origine, pour laquelle nous prendrons une fonction de Fourier-Bessel.

### 2.4.2 Ondes planes propagatives

Pour tout angle $\varphi \in [0, 2\pi)$, on définit l'onde plane propagative

$$ \mathrm{PW}_{\varphi}(\mathbf{x}) := e^{i\kappa \mathbf{d}(\varphi) \cdot \mathbf{x}}, \forall \mathbf{x} \in \mathbb{R}^2, \quad \text{où } \mathbf{d}(\varphi) := (\cos\varphi, \sin\varphi) \in \mathbb{R}^2. $$

Toute onde plane propagative satisfait l'équation de Helmholtz homogène, puisque $\mathbf{d}(\varphi) \cdot \mathbf{d}(\varphi) = 1$ pour tout angle $\varphi \in [0, 2\pi)$. En 2D, on pourra obtenir des approximations isotropes en utilisant des angles équidistant sur le cercle unité : pour $M \in \mathbb{N}$, l'espace d'approximation est défini par

$$ \Phi_M := \{M^{-1/2} \mathrm{PW}_{\varphi_{M,m}}\}_{m=1}^M, \quad \text{où } \varphi_{M,m} := \frac{2\pi m}{M}, \quad 1 \le m \le M. $$

Un des résultats importants de [13] (voir section 4.3) est le théorème 4.3 :

**Théorème 12**

Il n'existe pas de suite d'ensembles d'approximation constitués d'ondes planes propagatives qui soit une approximation stable pour l'espace des solutions de Helmholtz sur le disque.

### 2.4.3 Ondes planes évanescentes

Les ondes planes évanescentes sont des ondes planes dont les directions de propagation $\mathbf{d} \in \mathbb{C}^2$ sont à valeurs complexes, contrairement aux ondes propagatives dont les directions sont à valeurs réelles. A la différence des ondes propagatives, les ondes évanescentes n'ont pas d'amplitude constante.

Pour $\mathbf{y} := (\varphi, \zeta) \in [0, 2\pi) \times \mathbb{R}$, on définit une direction complexe

$$ \mathbf{d}(\mathbf{y}) := (\cos(\varphi + i\zeta), \sin(\varphi + i\zeta)) \in \mathbb{C}^2, $$

puis l'onde plane évanescente associée

$$ \mathrm{EW}_{\mathbf{y}}(\mathbf{x}) := e^{i\kappa \mathbf{d}(\mathbf{y}) \cdot \mathbf{x}}, \forall \mathbf{x} \in \mathbb{R}^2, \quad \text{où } \mathbf{d}(\mathbf{y}) := (\cos(\varphi + i\zeta), \sin(\varphi + i\zeta)) \in \mathbb{C}^2. $$

Les ondes planes évanescentes sont solutions de l'équation de Helmholtz homogène puisque $\mathbf{d}(\mathbf{y}) = 1$ pour tout $\mathbf{y} \in [0, 2\pi) \times \mathbb{R}$. Une autre expression d'une onde plane évanescente est la suivante :

$$ \mathrm{EW}_{\varphi, \zeta}(\mathbf{x}) = e^{i\kappa(\cosh\zeta)\mathbf{x} \cdot \mathbf{d}(\varphi) - \kappa(\sinh\zeta)\mathbf{x} \cdot \mathbf{d}^{\perp}(\varphi)} \quad \text{où } \mathbf{d}^{\perp}(\varphi) := (-\sin\varphi, \cos\varphi). $$
L'onde oscille, à fréquence $\kappa \cosh\zeta \ge \kappa$ dans la direction propagative $\mathbf{d}(\varphi) := (\cos\varphi, \sin\varphi)$ parallèlement à $\mathfrak{Re}(\mathbf{d}(\mathbf{y}))$, et décroit exponentiellement à vitesse $\kappa \sinh\zeta$ dans la direction orthogonale $\mathbf{d}^{\perp}$, parallèlement à $\mathfrak{Im}(\mathbf{d}(\mathbf{y}))$. En conséquence, le paramètre $\zeta \in \mathbb{R}$, qui contrôle la partie imaginaire de cet angle, est appelé paramètre d'évanescence.

Du fait que leur module décroit exponentiellement dans la direction $\mathfrak{Im}(\mathbf{d})$, les ondes évanescentes sont localisées dans des domaines physiques bornés mais également en Fourier, et sont donc des candidats naturels pour approcher le contenu haute-fréquence en Fourier de certaines solutions de l'équation de Helmholtz.

Un des principaux résultats de [13, section 6] est le Théorème 6.7, qui affirme que toute solution de l'équation de Helmholtz dans le disque unité peut être représentée sous la forme d'une superposition continue d'ondes planes évanescentes, qui a la propriété fondamentale d'être stable (c'est-à-dire de densité bornée dans un certain espace $\mathbb{L}^2$ pondéré).

### 2.4.4 Construction des espaces

Le choix de $\varphi$ et $\zeta$ est non trivial et nous ne le détaillerons pas. Il suit différentes distributions, auxquelles nous nous référerons quelques fois sous les noms `uniform_sampling`, `sobol_sampling`, et `random_sampling`.

La recette de construction des espaces d'approximation est basée sur deux paramètres :

- $N$ le nombre d'ondes par cellule de l'espace de Trefftz ($N_{\Omega_k}$ lorsque celui-ci sera différent pour chaque cellule $\Omega_k$ d'un maillage), qui est aussi le nombre de valeurs générées pour $\varphi$ et $\zeta$,

- $P$ le paramètre de troncation en Fourier, qui varie linéairement par rapport à $N$.

### 2.4.5 Borne supérieure sur $N$

Le problème est le suivant : lorsque $N$ devient grand, la recette de construction implique que $\zeta$ est choisi de plus en plus grand. La partie réelle de l'argument de l'exponentielle doit être bornée pour éviter un *overflow*. On le voit en exprimant le module de la partie réelle d'une onde évanescente : pour toute cellule $\Omega_k$ de rayon du cercle circonscrit $r_k$, on a, en utilisant la seconde expression de l'onde évanescente,

$$ 
\begin{aligned}
|\mathfrak{Re}(EW_{\varphi,\zeta}(\mathbf{x}))| &= |\mathfrak{Re}(e^{i\kappa(\cosh\zeta)\mathbf{x}\cdot\mathbf{d}(\varphi)} - \kappa(\sinh\zeta)\mathbf{x}\cdot\mathbf{d}^{\perp}(\varphi))|, \quad \text{où } \mathbf{d}^{\perp}(\varphi) := (-\sin\varphi, \cos\varphi) \\
&\le e^{\kappa r_k \sinh|\zeta|}, \quad \text{car } |\mathbf{d}^{\perp} \cdot \mathbf{x}| \le r_k.
\end{aligned}
$$

Au cours de nos expériences numériques nous ferons régulièrement face à ce problème d'overflow lorsque nous souhaiterons des valeurs trop importantes de $N$ (> 1000). Notre pratique de ces espaces nous conduira à choisir de ne conserver lors de leur construction que les ondes évanescentes dont le paramètre d'évanescence $\zeta$ satisfait $\kappa e^{|\zeta|} \le 700$.
## 2.5 Opérateur T

Au cours de notre recherche, nous avons à un moment souhaité trouver des conditions de transmission plus favorables que celle de Robin dans la formulation ultra-faible, afin d'améliorer l'efficacité de notre méthode. Nous avons donc étudié l'article [6], section 3.3 notamment, dont nous reproduisons ici les grandes lignes, afin de pouvoir ensuite dériver une formulation faible à partir du problème d'origine en incluant un opérateur T abstrait dans la partie de Dirichlet de la condition de transmission de Robin. Cette piste n'a pas été fructueuse, mais nous la consignons ici avant de passer au cœur de notre travail.

### 2.5.1 Présentation

On souhaite remplacer les quantités $(\partial_{\nu_k} \pm i\kappa)u_{|_{\Sigma_{kj}}}$ par des quantités

$$ (\partial_{\nu_k} \pm i\kappa T_{jk})u_{|_{\Sigma_{kj}}} \qquad (18) $$

où $T_{jk} = T_{kj}$ sont des opérateurs agissant sur des fonctions définies sur $\partial\Omega_k \cap \partial\Omega_j$.

Pour $j, k$ tels que $\partial\Omega_k \cap \partial\Omega_j \neq \emptyset$, on introduit l'opérateur $T_{jk}$ satisfaisant :

$$ T_{jk} = T_{kj} \in \mathcal{L}(\mathbb{H}^{1/2}(\partial\Omega_k \cap \partial\Omega_j), \mathbb{H}^{-1/2}(\partial\Omega_k \cap \partial\Omega_j)) \qquad (19) $$

$$ T \text{ est injective.} \qquad (20) $$

L'injectivité de $T_{jk}$ implique que le problème d'origine (1) est équivalent à une collection de problèmes locaux sur $\Omega_k$ :

$$ \left\{ \begin{array}{ll} -\Delta u_k - \kappa^2 u_k = f & \text{dans } \Omega_k \\ \partial_{\nu_k} u_k + i\kappa T_{jk}(u_k)|_{\Sigma_{kj}} = -\partial_{\nu_j} u_j + i\kappa T_{jk}(u_j)|_{\Sigma_{jk}} & \text{sur } \Sigma_{kj}, k \neq j \\ \partial_{\nu_k} u_k + i\kappa u_k = 0 & \text{sur } \Sigma_{kk}. \end{array} \right. \qquad (21) $$

Afin de s'assurer que chaque problème local (21) est bien posé, on fait les hypothèses suivantes sur $T_{jk}$, se ramenant à sa symétrie et sa positivité :

Il existe $s \in [0, (1/2)]$ (qui peut dépendre de $T_{jk}$) et un isomorphisme $\Lambda_{jk} : \mathbb{H}^s(\partial\Omega_k \cap \partial\Omega_j) \to \mathbb{L}^2(\partial\Omega_k \cap \partial\Omega_j)$ tel que $T_{jk} = \Lambda_{jk}^* \Lambda_{jk}$,

où l'opérateur $\Lambda_{jk}^*$, adjoint de $\Lambda_{jk}$, est un isomorphisme de $\mathbb{L}^2(\partial\Omega_k \cap \partial\Omega_j)$ dans $\mathbb{H}^{-s}(\partial\Omega_k \cap \partial\Omega_j)$.

On en déduit que $T_{jk}$ est un opérateur symétrique positif injectif de $\mathbb{H}^s(\partial\Omega_k \cap \partial\Omega_j)$ dans $\mathbb{H}^{-s}(\partial\Omega_k \cap \partial\Omega_j)$. On a

$$ \forall \varphi \in \mathbb{H}^{1/2}(\Sigma_{kj}), \quad \langle T_{jk} \varphi, \varphi \rangle_{\Sigma_{kj}} = \|\Lambda_{jk} \varphi\|_{L^2(\Sigma_{kj})}^2. \qquad (22) $$

On peut noter que l'opérateur identité correspond à $s = 0$ et $\Lambda_{jk} = \Lambda_{jk}^* = I$, le cas intéressant correspondant à $s = 1/2$.

Dans la suite il sera utile de définir $T_{kk} = I$, on a donc $s = 0$ pour $T_{kk}$.

### 2.5.2 Redéfinition des problèmes locaux

On considère à présent les problèmes locaux suivants :

$$ \left\{ \begin{array}{ll} (-\Delta - \omega^2)u_k = f_k & \text{in } \Omega_k \\ (\partial_{\nu_k} + i\kappa T_{jk})u_{|_{\Sigma_{kj}}} = g_{kj} & \text{on } \Sigma_{kj}. \end{array} \right. \qquad (23) $$
Si $g_{kj} \in \mathbb{H}^{-1/2}(\Sigma_{kj})$, puisque $T_{jk} \in \mathcal{L}(\mathbb{H}^s(\Sigma_{kj}), \mathbb{H}^{-s}(\Sigma_{kj}))$ pour $s \in [0, 1/2]$,
on a $\partial_{\nu_k}(u_k)|_{\Sigma_{kj}} \in \mathbb{H}^{-1/2}(\Sigma_{kj})$, et alors on a la propriété suivante :

$$ \forall \varphi \in \mathbb{H}^{1/2}(\Sigma_{kj}), \quad \langle \partial_{\nu_k} u_k, \varphi \rangle_{\partial \Omega_k} = \sum_j \langle \partial_{\nu_k}(u_k)|_{\Sigma_{kj}}, \varphi|_{\Sigma_{kj}} \rangle_{\Sigma_{kj}}. $$

Le problème local :

$$ \left\{ \begin{array}{l} w_k \text{ satisfait } (-\Delta - \kappa^2)w_k = 0 \\ (\partial_{\nu_k} + i\kappa)w_k = (\Pi x)_k \quad \text{sur } \partial\Omega_k \quad (\text{on rappelle que } x_k = (-\partial_{\nu_k} + i\kappa)u_{|\partial\Omega_k}) \end{array} \right. \qquad (24) $$

est alors équivalent au problème variationnel suivant :

$$ \left\{ \begin{array}{l} \text{Trouver } u_k \in \mathbb{H}^1(\Omega_k) \text{ tel que } \forall v_k \in \mathbb{H}^1(\Omega_k), \\ \displaystyle \int_{\Omega_k} (\nabla u_k \nabla v_k - \kappa^2 u_k v_k) dx + i\kappa \sum_j \langle T_{jk}(u_k)|_{\Sigma_{kj}}, (\bar{v}_k)|_{\Sigma_{kj}} \rangle_{\Sigma_{kj}} = \int_{\Omega_k} f v_k dx + i\kappa \sum_j \langle g_{kj}, (\bar{v}_k)|_{\Sigma_{kj}} \rangle_{\Sigma_{kj}}. \end{array} \right. \qquad (25) $$

### Théorème 13

- Supposons $f \in \mathbb{L}^2(\Omega_k)$ et $g_{kj} \in \mathbb{H}^{-1/2}(\Sigma_{kj})$ pour tout $j$, alors (23) admet une unique solution $u_k \in \mathbb{H}^1(\Omega_k)$.

- Si $g_{kj} \in \mathbb{H}^{-s}(\Sigma_{kj})$, les quantités :

$$ (\Lambda_{kj}^*)^{-1} (\partial_{\nu_k} (u_k)|_{\Sigma_{kj}}) \pm i\kappa \Lambda_{kj} ((u_k)|_{\Sigma_{kj}}) $$

appartiennent à $\mathbb{L}^2(\Sigma_{kj})$, et l'on a l'identité suivante pour chaque $k$ :

$$ \sum_j \left\| (\Lambda_{kj}^*)^{-1} (\partial_{\nu_k} (u_k)|_{\Sigma_{kj}}) + i\kappa \Lambda_{kj} ((u_k)|_{\Sigma_{kj}}) \right\|_{\mathbb{L}^2(\Sigma_{kj})}^2 \\ = \sum_j \left\| (\Lambda_{kj}^*)^{-1} (\partial_{\nu_k} (u_k)|_{\Sigma_{kj}}) - i\kappa \Lambda_{kj} ((u_k)|_{\Sigma_{kj}}) \right\|_{\mathbb{L}^2(\Sigma_{kj})}^2. $$

## 2.6 Dérivation de la formulation avec $\Lambda$.

Notre but est donc ici d'obtenir une version de la formulation ultra-faible incluant l'opérateur $\Lambda$. Effectuant les mêmes calculs que [5, eq.(1.8)] mais avec les opérateurs définis ci-dessus, on obtient l'égalité suivante pour chacune des $j$ interfaces de $\Omega_k$ :

$$ \begin{align} & \int_{\Sigma_{kj}} \left( -(\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa \Lambda_{kj} \right) u \overline{\left( -(\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa \Lambda_{kj} \right)} e_k \nonumber \\ &= \int_{\Sigma_{kj}} \left( (\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa \Lambda_{kj} \right) u \overline{\left( (\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa \Lambda_{kj} \right)} e_k \tag{26} \\ & \quad - 2i\kappa \int_{\Sigma_{kj}} \left( (\Lambda_{kj} u) \overline{((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} e_k)} - ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} u) \overline{(\Lambda_{kj} e_k)} \right). \nonumber \end{align} $$

On utilise [6, eq. 106 p. 15] :

$$ ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} (u_{k|\Sigma_{kj}}), \Lambda_{kj} (u_{k|\Sigma_{kj}}))_{L^2(\Sigma_{kj})} = \langle \partial_{\nu_k} u_k, u_k \rangle_{\Sigma_{kj}}. $$
Cela nous permet d'écrire :

$$
\int_{\Sigma_{kj}} (\Lambda_{kj} u) (\overline{(\Lambda_{kj}^*)^{-1} \partial_{\nu_k} e_k}) = ((\Lambda_{kj}(u|_{\Sigma_{kj}}), (\Lambda_{kj}^*)^{-1} \partial_{\nu_k} (e_{k|_{\Sigma_{kj}}}))_{L^2(\Sigma_{kj})} = \langle u, \partial_{\nu_k} e_k \rangle_{\Sigma_{kj}}
$$

et

$$
\int_{\Sigma_{kj}} ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} u) (\overline{\Lambda_{kj} e_k}) = ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} (u|_{\Sigma_{kj}}), \Lambda_{kj} (e_{k|_{\Sigma_{kj}}}))_{L^2(\Sigma_{kj})} = \langle \partial_{\nu_k} u, e_k \rangle_{\Sigma_{kj}},
$$

puis, en sommant sur les *j* faces de $\partial\Omega_k$, en utilisant [5, eq.1.9 p. 7] (avec $f=0$ dans notre cas) :

$$
\begin{cases}
(-\Delta - \kappa^2) \bar{e}_k = 0 & \text{dans } \Omega_k, \\
(-\Delta - \kappa^2) u = 0 & \text{dans } \Omega,
\end{cases}
$$

et la formule de Green (comme [6, eq.107 p. 16]), on aboutit aux identités suivantes :

$$
\sum_j ((\Lambda_{kj}(u|_{\Sigma_{kj}}), (\Lambda_{kj}^*)^{-1} \partial_{\nu_k} (e_{k|_{\Sigma_{kj}}}))_{L^2(\Sigma_{kj})} = \langle u, \partial_{\nu_k} e_k \rangle \partial_{\Omega_k} = \int_{\Omega_k} \nabla u \overline{\nabla e_k} - \kappa^2 u \overline{e_k}
$$

et

$$
\sum_j ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} (u|_{\Sigma_{kj}}), \Lambda_{kj} (e_{k|_{\Sigma_{kj}}}))_{L^2(\Sigma_{kj})} = \langle \partial_{\nu_k} u, e_k \rangle \partial_{\Omega_k} = \int_{\Omega_k} \nabla u \overline{\nabla e_k} - \kappa^2 u \overline{e_k}.
$$

On en déduit ainsi que :

$$
2i\kappa \sum_j \int_{\Sigma_{kj}} ((\Lambda_{kj} u) (\overline{(\Lambda_{kj}^*)^{-1} \partial_{\nu_k} e_k}) - ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} u) (\overline{\Lambda_{kj} e_k})) = 0.
$$

En sommant sur *j* l'égalité (26), et en utilisant, comme [5, eq.1.14 p. 8],

$$
\begin{align*}
& \sum_{k,j} \int_{\Sigma_{kj}} ((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_{kj}) u \overline{((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_{kj}) e_k} \\
& + \sum_k \int_{\gamma_k} ((\Lambda_k^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_k) u \overline{((\Lambda_k^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_k) e_k} \\
& = \sum_{k,j} \int_{\Sigma_{kj}} ((\Lambda_{jk}^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_{jk}) u \overline{((\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_{kj}) e_k} \\
& + \sum_k \int_{\gamma_k} g \overline{((\Lambda_k^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_k) e_k},
\end{align*}
$$

on obtient finalement une nouvelle formulation variationnelle ultra-faible similaire à (3), in-
tégrant les opérateurs $\Lambda$ :

$$
\begin{align}
& \sum_k \int_{\Sigma_k} \left( -(\Lambda_k^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_k \right) u \overline{\left( -(\Lambda_k^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_k \right) e_k} \nonumber \\
& - \sum_{k,j} \int_{\Sigma_{kj}} \left( (\Lambda_{jk}^*)^{-1} \partial_{\nu_j} + i\kappa\Lambda_{jk} \right) u \overline{\left( (\Lambda_{kj}^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_{kj} \right) e_k} \tag{27} \\
& = \sum_k \int_{\gamma_k} g \overline{\left( (\Lambda_k^*)^{-1} \partial_{\nu_k} + i\kappa\Lambda_k \right) e_k}. \nonumber
\end{align}
$$

Nous avons passé sans succès du temps à essayer de définir explicitement $\Lambda$ de sorte que
la convergence de nos solveurs itératifs soit accélérée. Notre premier espoir était le choix pour
$\Lambda$ du module $|d|$ des ondes planes évanescentes, par homogénéité avec $\partial_{\nu} u$ dans la condition
de transmission puisque le gradient d'une onde plane évanescente $u$ vaut $i\kappa du$, mais cela n'a
malheureusement pas fonctionné, pas plus que les autres tentatives que nous avons faites.
# 3 Implémentation du système linéaire

## 3.1 Implémentation des matrices

Nous présentons ici le pseudo-code, dans sa version séquentielle, permettant de générer les deux matrices $D$ et $C$, ainsi que le second membre $b$.

### 3.1.1 Matrice D

On rappelle ici l'expression des coefficients de $D$ : pour tout élément $\Omega_k$,

$$D_{k,j}^{l,m} = \delta_{kj} \int_{\partial \Omega_k} (-\partial_{\nu_k} + i\kappa) e_{jm} \overline{(-\partial_{\nu_k} + i\kappa) e_{kl}}.$$

Du fait que $\delta_{kj} = \begin{cases} 1 & k = j \\ 0 & k \neq j \end{cases}$, on observe ainsi que $D$ est purement locale. Elle est donc bloc-diagonale, composée de $K$ blocs (le nombre d'éléments du maillage), et chacun de ces blocs a pour taille $p \times p$, avec $p$ le nombre d'ondes dans la base locale d'approximation.

Nous pouvons donc réécrire ses coefficients

$$D_k^{l,m} = \int_{\partial \Omega_k} (-\partial_{\nu_k} + i\kappa) e_{km} \overline{(-\partial_{\nu_k} + i\kappa) e_{kl}} \quad (28)$$

$D$ est donc une matrice creuse, il serait absurde de l'implémenter pleine. Concernant le choix de la structure de données pour sa mise en forme, nous décidons donc de construire un vecteur de $K$ matrices denses, en parcourant successivement les bases locales de chacun des éléments d'Ω.

**Algorithm 1 Construction de $D$**

1: Initialisation du vecteur de matrices complexes denses $D_{vec}$

2: **for** Toutes les cellules $\Omega_K$ de $\Omega$ et toutes les bases locales associées $\Phi_K$ **do**

3: Initialisation de la matrice dense locale $D_K$

4: **for** Toutes les faces $F$ de $\Omega_K$ **do**

5: Calcul du vecteur normal $\vec{n}_F$

6: **for** Toutes les ondes planes $e_K m$ dans $\Phi_K$ **do**

7: Calculer $w_m = (-\partial_{\nu_K} + i\kappa)e_{Km}$

8: **for** Toutes les ondes planes $e_l$ dans $\Phi_K$ **do**

9: Calculer $w_l = (-\partial_{\nu_K} + i\kappa)e_{Kl}$

10: Calculer $D_K[l, m] = \int_F w_m e_m \overline{w_l e_l}$

11: **end for**

12: **end for**

13: **end for**

14: $D_{vec} \leftarrow D_K$

15: **end for**

### 3.1.2 Matrice C

Les coefficients de $C$ définis par $C_{k,j}^{l,m} = (\Pi z_{j,m}, F z_{k,l})_V$ sont

$$C_{k,j}^{l,m} = \int_{\Sigma_{kj}} (-\partial_{\nu_k} + i\kappa) e_{jm} \overline{(-\partial_{\nu_k} + i\kappa) e_{kl}} \quad (29)$$
Pour chaque interface de l'intérieur du domaine (i.e. du squelette), cette intégrale évalue le produit d'ondes planes des deux cellules partageant cette interface. Pour la construction de C, nous commençons par créer une structure de données pour le squelette, qui fournit pour chaque interface, outre l'arête commune, les deux cellules voisines. Nous itérons ensuite sur ce squelette, et calculons les coefficients en considérant le point du vue tour à tour depuis l'une des cellules puis de l'autre (en inversant le vecteur normal et conjuguant ce qui doit l'être), ce qui symétrise les rôles des ondes dans la boucle intérieure et nous permet d'optimiser le nombre de calculs, en n'évaluant l'intégrale qu'une fois sur deux.

Une question apparue ici a été de sélectionner la structure de données la plus adéquate pour C. En effet, cette matrice étant creuse par blocs, chaque bloc étant dense, nous avions plusieurs possibilités pour la définir.

Dans la première version, nous avons utilisé le format de matrice creuse compressed sparse column (Julia est un langage column major), construit par la définition des trois vecteurs I, J et V du format COO (pour coordinate format), stockant respectivement, uniquement pour les valeurs non nulles de C, les indices de ligne, de colonne, et les valeurs correspondantes, i.e. $C[I[i], J[i]] = V[i]$ si $V[i]$ est non nul, pour i compris entre 1 et le nombre total de coefficients non nuls de C. Nous devons donc gérer la correspondance entre les indices locaux des ondes des bases locales et leur indice global correspondant. L'objet final est donc une matrice creuse globale, de la taille des dimensions de l'espace de Trefftz.

<table><caption>**Algorithm 2** Construction de C creuse</caption><tr><td>1:</td><td>Initialisation des vecteurs COO I, J, V</td></tr><tr><td>2:</td><td>for Toutes les interfaces Σ dans le squelette d'Ω do</td></tr><tr><td>3:</td><td>Identifier les cellules J et K de part et d'autre de Σ</td></tr><tr><td>4:</td><td>for Toutes les ondes planes e<sub>m</sub> dans la base locale Φ<sub>J</sub> do</td></tr><tr><td>5:</td><td>for Toutes les ondes planes e<sub>l</sub> dans la base locale Φ<sub>K</sub> do</td></tr><tr><td>6:</td><td>Calculer In = ∫<sub>Σ</sub> e<sub>m</sub>e<sub>l</sub></td></tr><tr><td>7:</td><td>Calculer le vecteur normal n&#773;<sub>F,K</sub></td></tr><tr><td>8:</td><td>Calculer w<sub>m</sub> = (∂<sub>ν<sub>K</sub></sub> + iκ)e<sub>Jm</sub></td></tr><tr><td>9:</td><td>Calculer w<sub>l</sub> = (∂<sub>ν<sub>K</sub></sub> + iκ)e<sub>Kl</sub></td></tr><tr><td>10:</td><td>I ← indice global de l</td></tr><tr><td>11:</td><td>J ← indice global de m</td></tr><tr><td>12:</td><td>V ← w<sub>m</sub> · <span style="text-decoration: overline;">w&#773;<sub>l</sub></span> × In</td></tr><tr><td>13:</td><td>n&#773;<sub>F,J</sub> = -n&#773;<sub>F,K</sub></td></tr><tr><td>14:</td><td>Calculer w<sub>m</sub> = (∂<sub>ν<sub>J</sub></sub> + iκ)e<sub>Jm</sub></td></tr><tr><td>15:</td><td>Calculer w<sub>l</sub> = (∂<sub>ν<sub>J</sub></sub> + iκ)e<sub>Kl</sub></td></tr><tr><td>16:</td><td>I ← indice global de m</td></tr><tr><td>17:</td><td>J ← indice global de l</td></tr><tr><td>18:</td><td>V ← w<sub>l</sub> · <span style="text-decoration: overline;">w&#773;<sub>m</sub></span> × In</td></tr><tr><td>19:</td><td>end for</td></tr><tr><td>20:</td><td>end for</td></tr><tr><td>21:</td><td>end for</td></tr><tr><td>22:</td><td>return sparse(I,J,V)</td></tr></table>

Nous avons par la suite réalisé que ce format CSC scalaire n'était pas très adapté étant donnée la structure par blocs de notre matrice, et avons alors souhaité trouver une structure de données plus proche du problème : le bloc de C situé à la i-ème "ligne" et à la j-ème "colonne" (au sens où ce bloc est considéré comme un coefficient) correspond à l'échange à l'interface entre la i-ème cellule et la j-ème cellule, du point de vue de la i-ème cellule. Il nous paraissait donc logique
(et utile pour les versions parallélisées ultérieures) de regrouper les blocs communs aux échange entre une même cellule et ses voisines.

Nous avons construit plusieurs types de structures en ce sens, avant de nous fixer, au moment de la version multithread de notre code, sur la structure COO-block suivante : basée sur les informations du squelette du maillage, elle reprend le principe du format COO avec ses trois vecteurs I, J et V, mais les éléments du vecteur de valeurs V sont cette fois les blocs denses de C : V[n] = C[i,j pour I[n]= i et J[n]= j, avec n entre 1 et le nombre total de blocs denses dans C. Notons également que nous avons introduit ici la possibilité que les espaces de test et d'approximation soient différent.

<caption>**Algorithm 3 Construction de C creuse par blocs**</caption>

1: I, J, F ← skeleton_Idx(Ω)
2: Initialiser C, vecteur de matrices denses, chacune de taille les dimensions des bases locales
3: **for** Tous les indices i de I **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;K ← I[i]
5: &nbsp;&nbsp;&nbsp;&nbsp;J ← J[i]
6: &nbsp;&nbsp;&nbsp;&nbsp;F ← F[i]-ème face de la K-ème cellule
7: &nbsp;&nbsp;&nbsp;&nbsp;Calculer le vecteur normal $\vec{n}_{F,K}$
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** Toutes les ondes planes $e_m$ dans la base locale $\Phi_J$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$w_m$ ← $(\partial_{\nu_K} + i\kappa)e_{Jm}$
10: &nbsp;&nbsp;&nbsp;&nbsp;**for** Toutes les ondes planes $e_l$ dans la base locale $\Phi_K$ **do**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$w_l$ ← $(\partial_{\nu_K} + i\kappa)e_{Kl}$
12: &nbsp;&nbsp;&nbsp;&nbsp;C[i][l,m] ← $\int_F w_m e_m \overline{w_l e_l}$
13: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
14: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
15: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
16: **return** I, J, C

3.1.3 Vecteur second membre b

Les coefficients du second membre *b* sont définis par $b_{k,l} = \int_{\Gamma_k} g(\overline{\partial_{\nu_k} + i\kappa})e_{kl}$ où *k* représente l'indice de la cellule et *l* celui de l'onde plane dans la base locale de la cellule *k*.

La construction de *b* ne pose pas de difficulté : on crée deux vecteurs d'indices, où pour chaque face du bord du domaine le premier vecteur contient l'indice de cellule et l'autre contient l'indice de face correspondante. On parcourt ensuite ces deux vecteurs en faisant appel aux fonctions de bases locales adéquates pour évaluer les intégrales. Il faut simplement ici aussi veiller à gérer la correspondance entre indices locaux des ondes et indices globaux de *b* dans lesquels écrire les valeurs.

<caption>**Algorithm 4 Construction de b**</caption>

1: Initialisation à zéro de b, vecteur complexe de la dimension de HΦ
2: **for** Toutes les faces F sur le bord de Ω **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Identifier la cellule K correspondant à F
4: &nbsp;&nbsp;&nbsp;&nbsp;Calculer le vecteur normal $\vec{n}_{F,K}$
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** Toutes les ondes planes $e_l$ dans la base locale $\Phi_K$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$w_l$ = $(\partial_{\nu_K} + i\kappa)e_{Kl}$
7: &nbsp;&nbsp;&nbsp;&nbsp;$$b[\text{indice global de } l] = \overline{\int_F g \overline{w_l e_l}}$$
8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
9: **end for**
## 3.2 Vérification de la structure des matrices

On vérifie ici graphiquement les structures des matrices assemblées (voir la figure 2 pour les matrices construites à partir du maillage 1 à huit cellules) : pour $K$ le nombre de cellules du maillage et $p$ le nombre d'ondes planes par élément (on se restreint au cas où les espaces de test et d'approximation sont identiques), elles sont toutes de taille $pK \times pK$, et ont toutes une structure par blocs, avec $K$ lignes et $K$ colonnes de blocs, chacun étant de taille $p \times p$.

- La matrice $D$ est purement locale, donc bloc-diagonale,

- La matrice $C$ est un opérateur d'échange : les blocs non nuls sur la $k$-ème ligne de bloc (correspondant à la $k$-ème cellule du maillage) sont situés sur les colonnes dont les indices correspondent à ceux des cellules voisines de cette $k$-ème cellule.

- La matrice $D^{-1}$ étant bloc-diagonale comme $D$, la matrice $M = D^{-1}C$ conserve naturellement la même structure que $C$.

(voir [6] p.10).

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Structure des matrices du système linéaire, maillage à huit cellules. A gauche $D$, au centre $C$, à droite $I - D^{-1}C$.</figcaption></figure>
# 4 Etude spectrale

Motivons cette étude par le résultat très général suivant sur les méthodes itératives. Soit $A \in \mathbb{C}^{n \times n}$ une matrice inversible, décomposée en $A = M - N$ avec $M$ inversible. Pour la résolution du système linéaire $Ax = b$ avec $b \in \mathbb{C}^n$, on introduit la méthode itérative basée sur ce splitting

$$ \begin{cases} x_0 \in \mathbb{R}^n, \\ M x_{k+1} = N x_k + b, \end{cases} $$

Si la suite $(x_k)_{k \in \mathbb{N}}$ converge vers une limite $x_*$, alors $M x_* = N x_* + b$, de sorte que

$$ x_{k+1} - x_* = M^{-1}N(x_k - x_*), $$

la convergence de la suite $(x_k)_{k \in \mathbb{N}}$ dépend donc des propriétés spectrales de la matrice $M^{-1}N$.

On a le résultat fondamental de convergence suivant (voir par exemple [3, Lemme 13.1.26]) :

## Lemme 14

La méthode itérative définie ci-dessus converge pour tout $x_0 \in \mathbb{C}^n$ si et seulement si le rayon spectral $\rho(M^{-1}N)$ de la matrice d'itération $M^{-1}N$ vérifie $\rho(M^{-1}N) < 1$.

Dans le contexte de notre système linéaire, plutôt que $(D - C)X = b$, nous souhaiterons résoudre le système $(I - D^{-1}C)X = D^{-1}b$. Ce qui précède motive donc l'étude dans cette section des valeurs propres de $D$, puis celles de $M := D^{-1}C$. Cela nous conduit à introduire la SVD régularisée de $D$.

Notre étude va dans le même sens que celle réalisée dans [4], mais les espaces d'approximation utilisés ne sont pas les mêmes. En conséquence, l'efficacité de l'espace de Trefftz que nous utilisons nous permet de pousser plus loin les mesures, en taille des espaces locaux comme en nombre d'onde $\kappa$ du problème de Helmholtz considéré.

## 4.1 Propriétés de D

Citons [5, lemme 2.4 p. 16] :

### Lemme 15

Soit $D$ la matrice construite de fonctions $e_{kl}$ solutions d'un problème de Helmholtz homogène dans chaque élément $\Omega_k$ et linéairement indépendantes sur l'élément considéré. Alors $D$ est hermitienne définie positive.

Si les fonctions de la base dérivent d'ondes planes de directions $\vec{v}_{kl}$ dans l'élément $\Omega_k$, le certificate défini de $D$ est équivalent au fait que les directions des ondes planes $e_{kl}$ soient toutes distinctes (sur chaque élément), ou formellement, $\forall k, \forall (i, j), \vec{v}_{ki} \neq \vec{v}_{kj}$.

La base d'approximation de l'espace de Trefftz que nous utilisons vérifient les hypothèses de ce lemme. En effet, d'après [13, section 5.1 p. 12], les ondes planes évanescentes sont solutions de l'équation de Helmholtz homogène.

De plus d'après [13, lemmes 6.2 et 6.5 p. 14-16], il s'agit d'une base hilbertienne, ses éléments sont donc linéairement indépendants.

On vérifie que notre matrice $D$ construite par bloc est bien hermitienne, en l'assemblant sous forme d'une matrice creuse : l'appel à une fonction de test ishermitian renvoie bien true à tous nos essais, quels que soient les paramètres du problème et de l'espace de Trefftz construit.
Les tests nous montrent par contre que $D$ n'est pas définie positive dès que le nombre d'ondes planes par élément excède un certain seuil, en fonction de $k$ semble-t-il. Nous n'avons pas réellement compris sous quelles conditions précises du problème ce phénomène apparaissait (nous avons d'abord cru que c'était sous la condition que le nombre d'ondes planes soit au moins trois fois supérieur au nombre d'onde $k$, avant de reproduire cette situation sans que cette hypothèse soit vérifiée), mais penchons pour une instabilité numérique.

En effet, nous avons, pour différents paramètres du problème et de l'espace de Trefftz, affiché les valeurs propres négatives de $D$, et quoi que leur proportion était croissante par rapport à la taille globale de l'espace, celles-ci étaient toujours très proches de zéro. Le maillage est ici le carré unité à huit cellules, et les résultats sont similaires quelle que soit la distribution choisie pour l'espace d'approximation.

<table><thead><tr><th>k</th><th>Ondes/cellule</th><th>dim HΦ</th><th>#v.p.nég.</th><th>max(abs(v.p. nég.))</th></tr></thead><tbody><tr><td>5</td><td>10</td><td>80</td><td>0</td><td>/</td></tr><tr><td>5</td><td>25</td><td>200</td><td>8</td><td>2.5297e-14</td></tr><tr><td>5</td><td>100</td><td>800</td><td>245</td><td>2.3435e-13</td></tr><tr><td>50</td><td>300</td><td>2400</td><td>816</td><td>1.5772e-11</td></tr><tr><td>100</td><td>500</td><td>4000</td><td>1354</td><td>8.0978e-11</td></tr></tbody></table>

## 4.2 Valeurs propres de $D^{-1}C$

Plutôt que $(D - C)X = b$, nous souhaiterons résoudre le système $(I - D^{-1}C)X = D^{-1}b$. Il est donc naturel de s'intéresser aux propriétés de $D^{-1}C$. Citons [5, lemme 2.5 p. 17] :

**Lemme 16**

Les valeurs propres de $M = D^{-1}C$ sont situées dans le disque complexe unité, la valeur 1 étant exclue.

Tâchons de vérifier ce lemme dans le contexte de notre espace d'approximation. On utilise pour cela la version creuse de $C$, et l'on effectue une factorisation LU de $D$ pour construire $M$.

### 4.2.1 Affichage des valeurs propres

On écrit une fonction qui calcule et affiche dans le plan complexe les valeurs propres d'une matrice passée en paramètre, ainsi que son rayon spectral $\rho$ et le cercle unité.

### 4.2.2 tests

On se situe sur le maillage de référence 1 à huit cellules. On fixe le nombre d'onde de l'équation $k$ et l'on fait varier le nombre d'ondes planes par cellule. On procède au même test pour les différentes distributions disponibles de génération des paramètres d'ondes évanescentes $\varphi$ et $\zeta$. Les résultats de quelques-uns de nos tests sont montrés dans les figures 3 et 4.

Quels que soient les paramètres de nos tests, on observe qu'à partir d'un certain seuil d'ondes planes par cellule, les valeurs propres sortent du cercle unité, et le rayon spectral de $M$ ne cesse de croître à mesure que l'on augmente la dimension de l'espace d'approximation.

Ce phénomène contredit la théorie citée en début de partie, et motive l'introduction de la SVD régularisée comme solution pour y remédier.
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Valeurs propres de $D^{-1}C$, $\kappa = 5$, distribution de type uniform sampled.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Valeurs propres de $D^{-1}C$, distribution de type random sampled.
En haut, $\kappa = 5$ avec 20, 30 et 50 ondes. En bas, $\kappa = 50$ avec 40, 60 et 70 ondes.</figcaption></figure>
## 4.3 SVD régularisée

Une conséquence de nos observations précédentes est que $D^{-1}C$ et $I - D^{-1}C$, sont mal conditionnées, et leur conditionnement empire à mesure que l'on augmente la dimension de l'espace d'approximation. On présente l'évolution du conditionnement de ces deux matrices en fonction du nombre d'ondes par cellule à la figure 5.

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

Afin de résoudre ce problème, on s'appuie sur la méthode présentée dans [13, section 3.3 p. 7], mais en l'appliquant à $D$ :

— On effectue la décomposition en valeurs singulières de $D$ :

$$D = U\Sigma V^*,$$

dans laquelle les valeurs singulières $\sigma_{max} := \sigma_1 \ge \sigma_2 \ge \dots$ sont ordonnées en ordre décroissant.

— La régularisation consiste à remplacer par 0 les valeurs singulières inférieures à un seuil $\epsilon_{svd} \in (0, 1]$ défini par l'utilisateur (et que l'on prendra le plus souvent égal à $1e-10$), ce qui conduit à la création de la matrice diagonale $\Sigma_{\epsilon}$ correspondante, puis à l'approximation

$$D_{\epsilon} := U\Sigma_{\epsilon}V^*.$$

— On peut alors calculer $D_{\epsilon}^{-1}C$ de la façon suivante :

$$D_{\epsilon}^{-1}C = V\Sigma_{\epsilon}^{\dagger}U^{*}C, \qquad (30)$$

où $\Sigma_{\epsilon}^{\dagger}$ est la pseudo-inverse de $\Sigma_{\epsilon}$, c'est-à-dire la matrice diagonale, de coefficients

$$(\Sigma_{\epsilon}^{\dagger})_{j,j} = \begin{cases} (\Sigma_{j,j})^{-1} & \text{si } \Sigma_{j,j} \ge \epsilon_{svd}\sigma_{max} \\ 0 & \text{sinon} \end{cases}.$$

Relevons ici la remarque de fin de [13, section 3.3] : afin de calculer 30 correctement d'un point de vue numérique, il faut effectuer les produits de droite à gauche, c'est-à-dire $D^{-1}C = V(\Sigma_{\epsilon}^{\dagger}(U^{*}C))$. Nous avons expérimenté à nos dépens de ne pas suivre cette indication lors de notre première implémentation, ce qui avait pour conséquence de limiter la précision de l'approximation plusieurs ordres de grandeurs en-dessous de celle dont nos algorithmes étaient capables.
## 4.4 Valeurs propres de $D_{\epsilon}^{-1}C$

### 4.4.1 Implémentation de $D_{\epsilon}^{-1}$

Nous écrivons une fonction qui calcule $D_{\epsilon}^{-1}$, la SVD pseudo-inverse régularisée de $D$. La structure bloc-diagonale de $D$ permet de le faire en ne calculant que les SVD régularisées $\{(D_k)_{\epsilon}^{-1}\}_k$ de ses blocs denses $\{D_k\}_k$. Nous pouvons ensuite tirer profit des vecteurs d'indices de blocs I et J issus de la construction de $C$ pour procéder aux produits matriciels adéquats des blocs de $D_{\epsilon}^{-1}$ et $C$ et assembler le membre de gauche $I - D_{\epsilon}^{-1}C$ du système linéaire. La construction de $b' = D_{\epsilon}^{-1}b$ ne pose pas de difficulté par le même procédé.

### 4.4.2 Résultats

Nous réitérons alors nos tests d'affichage de valeurs propres, mais cette fois pour $D_{\epsilon}^{-1}C$. Quelques résultats sont présentés page suivante à la figure 7.

Le résultat est similaire quels que soient les choix de distributions pour $\varphi$ et $\zeta$, du nombre d'onde $k$, du nombre d'ondes planes par cellule, et de $\epsilon_{svd}$ : les valeurs propres restent à présent strictement inscrites dans le cercle unité, le lemme [5, lemme 2.5] est donc maintenant vérifié dans notre contexte. Nous présentons également sur deux exemples les comparaisons entre valeurs propres de la version régularisée et de la version non régularisée en figure 6.

<figure><img src="image_7.png" /><figcaption>(a) $k = 5$, distribution : uniformSampling</figcaption></figure>

<caption>(b) $k = 50$, distribution : sobolSampling</caption>

<caption>FIGURE 6 – Comparaisons des valeurs propres de $D^{-1}C$ et $D_{\epsilon}^{-1}C$.</caption>

On observe également sur les graphiques de la version régularisée les deux phénomènes suivants :
<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Valeurs propres de $D_{\epsilon}^{-1}C$ en fonction du nombre d'ondes par cellule, $\epsilon_{svd} = 1e-10$.</figcaption></figure>
— Plus le nombre d'ondes par élément est élevé, plus il y a de valeurs propres de module proche de 1,

— Il semble apparaître deux valeurs d'adhérence, diamétralement situées sur le cercle unité et se rapprochant de la droite réelle à mesure que l'on accroît le nombre d'ondes par élément.

On affiche à présent les valeurs propres de $I - D_{\epsilon}^{-1}C$ pour des valeurs croissantes d'ondes par éléments, et le conditionnement des matrices correspondantes (voir figure 8). Le conditionnement reste médiocre, ce qui nous incitera à chercher des méthodes de préconditionnement pour notre système, mais au moins à présent il semble ne plus croître indéfiniment avec le nombre d'ondes par cellule. On se doute que les valeurs propres s'accumulant à proximité du cercle unité, et particulièrement celles près de l'origine, ralentiront la convergence de nos algorithmes itératifs.

<figure><img src="image_9.png" /><figcaption>FIGURE 8 – Valeurs propres et conditionnement de $I - D_{\epsilon}^{-1}C$.</figcaption></figure>
# 5 Algorithmes séquentiels de résolution

On présente dans cette section les trois méthodes de résolution du système linéaire issu de la formulation ultra-faible

$$ (I - M)X = b' \quad \text{avec} \quad \begin{cases} M = D^{-1}C \\ b' = D^{-1}b \end{cases}. \qquad (31) $$

## 5.1 Résolution directe

La matrice du membre de gauche du système n'étant pas symétrique définie positive, nous ne pouvons pas utiliser une factorisation de Cholesky. De plus, la SVD n'est malheureusement pas disponible pour les matrices creuses exprimées sous la forme du type sparse du package SparseArrays.jl, telle que notre matrice. Nous sommes donc contraints d'utiliser une décomposition LU pour la résolution directe de notre système.

Cette résolution directe nous servira de référence, notamment pour la comparaison des erreurs en norme $H^1$.

## 5.2 Algorithme de Richardson

### 5.2.1 Présentation

La première méthode itérative que nous implémentons est celle, présentée sous deux variantes que nous reproduisons ici, dans [5], p.18 :

— i) Soit $\beta \in ]0.5; 1[$. L'algorithme est

$$ \begin{cases} X_1 = \beta b', \\ X_{n+1} = \beta b' + [(1 - \beta)I + \beta M]X_n. \end{cases} \qquad (32) $$

— Soit $\beta_n \in ]0.5; 1]$; $-\epsilon[, \epsilon > 0$ une suite de réels aléatoires de distribution uniforme. L'algorithme est

$$ \beta_n \in ]0.5; 1 - \epsilon[, \epsilon > 0 \begin{cases} X_1 = \beta_1 b', \\ X_{n+1} = \beta_n b' + [(1 - \beta_n)I + \beta_n M]X_n. \end{cases} \qquad (33) $$

On le réécrit sous la forme

$$ \begin{aligned} X_{n+1} &= \beta b' + [I - \beta(I - M)]X_n \\ &= X_n + \beta r_n, \quad \text{avec } r_n := b' - (I - M)X_n, \end{aligned} $$

et l'on obtient l'algorithme de Richardson :

— $X_0 \in \mathbb{C}^n$,

— $\begin{cases} r_n = b' - (I - M)X_n, \\ X_{n+1} = \beta r_n. \end{cases}$

### 5.2.2 Implémentation

On introduit une double condition d'arrêt à l'algorithme : d'une part un nombre maximum d'itérations, et d'autre part un seuil pour l'erreur relative du résidu $\frac{\|r_n\|}{\|b\|} < \epsilon$, avec $\epsilon$ un paramètre utilisateur.

Nous présentons ci-dessous le pseudo-code de cet algorithme, pour les deux versions décrites ci-dessus, pour la résolution d'un système linéaire général $AX = b$.
**Algorithm 5 Richardson**

1: $X_n \leftarrow \beta b$

2: $r_n \leftarrow b - AX_n$

3: $n \leftarrow 0$

4: **while** $r_n > \epsilon_{solver}$ and $n < n_{max}$ **do**

5:     $r_n \leftarrow b - AX_n$

6: **end while**

**Algorithm 6 Richardson avec $\beta$ aléatoire**

1: $\beta \leftarrow 0.5 + 0.5 * \text{rand()}$

2: $X_n \leftarrow \beta b$

3: $r_n \leftarrow b - AX_n$

4: $n \leftarrow 0$

5: **while** $r_n > \epsilon_{solver}$ and $n < n_{max}$ **do**

6:     $\beta \leftarrow 0.5 + 0.5 * \text{rand()}$

7:     $r_n \leftarrow b - AX_n$

8: **end while**

### 5.2.3 Propriétés

On vérifie ici la première partie de [5, remarque 8], la seconde sera testée dans la section suivante : "l'introduction du facteur correctif $\beta \in ]0; 1[$ tend à ramener globalement les valeurs propres de $D^{-1}C$ vers l'origine. Cela suggère une amélioration de la convergence d'un point de vue numérique".

<figure><img src="image_10.png" /><figcaption>FIGURE 9 - Valeurs propres de $D^{-1}C$ et de $\beta D^{-1}C$ pour $\beta=0.6$, pour Richardson.</figcaption></figure>

## 5.3 GMRES

### 5.3.1 Présentation

On considère le système linéaire $Ax = b$ avec $x, b \in \mathbb{R}^n$ et une matrice $A \in \mathbb{R}^{n \times n}$ non singulière.

Pour $r_0 \in \mathbb{R}^n$, on définit l'espace de Krylov $K(k; A; r_0) := \text{Vect}(r_0, Ar_0, \dots, A^{k-1}r_0)$. L'algorithme itératif Gmres (pour Generalized minimal residual) de résolution du système linéaire
ci-dessus construit une base orthonormée $V_k = [v_1, \dots, v_k]$ pour $K(k; A; r_0)$ par le procédé d'Arnoldi (analogue à celui de Gram-Schmidt, en général modifié pour une plus grande stabilité numérique). Dans l'algorithme Gmres complet, la dimension de l'espace de Krylov peut atteindre $n$, mais l'algorithme prend toujours fin en au plus $n$ itérations, en précision infinie tout du moins. Afin de réduire la quantité de mémoire nécessaire au stockage de $V_k$ dans le cas de systèmes de très grandes dimensions, une version restarted de Gmres a été développée, qui restreint la dimension de cet espace à un paramètre $m$, et redémarre toutes les $m$ itérations le procédé d'Arnoldi en utilisant la dernière itération $x_m$ comme vecteur initial. On présente ci-dessous l'algorithme de cette version restart, notée GMRES (m).

<caption>**Algorithm 7 GMRES(m)**</caption>

1: $\epsilon$ est la tolérance pour la norme du résidu.
2: On choisit $x_0$ (souvent nul ou égal à $b$).
3: convergence = false
4: $n \leftarrow 0$
5: **while** convergence==false **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;$r_0 = b - Ax_0$
7: &nbsp;&nbsp;&nbsp;&nbsp;$\beta = ||r_0||$
8: &nbsp;&nbsp;&nbsp;&nbsp;$v_1 = r_0/\beta$
9: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, \dots, m$ **do**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p = Av_j$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1, \dots, j$ **do**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$h_{ij} = v_i^T p$
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$p = p - h_{ij}v_i$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$h_{j+1,j} = ||p||$
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$v_{j+1} = p/h_{j+1,j}$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;calculer $s = ||b - Ax_j||$
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $s < \epsilon$ **then**
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;résoudre $\min_{y_j \in \mathbb{R}^m} ||\beta e_1 - \bar{H}_j y_j||$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_j = x_0 + V_j y_j$
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;convergence == true
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;exit
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;résoudre $\min_{y_m \in \mathbb{R}^m} ||\beta e_1 - \bar{H}_m y_m||$
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_m = x_0 + V_m y_m$
27: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $||b - Ax_m|| < \epsilon$ **then**
28: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;convergence == true
29: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**
30: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x_0 = x_m$
31: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
32: **end while**

Pour $1 \le k \le m$, la matrice $\bar{H}_k = (h_{ij})$ est une matrice de Hessenberg supérieure d'ordre $(k+1) \times k$, et l'on a la relation fondamentale $AV_k = V_{k+1}\bar{H}_k$. Gmres calcule $x = x_0 + V_k y_k$, où $y_k$ est solution du problème aux moindres carrés $\min_{y_m \in \mathbb{R}^m} ||\beta e_1 - \bar{H}_m y_m||$, d'où son nom. Généralement, une factorisation QR de $\bar{H}_k$ (générée par procédé de Householder ou par rotations de Givens) est utilisée pour résoudre ce problème aux moindres carrés. Nous renvoyons à [15]
pour une présentation très détaillée de l'algorithmme.

5.3.2 Version avec matrice assemblée

On utilise dans notre projet l'implémentation en Julia de cet algorithme issue du package
IterativeSolvers.jl. Dans son utilisation la plus élémentaire, il nous suffit de fournir au
solveur la matrice du membre de gauche et le vecteur second membre du système linéaire, dont
les assemblages ont été décrits plus haut. Nous avons à notre disposition les paramètres suivants :

— le paramètre *restart*, dont l'usage empirique veut que sa valeur par défaut est de 20, mais dont nous verrons qu'il faudra nous soucier,

— le seuil $\epsilon$ du résidu définissant la condition d'arrêt de l'algorithmme, et que nous laisserons à 1e-10 dans toutes nos simulations.

5.3.3 Version *matrix-free*

Nous allons profiter de la possibilité qu’offre l’implémentation que l’on utilise de ne pas avoir
besoin de construire la matrice du système : en effet l’objet qui doit être passé en argument au
solveur est d’un type opérateur linéaire abstrait, qui peut être aussi bien implémenté sous la
forme d’une matrice que d’une fonction exprimant une application linéaire.

La raison est la suivante : en fonction de la taille et de la densité du système, devoir en
construire la matrice peut être très coûteux, voire impossible, en mémoire comme en performance.
Au contraire, construire une fonction réalisant l’application linéaire correspondante pourra être
beaucoup plus léger. Il sera donc très avantageux, dans notre objectif d’aborder des problèmes
de grande taille, de tirer profit de cette possibilité.

Dans l'algorithmme Gmres, le seul usage de la matrice du système réside dans le produit
matrice-vecteur apparaissant lors de l'étape d'orthonormalisation d'Arnoldi. C'est donc ce pro-
duit matrice-vecteur que nous implémentons, sous la forme d'une fonction prenant deux vecteurs
X et Y en arguments, et retournant le résultat du produit Y = AX. Cette fonction étant ap-
pelée par un solveur itératif, elle le sera possiblement un grand nombre de fois, nous devons
donc avoir à l'esprit quelques principes de performance, et avant tout éviter toute allocation de
mémoire pendant les itérations. Dans notre première implémentation, si nous avions bien pris
soin de construire au préalable tous les vecteurs d'indices nécessaires pour effectuer les produits
"bloc-matrice"-"bloc-vecteur", nous n'avons pas alloué de mémoire aux résultats des produits
intermédiaires, ce qui résultait en un grand nombre d'allocations. Les mesures de performance
nous ont permis de corriger ce défaut au cours de notre travail (voir chapitre suivant).
# 6 Expériences numériques, algorithmes séquentiels

## 6.1 Tests de validité.

Nous commençons par effectuer une série de tests pour valider nos implémentations : la résolution directe nous permettra tout d'abord de nous assurer que les éléments du système sont corrects, puis nous comparerons nos méthodes itératives à celle-ci.

### 6.1.1 Erreur relative en norme $\mathbb{H}^1$

Pour évaluer nos solveurs, nous effectuons des tests pour référence avec des solutions de type Fourier-Bessel, présentées en section 2.4.1. Ceci nous permet de calculer l'erreur $e_h$ de nos approximations $u_h$ en norme $\mathbb{H}^1$ relative, c'est-à-dire

$$e_h = \frac{\|u - u_h\|_{\mathbb{H}^1(\Omega)}}{\|b\|_{\mathbb{H}^1(\Omega)}} = \frac{(\sum_k \|u - u_h\|_{\mathbb{H}^1(\Omega_k)})^{1/2}}{(\sum_k \|b\|_{\mathbb{H}^1(\Omega_k)})^{1/2}}.$$

Afin d'effectuer le calcul de cette erreur, nous nous appuyons sur une première fonction qui calcule un raffinement d'un maillage donné en argument. La précision de maillage dont nous aurons besoin ici dépend directement de la fréquence d'oscillation $k$ de notre solution, suivant les principes fondamentaux de l'échantillonnage. Le paramètre permettant d'influer sur cet précision du maillage est $lc$, la longueur caractéristique d'un élément, et pour le calcul d'erreur $\mathbb{H}^1$ nous choisissons de la définir selon la formule

$$lc = \frac{2\pi}{50k}.$$

Notons que pour les problèmes en régime de haute fréquence que nous essaierons de simuler, la valeur de $lc$ peut donc être très petite, de l'ordre de $10^{-4}$ pour $k = 1000$. Cela signifie que le temps de calcul du maillage raffiné va devenir prohibitif, de même que celui de l'erreur $\mathbb{H}^1$. Nous ferons l'expérience de ce que cette durée de calcul devienne rapidement dominante dans tous nos tests.

Une seconde fonction calcule ensuite l'erreur $\mathbb{H}^1$ relative sur la base du maillage raffiné par la fonction précédente, en créant d'abord les fonctions d'évaluation locales $u - u_h$ sur chaque cellule, puis en calculant les normes adéquates.

### 6.1.2 Visualisation

Pour tous nos tests, nous enregistrons des fichiers (de formats `.vtm` et `.vtk`) qui nous permettent de visualiser dans le logiciel `Paraview` à la fois :

- la partie réelle de notre approximation et sa valeur absolue,

- la partie réelle de la solution ainsi que sa valeur absolue,

- l'erreur ponctuelle entre les deux ainsi que sa valeur absolue,

toutes ces valeurs étant évaluées pour un raffinement du maillage donné, que nous prendrons au moins égal à $\frac{2\pi}{10\kappa}$ pour avoir au moins 10 points par longueur d'onde.

### 6.1.3 Résolution directe

On commence par tester la version directe de résolution du système. Comme expliqué précédemment, étant donné le format dans lequel est implémenté la matrice, le système est résolu via une factorisation LU.
<figure><img src="image_11.png" /><figcaption>(a) Erreurs ponctuelles en valeur absolue sur le maillage, pour N=20, 80, 160 et 320 : de [4e-9 à 2.5e-5] pour N=20 à [6.6e-16 à 2.5e-8] pour N=320</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Erreur $H^1$-relative</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Partie réelle de l'approximation pour $N = 320$</figcaption></figure>

FIGURE 10 – Test de résolution directe (LU) pour $k = 10$.

Pour $k = 10$, le maillage de référence $1$, $\epsilon_{svd} = 1e-10$ et avec une solution de type Fourier-Bessel, on effectue des simulations pour un nombre d'ondes planes par cellule $N \in [20, 80, 160, 320]$. On calcule l'erreur $H^1$ relative décrite ci-dessus, et l'on enregistre les fichiers pour visualisation dans Paraview. Les résultats sont présentés en figure 10. Ils nous confortent dans le fait que notre implémentation est correcte : la visualisation de la partie réelle des approximations est conforme à nos attentes, et les erreurs diminuent avec l'augmentation du nombre d'ondes planes par cellule. On observe déjà, comme ce sera le cas partout, que les erreurs les plus importantes se concentrent aux interfaces entre les cellules, ce qui motivera la recherche d'un opérateur d'échange plus efficace que celui utilisé ici dans la formulation variationnelle ultra-faible.

### 6.1.4 Richardson

Nous testons à présent notre premier solveur itératif, et le comparons à la résolution directe. Nous prenons les mêmes paramètres que dans le test précédent et choisissons un seuil d'arrêt $\epsilon_{solver} = 1e-10$, le paramètre de relaxation $\beta = 0.7$, et l'on affiche en plus l'historique du résidu. Nous présentons les résultats en figure 11.
<figure><img src="image_14.png" /><figcaption>(a) Erreurs ponctuelles en valeur absolue sur le maillage, pour $N=20, 80, 160$ et $320$ : de $[4e-9]$ à $2.5e-5]$ pour $N=20$ à $[1.6e-15]$ à $3.5e-8]$ pour $N=320$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Test de résolution par Richardson pour $k = 10$, 
erreur en valeur absolue sur le maillage, erreur $H^1$ 
relative.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Test de résolution par Richardson, historique de
convergence.</figcaption></figure>

<caption>FIGURE 11 – Test de résolution par Richardson pour $k = 10$.</caption>

La comparaison avec la version directe nous permet de nous assurer que notre implémentation est correcte, mais on soulève deux points à la lecture de l'historique de convergence :

— L’algorithme semble être de plus en plus lent à converger à mesure que la taille u système augmente, et l’on observe que pour seulement 320 ondes planes par élément on a dépassé le nombre maximal d’itérations, alors que l’on souhaite pouvoir effectuer des simulations de taille bien supérieure. Cela préfigure une difficulté à passer à plus grande échelle.

— La théorie voudrait que l'évolution de la norme du résidu soit monotone, or on constate que cela n'est pas le cas : il y a une première phase pendant laquelle le résidu oscille de façon erratique, quoi que globalement décroissante, puis une seconde phase durant laquelle il est bien monotone comme attendu. Nous avons effectué d'autres tests en construisant l'espace de Trefftz avec des ondes planes propagatives uniquement, et non évanescentes. Le résultat est le suivant : outre une courte amorce erratique, la norme du résidu adopte très rapidement, et jusqu'à la fin de l'algorithme, le comportement monotone attendu. Cette courte amorce reste mystérieuse, mais le reste nous éclaire sur le fait que ce comportement est dû à la nature complexe des coefficients de la base.

Nous effectuons également des tests de comparaison entre la version de l'algorithme de Richardson dans laquelle on choisit le paramètre $\beta$, et celle dans laquelle il est choisi au hasard dans l'intervalle $[0.5, 1[$, afin de vérifier la remarque 9 p.18 de [5] affirmant une réduction de 10 à 20% du nombre d'itérations. Nos tests ont montré que dans notre contexte le nombre d'itérations est absolument identique dans les deux versions, nous utiliserons donc la version avec $\beta$ fixe dans la suite.
### 6.1.5 Gmres

<figure><img src="image_17.png" /><figcaption>FIGURE 12 – Gmres, erreurs $H^1$, historiques de convergence du résidu, et erreurs ponctuelles en valeur absolue sur le maillage pour 20, 80, 160 et 320 ondes.</figcaption></figure>

Nous testons à présent les deux versions de gmres, celle avec matrice assemblée et celle avec produit matrice-vecteur, avec les mêmes paramètres que les tests des autres méthodes. Comme on le voit sur la figure 12, les erreurs $H^1$ et en valeur absolue sur le maillage sont similaires et conformes à la méthode de résolution directe, et les deux versions convergent en un nombre similaire d'itérations. La seule différence notable entre les deux est dans le comportement du la norme du résidu sur l'historique de convergence, qui n'est pas monotone sur la fin de l'algorithme pour la version avec produit matrice-vecteur. On explique ceci par le paramètre restart de l'algorithme, dont la valeur par défaut est 20. Les tests de temps d'exécution nous permettront de les départager.

## 6.2 Tests de convergence en N.

Etant confiant dans ces premiers résultats, nous pouvons débuter les premiers tests de convergence et de performance.
### 6.2.1 Test comparatif des trois méthodes.

Pour ce test, on fixe les paramètres suivants :

- $k = 10$ et une solution de type Fourier-Bessel sur le maillage de référence 1,

- Distribution `uniformSampling` (sauf pour les erreurs $H^1$),

- $\epsilon_{svd} = 1e-10$ et $\epsilon_{solver} = 1e-10$,

et l'on teste les trois algorithmes pour un nombre croissant $N$ d'ondes planes par élément, multiple de 40, de 40 à 600.

Nous commençons par comparer les erreurs $H^1$ (voir figure 13), pour les trois méthodes de sampling de $\varphi$ et $\zeta$. Celles-ci sont identiques pour les trois solveurs, sauf pour $N$ grand où Richardson se différencie (voir ci-dessous), nous obtenons jusque là donc la même précision pour chacune d'elles. La convergence semble être de l'ordre de $O(1/N^2)$, quelle que soit la méthode de sampling utilisée pour les espaces d'approximation.

L'amélioration de la précision semble cependant stagner pour les plus grandes valeurs de $N$. Nous voyons une corrélation entre cela et un phénomène lié au mode de construction des espaces d'approximation eux-même : plus $N$ augmente, plus il y a dans les bases locales des ondes dont les directions de propagation sont proches. En contexte numérique d'arithmétique finie, cela signifie que l'on perd l'indépendance linaire des fonctions de base, entraînant une instabilité numérique qui vient perturber la convergence de l'approximation.

<figure><img src="image_18.png" /><figcaption>FIGURE 13 – Comparaison des erreurs $H^1$ finales : `uniformSampling` à gauche, `sobolSampling` au centre, `randomSampling` à droite.</figcaption></figure>

La différence qui apparaît entre les erreurs $H^1$ de Richardson et celles des autres méthodes pour $N$ grand est dûe au fait que Richardson ne converge pas dans ces cas-là (c'est-à-dire atteint le nombre maximal d'itérations autorisées, que nous avons fixé à 10000, sans obtenir le seuil de précision souhaité). On le voit sur la figure 14, qui affiche le nombre d'itérations nécessaire à leur convergence pour chacune des méthodes itératives en fonction de $N$ : celui des deux versions de $Gmres$ ne croît que modérément comparativement à la croissance de la taille du système (asymptotiquement en $O(N)$), tandis que celui de Richardson croît très rapidement, de l'ordre de $N^2$, jusqu'à atteindre sans converger le nombre maximum d'itérations autorisé.
<figure><img src="image_19.png" /><figcaption>FIGURE 14 – Test 4 - Méthodes itératives, comparaison du nombre d'itérations.</figcaption></figure>

Nous mesurons également les temps de résolution pour les trois méthodes, à l'aide du package Julia BenchmarkTools.jl. Nous excluons dans cette comparaison les temps d'assemblage des matrices, vecteurs, et produit matrice-vecteur, afin de comparer uniquement les temps de résolution, et présentons ces mesures en figure 15. La croissance asymptotique de Richardson semble en $O(N^2)$, tandis que celles de la résolution directe et de Gmres semblent être entre $O(N^2)$ et $O(N^3)$. Ces résultats correspondent bien à ce à quoi on s'attend :

- Gmres croît entre $O(N^2)$ (coût du produit matrice-vecteur dominant) et $O(N^3)$, en fonction de la densité de la matrice du système. Pour le maillage à huit cellules de ce test, celle-ci est plutôt pleine, puisque chaque bloc de lignes est rempli entre un huitième et trois huitièmes, en fonction de son nombre de voisins. Il est donc normal que nous soyons entre ces deux ordres.

- Pour une matrice dense, la factorisation LU est de l'ordre de $O(N^3)$, et la résolution est ensuite en $O(N^2)$. Là aussi, le résultat obtenu est donc conforme à nos attentes.

Ces remarques nous conduiront à refaire des tests comparatifs avec des matrices très creuses, c'est-à-dire pour des maillages très raffinés, afin d'observer l'évolution du comportement de ces deux méthodes.

Dans le cas présent, nous observons que LU est un peu plus efficace que gmres pour ces systèmes de petite taille (entre 320 x 320 et 4800 x 4800), tandis que l'algorithme de Richardson est de loin le moins performant. Ces tests nous incitent à l'exclure de nos travaux.
<figure><img src="image_20.png" /><figcaption>FIGURE 15 – Comparaison des temps de résolution des trois méthodes.</figcaption></figure>

### 6.2.2 Gmres sparse versus Gmres mvp.

Pour les mêmes paramètres de test, nous comparons à présent les performances de la version de Gmres avec matrice creuse assemblée (Gmres sparse) et de sa version avec produit matrice-vecteur (Gmres mvp). La mesure des temps d'exécution nous permet de constater que la version produit matrice-vecteur est la plus rapide, quoi que de peu (voir figure 16).

<figure><img src="image_21.png" /><figcaption>FIGURE 16 – Gmres sparse vs Gmres mvp, comparaison des temps d'exécution.</figcaption></figure>
Ces résultats peuvent étonner, au regard de la comparaison du nombre d'itérations nécessaires pour atteindre le seuil de résidu désiré (voir figure 17). On y voit clairement que *Gmres mvp* nécessite plus d'itérations que *Gmres sparse* pour y parvenir, la croissance y étant linéaire avec la taille du système. Cette différence peut s'expliquer par un ordre différent dans les opérations effectuées au sein des itérations, ce qui a une incidence sur le paramètre restart de l'algorithme. Celui-ci entraînant une perte d'informations une fois le nombres de cycles correspondant atteint, cela retarde la convergence de l'algorithme dans la version *Gmres mvp*. On peut observer ce phénomène sur l'évolution des résidus dans les deux versions, figure 18. Le temps global d'exécution de la version produit matrice-vecteur étant pourtant meilleur, on déduit que le coût de chaque produit matrice-vecteur est bien moindre que celui de la version *sparse*.

<figure><img src="image_22.png" /><figcaption>FIGURE 17 – *Gmres sparse* vs *Gmres mvp*, comparaison des nombres d'itérations.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIGURE 18 – *Gmres sparse* vs *Gmres mvp*, comparaison des historiques de convergence du résidu.</figcaption></figure>

Nous poursuivons notre comparaison des performances de ces deux versions en affichant (voir figure 19) la quantité totale de mémoire allouée au cours des simulations (comme pour les autres tests, il s'agit de la médiane prise sur l'ensemble des occurences des tests effectués par *BenchmarkTools*).
<figure><img src="image_24.png" /><figcaption>FIGURE 19 – *Gmres sparse* vs *Gmres mvp*, quantité totale de mémoire allouée durant les simulations, première version.</figcaption></figure>

Cela nous permet de nous rendre compte d'un problème dans la fonction réalisant le produit matrice-vecteur de notre code : celle-ci réalise un très grand nombre d'allocations mémoire, alors que dans l'idéal il ne devrait y en avoir aucune. Nous réimplémentons donc notre fonction, d'une part en préallouant toute les vecteurs nécessaires aux calculs intermédiaires, et d'autre part en faisant directement appel aux routines BLAS2 gemv!. Nous parvenons ainsi à éliminer la quasi-totalité des allocations mémoire, et obtenons des performances similaires à la version sparse (voir figure 20).

<figure><img src="image_25.png" /><figcaption>FIGURE 20 – *Gmres sparse* vs *Gmres mvp*, quantité totale de mémoire allouée durant les simulations, version *Gmres mvp* avec BLAS et préallocations.</figcaption></figure>

Malgré ce gain, nous avons été surpris de ne pas constater d'amélioration des temps d'exécution dans la version produit matrice-vecteur. Ceci dit, puisque les deux versions ont à présent le même coût mémoire et que *Gmres mvp* est plus rapide, c'est celle-ci que nous choisissons de
conserver pour la suite de nos travaux.

Nous comparons ainsi son coût total en allocations mémoire à celui de la méthode LU (voir figure 21). On y voit que le coût de gmres croît en $O(N)$ avec l'augmentation de la taille du système, tandis que celui de LU croît en $O(N^2)$, Gmres semble donc plus adapté à traiter de grands systèmes.

<figure><img src="image_26.png" /><figcaption>FIGURE 21 – LU vs Gmres mvp, quantité totale de mémoire allouée durant les simulations, version Gmres mvp avec BLAS et préallocations.</figcaption></figure>

### 6.2.3 Second test ($\kappa = 100$).

Notre objectif étant d'atteindre le régime haute fréquence, dans lequel la longueur d'onde de l'équation est très petite devant la taille du domaine, nous poursuivons notre étude par un test comparatif entre les deux méthodes conservées précédemment, LU et Gmres mvp, cette fois avec un nombre d'onde plus élevé fixé à $\kappa = 100$. Les autres paramètres de ce test restent inchangés par rapport au test précédent (voir en début de section précédente). Nous en profitons aussi pour comparer un peu plus les différents modes de sampling utilisés pour construire les espaces d'approximation.

En observant l'erreur $H^1$ finale en fonction du nombre d'ondes $N$ par cellule, on remarque qu'il semble nécessaire d'avoir $N \ge 2\kappa$ pour espérer obtenir une précision correcte, quel que soit le type de sampling (voir figure 22). On observe aussi que, pour un nombre moyen d'ondes par cellule, la précision atteinte est meilleure avec randomSampling, mais que les méthodes atteignent le même ordre de précision pour les grandes valeurs de $N$.

Les nombres d'itérations de Gmres croissent moins vite que linéairement avec l'augmentation du nombre d'ondes par élément, entre $O(N^{1/2})$ pour randomSampling et $O(N^{2/3})$ pour uniformSampling (voir figure 23). Le nombre d'itérations nécessaire à la convergence est cependant asymptotiquement le même pour les différents types de sampling.
<figure><img src="image_27.png" /><figcaption>FIGURE 22 – Erreurs $H^1$ finales, uniform sampling (à gauche), random sampling (à droite).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 23 – Nombre d’itérations de Gmres, uniform_sampled (à gauche), random_sampled (à droite).</figcaption></figure>

On présente enfin les performances de LU et Gmres, en temps d'exécution et quantité totale de mémoire allouée, pour la version *uniform_sampled* (voir figure 24). Les performances de la version *random_sampled* sont équivalentes, très légèrement meilleures pour les deux solveurs et les deux mesures. Nous voyons que pour ce nombre d'ondes $\kappa = 100$, Gmres prend le dessus sur LU en temps d'exécution, quoi que les temps des deux méthodes croissent de même entre $O(N^2)$ et $O(N^3)$ comme auparavant, ce qui est normal puisque les tailles des systèmes sont similaires au test précédent. Pour cette même raison, les allocations mémoire sont identiques au test précédent, au net avantage de Gmres, qui reste donc notre meilleur candidat pour nous attaquer à des problèmes plus difficiles.
<figure><img src="image_29.png" /><figcaption>(a) Temps d'exécution
FIGURE 24 – $\kappa = 100$, performance de LU vs Gmres.</figcaption></figure>

## 6.3 Tests de convergence en h.

### 6.3.1 Présentation

On souhaite ici observer le comportement des solveurs lorsque l'on raffine le maillage. Pour cette version séquentielle, on se limite à raffiner trois fois le maillage de référence `ref_msh`, le nombre de cellules étant multiplié par quatre à chaque raffinement, comme montré sur la figure 25.

<figure><img src="image_30.png" /><figcaption>FIGURE 25 – Raffinement du maillage : 8, 32, 128 et 512 cellules.</figcaption></figure>

### 6.3.2 Paramètre restart de Gmres

Une de nos difficultés lors de ces tests a été que, dès le maillage raffiné une première fois, Gmres stagnait parfois. Plusieurs expériences nous ont permis de réaliser que le paramètre `restart` de gmres était en cause dans cette stagnation, et qu'il semblait être nécessaire, pour que l'algorithme converge, que ce paramètre soit croissant proportionnellement au raffinement. Sa valeur par défaut est de 20, et nous avons abouti à la formule heuristique suivante qui fonctionnait dans tous nos tests, où sa croissance est proportionnelle à la racine carrée de $lc$, la longueur caractéristique du maillage : si $\Omega_1$ est notre maillage de référence à huit cellules et $\Omega_k$ un raffinement de $\Omega_1$, avec donc $lc(\Omega_k) < lc(\Omega_1)$, alors on définit `restart` par

$$ restart = \text{Int}\left(20 \times \frac{\sqrt{lc(\Omega_1)}}{\sqrt{lc(\Omega_k)}}\right). $$

Pour nos quatre maillages successivement raffinés, cela donne les valeurs présentées ci-dessous.
<table><thead><tr><td>Nombre de cellules</td><td>lc</td><td>restart</td></tr></thead><tbody><tr><td>8</td><td>0.5</td><td>20</td></tr><tr><td>32</td><td>0.25</td><td>28</td></tr><tr><td>128</td><td>0.125</td><td>40</td></tr><tr><td>512</td><td>0.0625</td><td>57</td></tr></tbody></table>

### 6.3.3 Test.

Les paramètres pour ce test sont les suivants :

- $k = 100$ et une solution de type Fourier-Bessel sur les maillages successivement raffinés décrits ci-dessus,

- Le nombre d'ondes planes par cellule est fixé à $N = 200$,

- `smpl_type` = `uniform_sampling` (sauf pour les erreurs $H^1$),

- $\epsilon_{svd} = 1e-10$ et $\epsilon_{solver} = 1e-10$.

L'affichage de l'erreur $H^1$ (voir figure 26) donne la tendance suivante : à part pour le plus petit nombre de cellules, l'erreur est divisée par 4 à 8 lorsque l'on multiplie le nombre de cellules par 4, donc mieux que linéairement. Et l'on remarque que l'erreur obtenue avec la plus grande résolution est légèrement meilleure que la meilleure de notre test de convergence en $N$. On souhaiterait savoir si le degré de précision pourrait être encore mailleur en raffinant d'avantage, mais multiplier de nouveau par 4 le nombre de cellules conduirait à une taille de système, même très creux, qui devient trop grand pour un traitement séquentiel.

<figure><img src="image_31.png" /><figcaption>FIGURE 26 - Raffinement, LU vs Gmres, erreurs $H^1$.</figcaption></figure>

La comparaison des temps d'exécution ainsi que des quantités de mémoire allouée (figure 27) reste à l'avantage de Gmres. Les temps croissent à peu près à même vitesse, un peu moins que $O(N^{3/2})$, mais LU semble croître plus vite pour les plus grands systèmes. La croissance des quantités de mémoire allouée est légèrement moins importante pour Gmres (complètement linéaire) que pour LU.
<figure><img src="image_32.png" /><figcaption>FIGURE 27 – Test de raffinement, temps de résolution et mémoire totale allouée, LU vs Gmres.</figcaption></figure>

Pour ce qui est du nombre d'itérations de Gmres (figure 28), celui-ci croît bien moins vite que linéairement avec l'accroissement du nombre de cellules, mieux qu'en $O(N^{1/2})$ et un peu moins bien qu'en $O(N^{1/3})$. Lorsque l'on fait un produit matrice-vecteur, l'information ne se propage que par cellules contiguës, on s'attend donc à ce que le nombre d'itérations soit proportionnel au nombre de cellules. Or en divisant $lc$ par 2, on multiplie par 4 le nombre de cellules, on s'attendrait donc à une augmentation de l'ordre de $O(N^{1/2})$, et le résultat est un peu meilleur.

<figure><img src="image_33.png" /><figcaption>FIGURE 28 – Raffinement, Gmres, nombre d'itérations.</figcaption></figure>

Le nombre d'onde $\kappa = 100$ étant le même dans les deux derniers tests, nous comparons les résultats de celui avec convergence en $h$ avec celui du précédent en $N$ : à degré de précision équivalent (i.e. d'erreur $H^1$), comment se comportent les deux ?

On affiche ce comparatif dans le tableau suivant, en se basant sur les erreurs $H^1$ obtenues dans le présent test et en tâchant de trouver les valeurs équivalentes dans le précédent, avec $K$ le nombre de cellules du maillage et $N$ le nombre d'ondes par cellules.
<table><thead><tr><th></th><th>K</th><th>N</th><th>Taille du système (K × N)</th><th>Temps de résolution</th></tr></thead><tbody><tr><td>5.8e-7</td><td>32</td><td>200</td><td>6400</td><td>≈750s</td></tr><tr><td>5.8e-7</td><td>8</td><td>240</td><td>1920</td><td>≈150s</td></tr><tr><td>1.12e-7</td><td>128</td><td>200</td><td>25600</td><td>≈5200s</td></tr><tr><td>1.12e-7</td><td>8</td><td>360</td><td>2880</td><td>≈400s</td></tr><tr><td>1.4e-8</td><td>512</td><td>200</td><td>102400</td><td>≈39000s</td></tr><tr><td>3.3e-8</td><td>8</td><td>600</td><td>4800</td><td>≈1500s</td></tr></tbody></table>

Il semble que, bien que la matrice du système linéaire devienne plus creuse à mesure que le nombre de cellules augmente, le choix d'un domaine plus grossier mais avec plus d'ondes par cellule soit toujours le plus performant... en tout cas tant que l'on ne souhaite pas obtenir une précision trop importante : dans notre test à huit cellules, nous n'avons pas atteint le même degré de précision que celui obtenu avec le maillage le plus raffiné. Une piste de travail serait de chercher les paramètres N et lc permettant d'obtenir le degré de précision optimal pour un κ donné.

<figure><img src="image_34.png" /><figcaption>FIGURE 29 – k = 100. Test sur un raffinement de maillage à 8, 64 et 256 cellules - Gmres, visualisation des erreurs ponctuelles en valeur absolue sur les maillages et de la partie réelle de l'approximation.</figcaption></figure>

## 6.4 Vers le régime haute fréquence : Gmres, nombre d'ondes planes fixé, k croissant

On se concentre ici sur Gmres, toujours en version produit matrice-vecteur. On raffine le maillage de notre carré unité, de sorte qu'il soit composé de 256 cellules, on fixe le nombre d'ondes planes par cellules à Nb = 400, on conserve comme précédemment $\epsilon_{svd}$ et $\epsilon_{solver}$ à 1e-10,
et l'on fait cette fois varier $\kappa \in [200, 400, 600, 800, 1000]$ (et en conséquence la solution test, de type Fourier-Bessel). Ceci nous permettra pour la première fois d'atteindre un régime haute-fréquence, mais aussi d'atteindre les limites de la version séquentielle de notre code en terme de temps de calcul, mettant en évidence la nécessité d'une version parallélisée.

Soulignons que c'est le calcul de l'erreur $\mathbb{H}^1$ qui devient le coût dominant des simulations, puisque, en conservant la formule $lc = \frac{2\pi}{50\kappa}$ afin de conserver une précision satisfaisante d'estimation de l'erreur par rapport à la longueur d'onde du problème, on obtient une résolution de $lc = 1e-4$ pour $\kappa = 1000$, donc un maillage possiblement de l'ordre du million de cellules pour ce calcul d'erreur.

Concernant l'erreur $\mathbb{H}^1$, sans grande surprise nous perdons en précision à mesure que $\kappa$ devient grand devant $N$ (voir figure 30).

<figure><img src="image_35.png" /><figcaption>FIGURE 30 – Test vers le régime haute-fréquence, erreurs $\mathbb{H}^1$.</figcaption></figure>

A seuil de résidu $\epsilon_{solver}$ fixé, le nombre d'itérations reste cependant étonnamment stable en augmentant $\kappa$, de même que les temps de résolution (voir figure 31). Les systèmes linéaires étant de taille exactement identiques, les temps d'assemblages (ainsi que les quantités de mémoire allouée) restent rigoureusement identiques lorsque $\kappa$ change.

<figure><img src="image_36.png" /><figcaption>FIGURE 31 – Test vers le régime haute-fréquence, nombre d'itération et performance de gmres.</figcaption></figure>

Nous finissons cette section par montrer quelques représentations que nous avons pu réaliser dans *Paraview* des approximations réalisées à ces régimes haute fréquence. Remarquons que, du fait de la grande résolution nécessaire à un affichage précis de ces petites oscillations, une carte graphique puissante est requise pour parvenir à ces visualisations.
<figure><img src="image_37.png" /><figcaption>FIGURE 32 – $\kappa$=400, 256 cellules, Gmres, approximation et erreur en valeur absolue.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIGURE 33 – $\kappa = 700$, approximation par gmres.</figcaption></figure>
# 7 Algorithmes parallèles de résolution

## 7.1 Présentation

Dans le but d'améliorer la performance de notre algorithme, on souhaite exploiter les possibilités liées aux architectures multi-cœurs à mémoire partagée, tel que le cluster Leto du Ljll auquel nous avons accès durant notre stage. On implémente donc une version multithread, aussi bien de nos fonctions d'assemblage du système que du produit matrice-vecteur que l'on fournit à gmres.

L'idée générale du multithreading est de répartir les calculs en des tâches qui s'exécutent simultanément. Dans le contexte de notre problème, on effectue d'abord une partition de notre maillage en groupes de cellules, et l'on assigne à chacun de nos threads les calculs liés à ces cellules. Ceci est rendu possible par le caractère grandement local des calculs que nous devons effectuer, aussi bien lors de l'étape de construction du système que lors de celle de sa résolution. En effet, pour ce qui est de la construction, comme nous l'avons déjà vu :

- la matrice *D* est une matrice bloc-diagonale purement locale et peut donc être construite de manière entièrement locale,

- la construction de la matrice *C* ne nécessite que la connaissance des cellules voisines de chaque cellule,

- la construction du membre de droite *b* est également locale à chaque cellule.

De même, pour ce qui est de la résolution du système, le produit matrice-vecteur que l'on fournit à gmres peut calculer en parallèle les différentes parties du vecteur résultant de ce produit dans l'étape d'algorithmme d'Arnoldi dans lequel il intervient : dans cette situation de mémoire partagée on récupère directement le vecteur mis à jour à chaque itération, sans communication nécessaire comme dans un contexte de mémoire distribuée.

### 7.1.1 Data races

La difficulté principale d'implémentation relative au multithreading réside dans la survenue possible de data races : plusieurs threads essaient d'écrire en même temps au même emplacement de la mémoire (typiquement au(x) même(s) indice(s) d'une structure de données), résultant en un undefined behavior qui fait planter le programme. Dans nos algorithmes, outre cette nécessité de répartir les blocs mémoire auxquels les différents threads doivent écrire, les seules instructions qui peuvent conduire à cette situation sont les *push* ! : la tentative simultanée de plusieurs threads d'ajouter un élément en fin du même *Vector* (et donc de potentiellement effectuer une réallocation de mémoire) est un cas typique de data race.

Afin de contourner ce problème, le plus simple est d'allouer la mémoire nécessaire à toutes nos structures de données au moment de leur initialisation, puis, durant les boucles, d'assigner les valeurs au moyen d'indexation dans les emplacement alloués plutôt que d'utiliser *push!* . En effet, nous connaissons exactement toutes les dimensions nécessaires à ces allocations dès la création de l'espace de Trefftz et du maillage sur lequel il se définit, à savoir :

- le nombre total de cellules,

- les cellules voisines de chaque cellule,

- la dimension de l'espace d'approximation local à chaque cellule.

### 7.1.2 Exemple

Nous donnons ici l'exemple du maillage à huit cellules du carré unité 1, de pair avec l'utilisation de quatre threads parallèles. On illustre la répartition du maillage et des éléments du système dans la figure 34.
<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>FIGURE 34 – Exemple de la répartition des éléments du système entre quatre threads, maillage 1 à huit cellules.</figcaption></figure>

## 7.2 Construction des matrices et vecteurs

### 7.2.1 Construction de $D$

L'écriture de la version parallélisée de $D$ ne pose pas de problème, du fait de sa nature totalement locale, donc bloc-diagonale. On commence par initialiser un vecteur de matrices, de dimensions les tailles locales des espaces d'approximation, puis chaque thread s'occupe indépendamment d'assigner les valeurs des blocs auquel il se voit assigné, la répartition entre les threads pouvant se faire simplement lors de l'initialisation de la boucle d'itération sur les cellules du maillage.

### 7.2.2 Construction de $C$

On commence par définir une fonction qui parcourt le squelette d'Ω afin de définir les trois vecteurs d'indices $I_C, J_C, F_C$, de taille identique, tels que :

la $I_C[i]$-ème cellule partage sa $F_C[i]$-ème arête avec la $J_C[i]$-ème cellule.

On prend soin d'ordonner nos vecteurs par indice croissant dans $I_C$, de sorte que les voisines de la $I_C[i]$-ème cellule soient contiguës dans $J_C$ et $F_C$.

Ces vecteurs nous donnent donc les informations nécessaires pour connaître, pour chaque cellule du maillage, le nombre et les indices de ses cellules voisines. Faisant appel à cette fonction, nous pouvons initialiser le vecteur $C_K$ (de taille le nombre d'indices de $I_C$, c'est-à-dire de blocs denses de $C$) de matrices locales, chacune de dimensions locales adéquates. Chaque thread peut ensuite remplir les blocs qu'il s'est vu assigner, sans risque de data race de même que pour $D$, et nous obtenons une structure COO par blocs ($I_C, J_C, C_K$) identique à la version séquentielle.

### 7.2.3 Construction de $b$

Dans le même but que pour $C$, on définit d'abord une fonction qui parcourt le bord d'Ω afin de définir les deux vecteurs d'indices $I_b$ et $F_b$ de taille identique, tels que la $F_b[i]$-ème face
de la I_b[i]-ème cellule est un bord, pour $1 \le i \le N$. On ordonne également ces deux vecteurs tel que I_b soit par ordre croissant.

La situation est toutefois ici un peu différente que pour C : une même cellule peut contenir plusieurs bords (les cellules 4 et 5 du maillage 1 par exemple), et chacun de ceux-ci va contribuer au même bloc dans b. Il n'est pas question de permettre à différents threads de tenter d'écrire simultanément sur ce même bloc, nous devons donc répartir les calculs par cellule, et non par face du bord.

Pour cela, nous générons d'abord un vecteur qui contient, pour chacune des cellules du maillage, le nombre de ses arêtes de bord. Nous parcourons ensuite le vecteur ainsi créé, en répartissant entre les threads les calculs par cellule pour remplir b, de sorte que si une cellule contient plusieurs faces de bord c'est le même thread qui se charge de tous les calculs qui y sont liés.

## 7.3 Assemblage du système linéaire

### 7.3.1 SVD multithread

Une dernière fonction préalable à l'assemblage du système est celle qui réalise les décompositions en valeurs singulières des blocs diagonaux de D et calcule leur pseudo-inverse. De même que pour D, il n'y a ici aucun risque de data race, et il nous suffit d'initialiser un vecteur de la taille de D pour que chaque thread le remplisse ensuite des SVD régularisées relatives aux blocs de D qu'il s'est vu assigner.

### 7.3.2 Construction sparse de $Id - D^{-1}C$

Nous souhaitons implémenter l'assemblage complet du système dans sa version multithread, avec $M = Id - D^{-1}C$ prenant la forme d'une matrice creuse du package SparseArrays.jl. Elle doit être initialisée avec le format COO : nous devons donc créer les vecteurs I_M, J_M, V_M, correspondant respectivement aux indices de lignes, de colonnes, et aux valeurs correspondantes dans M.

Les difficultés résident ici dans le fait d'écrire de façon concurrente dans chacun de ces trois vecteurs depuis nos différents threads, et dans les relations entre les structures par blocs locaux des objets créés précédemment et les vecteurs globaux à créer, notamment concernant la gestion des indices.

Il nous faut d'abord, comme auparavant, nous assurer que l'allocation de mémoire de nos trois vecteurs soit réalisée avant le produit parallélisé $D^{-1}C$, puis que chaque indice de ces vecteurs ne soit accessible qu'à un seul thread.

Pour allouer correctement la mémoire, puisque la matrice $D^{-1}C$ possède exactement la même structure que la matrice C, on s'appuie sur les vecteurs I_C et J_C provenant de la construction de C, ainsi que sur le vecteur des dimensions locales des espaces d'approximation décrit plus haut : on somme les dimensions de tous les blocs denses de C pour définir la taille totale nécessaire, afin d'allouer la mémoire nécessaire à nos trois vecteurs I_M, J_M, V_M et les initialiser à zéro.

Pour s'assurer de la répartition correcte des calculs entre les threads, on commence par effectuer une partition des vecteurs J_C et C_K, afin de regrouper les informations relatives à chaque cellule : les indices des cellules voisines à la i-ème cellule sont regroupées dans le i-ème bloc de J_C, et les blocs de matrices correspondantes dans le i-ème bloc de C_K. On doit aussi gérer l'offset des indices d'écriture dans les vecteurs COO globaux, et pour cela compter le nombre de cellules voisines pour chaque cellule.

On parcourt alors de façon parallèle l'ensemble des cellules du maillage, et chaque thread, se voyant ainsi confier les calculs d'une ou plusieurs cellules, se charge des calculs liés à l'ensemble de ses voisines. Pour chacune d'elles on calcule la matrice dense correspondante, et l'on répartit
les indices et la valeurs dans les trois vecteurs COO aux indices adéquats obtenus après calcul
de l'offset.

Ayant obtenu les vecteurs COO de la matrice $D^{-1}C$, nous créons une matrice identité creuse
à laquelle nous soustrayons la matrice creuse construite par nos vecteurs, pour obtenir finalement
$I - D^{-1}C$.

La version multithread du membre de droite $D^{-1}b$ se construit sans problème, en répartissant
les calculs bloc-diagonaux sur les threads sans risques de *data race*.

7.3.3 Construction du produit matrice-vecteur pour gmres

On décrit ici notre construction de la fonction implémentant le produit matrice-vecteur pour
la version *matrix-free* de gmres, c'est-à-dire renvoyant, à partir de la donnée de vecteurs *X* et *Y*,
le vecteur *Y* = (Id - D<sup>-1</sup>C)*X*.

La problématique est ici encore la nécessité que le même thread soit le seul à écrire sur même
bloc de Y. Pour cela il doit donc effectuer tous les produits matrice-vecteur des blocs de D⁻¹C
correspondant à ce bloc de Y, c'est-à-dire, pour le i-ème bloc de Y, tous les produits relatifs à
la i-ème cellule.

La méthode choisie est la même que dans la version creuse précédente : on commence par
partitionner les informations relatives à chaque cellule dans les vecteurs J_C et C_K. Au sein de
notre fonction réalisant le produit matrice-vecteur, on implémente la partie identité du produit,
puis on répartit les cellules entre les threads, et on effectue une boucle interne sur les n cellules
voisines : pour chaque i-ème cellule on soustrait successivement les 1 ≤ j ≤ n produits Dᵢ⁻¹Cᵢ,jXⱼ
locaux au i-ème bloc de Y.
# 8 Expériences numériques, algorithmes parallèles

## 8.1 Test de validation

On effectue un premier test avec les mêmes paramètres que le second test de convergence en $N$ de la version séquentielle, afin d'en comparer les résultats : on fixe donc $\kappa=100$ et l'on fait varier le nombre d'ondes planes des espaces locaux entre 40 et 600, avec une solution de type Fourier-Bessel sur le maillage 1. On définit ici à 8 le nombre threads disponibles, car le maillage est composé de huit cellules, et c'est la boucle sur les cellules du maillage qui est parallélisée dans la fonction du produit matrice-vecteur.

Les résultats nous montrent que la version multithread atteint exactement les mêmes erreurs $H^1$ en les mêmes nombres d'itérations que la version séquentielle, pour toutes les valeurs de $N$. Cela nous conforte sur la validité de notre code. Les informations intéressantes sont les différences de performance, que l'on affiche en figure 35. La version parallélisée est globalement plus rapide que la version gmres séquentielle. On constate que pour un nombre faible d'ondes par cellule, la différence est faible, le gain du multithread étant probablement contrebalancé par l'overhead lié à la construction et destruction des threads. Il semble y avoir un régime de gain optimal entre les deux versions pour 160 à 360 ondes environ, puis asymptotiquement le gain se fait moindre par rapport à la taille du système. La quantité totale de mémoire allouée est plus importante dans la version multithread, mais croît plus faiblement que la version séquentielle avec $N$.

<figure><img src="image_41.png" /><figcaption>FIGURE 35 – Test de validation de la version multithread, performance par rapport à la version séquentielle.</figcaption></figure>

## 8.2 Weak scaling

Dans ce contexte de multithreading, le principe du weak scaling est d'augmenter à la fois le nombre de threads et la taille du problème, de sorte que la charge de travail par thread reste constante : d'après la page HPC Wiki sur le sujet ([12]), "pour qu'une application passe à l'échelle parfaitement de manière faible, le travail réalisé par chaque cœur reste le même tandis que la taille de la machine croît, ce qui signifie que l'on résoud progressivement de plus gros problèmes dans le même temps nécessaire à résoudre de plus petits problèmes sur une plus petite machine".

Il n'existe pour le moment pas de méthode à deux niveaux (de type genEO, voir [8])) pour notre méthode de Trefftz. Nous savons que le temps pris pour que l'information se propage est proportionnel à la racine carrée du nombre de cellules (en 2D) : en effet, à chaque produit matrice-vector, elle ne se propage que par cellules voisines. Nous nous doutons donc qu'en l'absence de deuxième niveau le passage à l'échelle ne sera pas bon.
On fixe le nombre d'ondes par élément à 400, et l'on fait varier le nombre de threads et avec
lui le nombre de cellules du maillage et la fréquence du problème, de sorte que d'une simulation
à la suivante les deux premiers soient multipliés par 4 et le dernier par 2. Nous présentons les
grandeurs associées à ces simulations dans les deux tableaux suivants.

<table><thead><tr><th>Threads</th><th>κ</th><th>nb cellules</th><th>taille système</th><th>nb blocs</th><th>nb coefs I - D<sup>-1</sup>C</th><th>% dense</th></tr></thead><tbody><tr><td>1</td><td>50</td><td>8</td><td>3200<sup>2</sup></td><td>16</td><td>2.5632+e6</td><td>≈25%</td></tr><tr><td>4</td><td>100</td><td>32</td><td>12800<sup>2</sup></td><td>80</td><td>1.28128+e7</td><td>≈7.81%</td></tr><tr><td>16</td><td>200</td><td>128</td><td>51200<sup>2</sup></td><td>352</td><td>5.63712+e7</td><td>≈2.15%</td></tr><tr><td>64</td><td>400</td><td>512</td><td>204800<sup>2</sup></td><td>1472</td><td>2.357248+e8</td><td>≈0.56%</td></tr></tbody></table>

<table><thead><tr><th>Threads</th><th>nb coefs par thread</th><th>ratio par thread</th></tr></thead><tbody><tr><td>1</td><td>2.5632+e6</td><td>×1</td></tr><tr><td>4</td><td>3.2032+e6</td><td>≈ ×1.25</td></tr><tr><td>16</td><td>3.5232+e6</td><td>≈ ×1.37</td></tr><tr><td>64</td><td>3.6832+e6</td><td>≈ ×1.44</td></tr></tbody></table>

La connaissance des nombres de blocs denses de $D^{-1}C$ en fonction des maillages ainsi que la connaissance de $N^2 = 400^2$ et des tailles totales des systèmes nous a permis de calculer le nombre total de coefficients non nuls de $I - D^{-1}C$ pour chacun des systèmes, ce qui nous semble la mesure la plus juste dans ce contexte de matrice creuse par blocs. Nous avons ainsi pu obtenir, pour chaque simulation, la charge par thread en nombre de coefficients non nuls, puis son rapport à celle de la première simulation avec un seul thread, qui devrait en théorie rester égal à 1 mais dont on voit ici qu'il est supérieur à 1.

Nous avons donc pondéré nos mesures temporelles des durées de construction et de résolution par ces facteurs pour se ramener à ce cas idéal, ce qui avantage un peu les résultats que nous présentons en figure 36.

Sans surprise puisque nous sommes dans une méthode à un niveau, le passage à l'échelle n'est pas bon. Concernant les temps d'assemblage, les croissances sont globalement linéaires avec celle du nombre de coefficients du système, hormis pour la construction du membre de droite où cette croissance est en racine carrée, ce qui nous paraît logique. Le coût dominant est celui de la construction de $D$, $C$, $b$ et $D_\epsilon$. On observe l'avantage indéniable de la version produit matrice-vecteur sur la version matrice creuse, même si nous n'expliquons pas la pente élevée de la construction du produit matrice-vecteur. De même, il nous semble étrange que l'allocation mémoire grandisse plus que linéairement pour 64 threads, et penchons pour un défaut d'implémentation de notre part. La croissance des temps de résolution est également asymptotiquement linéaire par rapport au nombre de coefficients non nuls de la matrice et ne passe pas non plus faiblement à l'échelle, même si la rapidité de résolution par Gmres nous semble dans la pratique très satisfaisante.
<figure><img src="image_42.png" /><figcaption>FIGURE 36 – *Weak scaling*, temps de construction et de résolution du système en fonction du nombre de coefficients du système.</figcaption></figure>

## 8.3 strong scaling

Nous relançons une série de tests en variant le nombre de threads et le raffinement du maillage exatement comme précédemment, mais cette fois en fixant $\kappa = 200$ en plus de $N = 400$, nous plaçant dans les conditions de *strong scaling*. Comme précédemment, les temps de résolution (voir figure 37) nous montrent que notre méthode à un niveau ne passe pas fortement à l'échelle. Il est intéressant de remarquer que les temps de résolution sont très proches de ceux du test précédent, ne dépendant que très peu de $\kappa$. Nous affichons également les erreurs $H^1$ à droite de la même figure. C'est finalement la mesure la plus intéressante de ce test : lors de la simulation à 64 threads et le maillage le plus raffiné à 512 cellules, l'erreur n'est que de 3.92e-9. Pour la première fois de toutes nos expériences numériques relatives à notre méthode de Trefftz, nous obtenons une erreur $H^1$ inférieure à 1e-8. Cela nous incite bien sûr à creuser cette piste.

<figure><img src="image_43.png" /><figcaption>FIGURE 37 – *Strong scaling*, temps de résolution du système et erreurs $H^1$.</figcaption></figure>
## 8.4 Une recherche de précision maximale

Fort du résultat précédent, nous souhaitons reprendre la question, laissée en suspens à la fin de la partie séquentielle, du plus haut degré de précision atteignable avec notre méthode de Trefftz. Que le nombre d'onde soit ici fixé à $\kappa = 200$ seulement nous donne une marge de manœuvre suffisante pour augmenter $N$ un certain nombre de fois, avant qu'il ne vérifie la condition qui exclut, lors de la construction des espaces locaux, les ondes planes qui entraîneraient un overflow du fait d'un exposant trop grand dans leur exponentielle.

Nous conservons donc $\kappa = 200$ et le maillage raffiné à 512 du test précédent avec 64 threads, et procédons à plusieurs tests en augmentant $N$ jusqu'à atteindre cette limite. Les effets de seuil (liés à la distribution utilisée dans la sélection des paramètres d'ondes planes $\varphi$ et $\zeta$) ainsi que la condition que nous venons d'évoquer dans la formule de construction de l'espace de Trefftz conduisent aux dimensions suivantes pour nos tests :

<table><thead><tr><th>$N$ demandé</th><th>$N$ obtenu</th><th>taille système</th></tr></thead><tbody><tr><td>400</td><td>400</td><td>204800 × 204800</td></tr><tr><td>500</td><td>529</td><td>270848 × 270848</td></tr><tr><td>600</td><td>625</td><td>320000 × 320000</td></tr><tr><td>700</td><td>729</td><td>373248 × 373248</td></tr><tr><td>800</td><td>841</td><td>430592 × 430592</td></tr><tr><td>900</td><td>900</td><td>460800 × 460800</td></tr><tr><td>1000</td><td>1024</td><td>524288 × 524288</td></tr><tr><td>1100</td><td>1088</td><td>557056 × 557056</td></tr><tr><td>1200</td><td>1155</td><td>591360 × 591360</td></tr><tr><td>1300</td><td>1184</td><td>606208 × 606208</td></tr><tr><td>1400</td><td>1216</td><td>622592 × 622592</td></tr></tbody></table>

On constate ici que pour les plus grandes valeurs souhaitées, on atteint la limite des tailles de nos espaces de Trefftz (dépendant de ce maillage et de cette valeur de $\kappa$), ce qui est possiblement un frein à l'obtention d'un degré de précision arbitraire.

Nous affichons les erreurs $H^1$ obtenues pour ces approximations (figure 38), et observons avec plaisir qu'en nous plaçant dans ce régime de taille maximale des espaces locaux pour ce $\kappa = 200$ et ce raffinement du maillage, nous atteignons pour la première fois un ordre de précision de 1e-10.

<figure><img src="image_44.png" /><figcaption>FIGURE 38 – Recherche du degré maximal de précision pour $\kappa = 200$, et 512 cellules. Erreurs $H^1$ et nombres d'itérations en fonction de $N$.</figcaption></figure>

On y voit aussi (à droite) que l'augmentation du nombre de degrés de liberté dans les espaces locaux a peu d'influence sur le nombre d'itérations dont gmres a besoin pour converger, à
problème fixe (on multiplie par trois la taille des espaces locaux entre les premier et dernier test, tandis que le nombre d'itérations croît seulement d'un facteur 1.2 environ).

Il est enfin intéressant d'observer les temps de résolution par gmres (en figure 39) comparativement à ceux pris par sa version séquentielle sur le maillage à huit cellules (figures 15 et 16). Du fait du raffinement, la matrice du système linéaire est très creuse (0.56% environ de densité), et nous nous trouvons enfin dans le régime optimal de cet algorithme dans lequel le temps de résolution croît de manière très proche de $O(N^2)$.

<figure><img src="image_45.png" /><figcaption>FIGURE 39 – Recherche du degré maximal de précision pour $\kappa = 200$, et 512 cellules. Temps de résolution de Gmres en fonction de N.</figcaption></figure>

Ayant trouvé notre degré de précision maximal avec ces paramètres du problème, la question s'est posée de savoir si, en essayant d'approcher un problème un peu plus simple, c'est-à-dire en diminuant $\kappa$, on pouvait atteindre un degré encore supérieur de précision.

Pour le même maillage à 512 cellules, on a donc fixé $\kappa = 100$ afin nous donner une marge supplémentaire pour obtenir plus de degrés de liberté dans nos espaces locaux. Nous avons pu construire un espace avec $N=1225$, mais avons malheureusement obtenu une erreur d'overflow dans le calcul de la svd régularisée de $D$, montrant que la borne supérieure que nous avions fixée comme condition pour empêcher que cela n'advienne était trop haute. On a recommencé avec $\kappa = 50$, mais avons obtenu la même erreur...

Nous n'avons pas poursuivi nos recherches dans ce sens. Si la question reste ouverte d'en généraliser les conditions précises en fonction de $\kappa$, nos tâtonnements auront au moins permis de montrer l'existence d'un ensemble de paramètres avec lesquels il est possible d'atteindre une erreur $H^1$ de 1e-10 avec notre méthode de Trefftz.

## 8.5 Somme aléatoire de Fourier-Bessels

A chaque simulation, nous nous sommes jusque là contentés d'essayer d'approcher une seule fonction de Fourier-Bessel. On construit à présent une solution un peu plus riche, en générant une somme de fonctions de Fourier-Bessel, dont les origines sont disposées aléatoirement dans notre domaine. De même, leur ordre (que l'on garde entier ici pour éviter les singularités) et leur poids sont générés aléatoirement (on utilise une distribution uniforme pour tous les paramètres).

Nous gardons $N=1000$ et le maillage à 512 cellules tour les deux simulations, pour rester dans le régime vraisemblablement optimal de notre méthode. Les mesures sont présentées ci-dessous, ainsi que des représentations des approximations obtenues, et l'on constate avec satisfaction que le degré de précision atteint reste le même que dans nos simulations précédentes.
<table><thead><tr><th>κ</th><th>taille du système</th><th>nb de F-B</th><th>erreur H<sup>1</sup></th><th>nb d'itérations</th><th>temps de résolution</th></tr></thead><tbody><tr><td>200</td><td>524288×524288</td><td>100</td><td>7.38e-10</td><td>436</td><td>107s</td></tr><tr><td>400</td><td>491520×491520</td><td>200</td><td>2.23e-9</td><td>368</td><td>76s</td></tr></tbody></table>

<figure><img src="image_46.png" /><figcaption>FIGURE 40 – κ = 200 et 512 cellules, approximation d'une somme de 100 fonctions de Fourier-Bessel de coefficients aléatoires. Partie réelle à gauche, erreur ponctuelle à droite (comprise entre 3.3e-12 et 8.6e-9).</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>FIGURE 41 – κ = 400 et 512 cellules, approximation d'une somme de 200 fonctions de Fourier-Bessel de coefficients aléatoires. Partie réelle à gauche, erreur ponctuelle à droite (comprise entre 9.2e-12 et 1.9e-7).</figcaption></figure>
# 9 Problème de diffraction sur un écran

## 9.1 Présentation - problème aux limites et formulation variationnelle

On réfère ici aux présentations faites dans [7, chapitre 2] ou [11, section 2]. On place un écran plat $\Gamma_s$ (c'est-à-dire un segment) dans notre domaine $\Omega \subset \mathbb{R}^2$ de référence. Soient une onde incidente $u_i \notin L^2(\mathbb{R}^2)$, et $u_s$ l'onde inconnue réfléchie sur $\Gamma_s$ (s pour scattered). La solution totale $u_t = u_i + u_s$, satisfait le problème de Dirichlet suivant :

$$ egin{cases} -\Delta u_t - \kappa^2 u_t = 0 & \text{dans } \mathbb{R}^2 \setminus \Gamma_s, \\ u_t = 0 & \text{sur } \Gamma_s, \\ \lim_{|\mathbf{x}| \to \infty} |\mathbf{x}| \left( \frac{\partial u_s(\mathbf{x})}{\partial |\mathbf{x}|} + i\kappa u_s(\mathbf{x}) \right) = 0, \end{cases} $$

où la dernière condition sur $u_s$, appelée condition de Sommerfeld, exprime l'idée que "$u_s$ sort du domaine" et assure l'unicité de la solution, et peut être approchée par une condition de Robin homogène $(\partial_n + i\kappa)u_s = 0$. L'onde incidente $u_i$ ne satisfait pas cette condition mais est solution de l'équation de Helmholtz. La condition (de type sound-soft) qui impose que la pression de l'onde totale est nulle sur l'écran implique que $u_s = -u_i$ sur $\Gamma_s$. Nous cherchons donc à présent la solution $u_s$ au problème

$$ egin{cases} -\Delta u_s - \kappa^2 u_s = 0 & \text{dans } \Omega, \\ u_s = -u_i & \text{sur } \Gamma_s, \\ (\partial_n + i\kappa)u_s = 0 & \text{sur } \partial\Omega. \end{cases} $$

Nous formalisons ce problème à l'aide du problème équivalent d'origine (0.1) de [5] :

$$ egin{cases} -\Delta u - \kappa^2 u = 0 & \text{dans } \Omega, \\ (\partial_\nu + i\kappa)u = t(-\partial_\nu + i\kappa)u + g & \text{sur } \Gamma, \end{cases} $$

avec des choix particuliers de $t$ et $g$ qui nous permettent de satisfaire le problème de diffraction que l'on cherche à modéliser, tout en rendant son implémentation aisée. En effet,

$$ \text{— sur } \Gamma = \partial\Omega, \begin{cases} t = 0 \\ g = 0 \end{cases} \Rightarrow (\partial_n + i\kappa)u_s = 0, $$

$$ \text{— sur } \Gamma = \Gamma_s, \begin{cases} t = -1 \\ g = -2i\kappa u_i \end{cases} \Rightarrow (\partial_n + i\kappa)u_s = (\partial_n - i\kappa)u_s - 2i\kappa u_i \Rightarrow u_s = -u_i, $$

et les deux conditions, sur l'écran et le bord extérieur, sont alors vérifiées.

Nous devons simplement modifier en conséquence les intégrales définissant les coefficients de la matrice C et du vecteur second membre b :

- Le coefficient $C_{k,j}^{l,m}$ situé à la l-ème ligne et m-ème colonne du bloc situé à la k-ème ligne de bloc et j-ème colonne de bloc de C s'écrit à présent

$$ C_{k,j}^{l,m} = \int_{\Sigma_{kj}} (\partial_{\nu_k} + i\kappa) e_{jm} \overline{(\partial_{\nu_k} + i\kappa) e_{kl}} + \int_{\Gamma_k} t(-\partial_{\nu_k} + i\kappa) e_{km} \overline{(\partial_{\nu_k} + i\kappa) e_{kl}}, $$
avec $$ \begin{cases} t = 0 & \text{si } \Gamma_k \subset \partial\Omega, \\ t = -1 & \text{si } \Gamma_k \subset \Gamma_s, \end{cases} $$
— la contribution $b_{k,l}$ au coefficient $b_k$ situé en $k$-ème indice du vecteur $b$, s'écrit

$$b_{k,l} = \int_{\Gamma_k} g(\overline{\partial_{\nu_k} + i\kappa})e_{kl}, \text{ avec } \begin{cases} g = 0 & \text{si } \Gamma_k \subset \partial\Omega, \\ g = -2i\kappa u_i & \text{si } \Gamma_k \subset \Gamma_s. \end{cases}$$

## 9.2 Maillages et espaces locaux

### 9.2.1 Génération et raffinement local des maillages

Pour générer les maillages liés à ce problème, nous avons utilisé l'API Julia de GMSH. Outre le fait de générer un maillage à partir d'une géométrie, GMSH permet d'effectuer des raffinements locaux en des points, arêtes ou surfaces choisis dans les maillages, avec des modes de raffinement pleinement définissables par l'utilisateur.

Nous nous concentrons pour nos tests sur un algorithme de raffinement local de GMSH appelé Distance, qui interpole linéairement la taille d'une cellule entre une taille minimale ($lc_min$) et une taille maximale ($lc_max$) en fonction de sa distance à un point : si celle-ci est inférieure à $dist_min$ la taille sera de $lc_min$ et si elle est supérieure à ($dist_max$), elle sera de $lc_max$.

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

## 9.3 Implémentation et premiers tests de validation

Nous avons réalisé une nouvelle version de nos codes d'assemblage et de résolution du système linéaire, de sorte qu'il prenne en compte les modifications de la formulation variationnelle, et qu'il fonctionne avec des tailles variables d'espaces locaux. Pour prendre en compte la formulation variationnelle, nous avons ajouté à nos fonctions générant les vecteurs d'indices du squelette et du bord (respectivement I, J, F et I, F, décrites dans la partie d'implémentation multithread), la possibilité de générer également un vecteur T de flag correspondant au t de la formulation variationnelle, indiquant ainsi si la $F[i]$-ème face de la $I[i]$-ème cellule fait partie de l'écran (si $T[i] = -1$), ou pas (si $T[i] = 0$). Nous nous servons ensuite de cette information pour construire $C$ et $b$ de façon adéquate.

Nous effectuons des premiers essais de validation de nos codes de construction des maillages et du système linéaire avec écran. Nous choisissons pour ces tests une onde incidente de type onde plane propagative (le vecteur de direction est à coefficients réels). La direction choisie ici est de $\frac{5}{8}\pi$.

Le premier essai est effectué sur un simple écran horizontal, situé au centre de notre carré entre les points d'abscisse 0.25 et 0.75, qui nous servira plus loin dans nos tests de convergence. Nous fixons $\kappa = 400$ et 400 ondes par cellules, et validons nos fonctions d'assemblage par résolution directe (LU), comme nous l'avons fait précédemment.

A la différence des tests des sections précédentes, nous ne connaissons pas ici la solution exacte, et ne pouvons pour le moment qu'essayer de déduire de l'observation de l'approximation obtenue si notre code est correct. Le résultat obtenu est présenté en figure 43. On y voit à
gauche l'onde incidente, c'est-à-dire l'onde plane servant de terme source dans la formulation variationnelle, au centre l'approximation obtenue de l'onde réfléchie, et à droite la somme des deux (partie réelle). Il nous semble que cette onde totale reproduit correctement le comportement attendu : l'onde incidente, venant de la droite, se reflète sur la partie supérieure de l'écran et rebondit vers la partie supérieure gauche du domaine, on y voit bien la superposition des ondes incidente et réfléchie. Sous l'écran au contraire on observe bien "l'ombre" attendue, c'est-à-dire la zone que l'onde incidente n'atteint pas puisqu'elle ne traverse pas l'écran.

<figure><img src="image_49.png" /><figcaption>FIGURE 43 – Ecran simple, premier test (résolution directe), $\kappa = 400$, 400 ondes par élément. Onde incidente, approximation de l'onde réfléchie, et somme des deux.</figcaption></figure>

Pour notre second essai, nous créons un écran en forme de spirale, afin d'essayer de modéliser le phénomène de *trapping*. Le principe en est que des ondes sont prises entre deux éléments réfléchissants parallèles, et rebondissent indéfiniment. Le maillage créé est montré en figure 45. Nous utilisons cette géométrie pour tester cette fois une résolution par Gmres. Le nombre d'onde du problème est toujours de $\kappa=400$, avec 400 ondes par élément. Il est intéressant de noter que Gmres converge en tout de même 918 itérations, soit bien plus (environ trois fois) que pour tous les tests des sections précédentes. Ceci confirme la difficulté de résolution dans ce contexte. Nous montrons en figure 46 des visualisations de l'approximation obtenue, en valeur absolue, et de la partie réelle de la somme totale. Comme dans le premier test on y voit bien l'ombre pauvre en ondes à gauche de l'écran, et les multiples réflexions dans la spirale, dont l'addition provoque une augmentation de l'amplitude totale (dans les zones aux teintes plus sombres).

<figure><img src="image_50.png" /><figcaption>FIGURE 44 – Ecran simple, approximation obtenue pour $\kappa = 16$. Vue sur la singularité du gradient.</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>FIGURE 45 – Ecran spirale, maillage.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>FIGURE 46 – Ecran spirale, $κ = 400, 400$ ondes par élément. Onde réfléchie en valeur absolue (à gauche), partie réelle de la somme de l'onde incidente et de l'onde réfléchie (à droite).</figcaption></figure>

## 9.4 Tests de convergence

Le but des tests suivants est d'estimer s'il est possible d'obtenir avec notre méthode de Trefftz un degré de précision arbitraire à un problème de diffusion avec écran, dans la configuration suivante : nous prenons le premier écran horizontal simple, et définissons l'angle de l'onde incidente égal à $π$. L'onde arrive donc de la droite du domaine, exactement dans l'axe de l'écran, et ne se réfléchit que sur son point d'extrémité droite. Ce problème est réputé difficile du fait cette extrême concentration en ce point, et de la singularité du gradient de l'approximation qui en résulte au contact de l'onde avec l'extrémité de l'écran. En conséquence, ce sont les zones autour
des points correspondant aux extrémités de l'écran que nous souhaiterons raffiner.

<figure><img src="image_53.png" /><figcaption>FIGURE 47 – Test 1, visualisation du résultat, approximation de référence. A gauche l'onde incidente, au centre l'approximation par gmres, à droite la somme des deux (parties réelles).</figcaption></figure>

### 9.4.1 Méthode d'évaluation de la convergence

Comme nous l'avons déjà évoqué, nous ne connaissons pas la solution exacte à ce problème, il nous faut donc construire une approximation de référence, pour laquelle nous choisirons un nombre élevé d'ondes par cellule et un grand raffinement du maillage. L'estimation de convergence, c'est-à-dire l'erreur, sera ensuite réalisée pour chaque approximation par rapport à cette référence.

### 9.4.2 Maillages différents et plus proches voisins

Une difficulté est qu'ici un maillage différent est généré pour chaque approximation. Ainsi, pour pouvoir évaluer l'erreur, nous devons, pour chaque point d'évaluation, trouver à quelle cellule il appartient dans l'approximation à évaluer d'une part, et dans l'approximation de référence d'autre part, afin d'utiliser les espaces locaux d'approximation adéquats aux deux cellules correspondantes, et donc d'évaluer l'erreur entre les bonnes solutions locales correspondant à ce point.

Du point de vue de l'implémentation pratique de ce problème, nous utilisons un KD-Tree, une structure de graphe qui permet de trouver efficacement les plus proches voisins d'un point dans un nuage de points : on commence par construire une structure dictionnaire, dont les clés sont les nœuds du maillage et les valeurs sont les indices des cellules du maillage dont le nœud correspondant est un sommet. On passe l'objet construit à une fonction qui, à partir d'un nuage de points et de ce dictionnaire, retourne les indices des cellules du maillage dans lesquelles les points de ce nuage se situe. Pour chacun de ces points, on peut ensuite évaluer l'approximation et la référence, en utilisant les espaces locaux adéquats correspondant aux cellules respectives des maillages, puis calculer enfin les erreurs.

### 9.4.3 Erreurs

Nous disons les erreurs car, au lieu de calculer l'erreur $\mathbb{H}^1$ dont nous avons vu qu'elle était très coûteuse à calculer, nous nous contenterons de calculer les erreurs $||u_{ref}(x)| - |u_{approx}(x)||$ ainsi que $||\nabla u_{ref}(x)| - |\nabla u_{approx}(x)||$.

Ce calcul d'erreur sera de plus effectué en nous concentrant sur un segment dans l'axe de l'écran, en partant de l'extrémité droite du domaine et en nous approchant du bout de l'écran : nous prenons pour cela 1000 points, dont les abscisses sont représentées sur la figure 48. On y
montre aussi un zoom progressif sur les cellules du maillage de référence sur lesquels se trouvent ces points d'évaluation.

<figure><img src="image_54.png" /><figcaption>FIGURE 48 – Abscisses des points d'évaluation des erreurs (à gauche), Zooms sur les cellules dans lesquels se trouvent les points d'évaluation des erreurs (à droite).</figcaption></figure>

### 9.4.4 Tailles des espaces locaux d'approximation

Comme nous l'avons déjà évoqué, le nombre d'ondes par cellule de nos bases de Trefftz est limité de manière intrinsèque par le rayon du cercle circonscrit à la cellule (tout en dépendant également du nombre d'onde de l'équation $\kappa$), du fait même la recette de construction des espaces d'approximation.

Afin de pouvoir raffiner autant que souhaité tout en maximisant le nombre d'ondes disponible pour chaque cellule, nous avons donc adapté le code de construction de l'espace de Trefftz : au lieu de définir un nombre $N$ d'ondes par cellule fixe, on décide d'adapter la formule existante, et de définir $N_{\Omega_k}$ pour chaque cellule $\Omega_k$ en fonction du rayon $r_{\Omega_k}$ du cercle qui lui est circonscrit. Nous verrons que des choix différents pour cette dépendance amèneront des résultats très différents à nos tests.

## 9.5 Test 1

### 9.5.1 Paramètres

Pour cette première expérience, nous raffinons, en utilisant la méthode Distance de GMSH décrite plus haut, en fixant $lc = 0.3$ la taille maximale des cellules, $dist\_max = 0.25$, en définissant $dist\_min = lc\_min$, avec $lc\_min$ que nous réduisons progressivement de presque deux ordres de grandeur, avec pour $0 \le n \le 4$, $lc\_min = \frac{0.3}{(2\sqrt{2})^n}$, ce qui donne les tailles de maillages et de systèmes suivants.

<table><thead><tr><th>simulation</th><th>lc<sub>min</sub></th><th>nb cellules</th><th>taille système</th></tr></thead><tbody><tr><td>1</td><td>0.3</td><td>36</td><td>7056×7056</td></tr><tr><td>2</td><td>0.1061</td><td>70</td><td>10096×10096</td></tr><tr><td>3</td><td>0.0375</td><td>114</td><td>11956×11956</td></tr><tr><td>4</td><td>0.0133</td><td>216</td><td>15288×15288</td></tr><tr><td>réf</td><td>0.0047</td><td>272</td><td>19482×19482</td></tr></tbody></table>

Concernant la dimension des espaces locaux d'approximation, notre premier choix est que, pour chaque cellule $\Omega_k$, la dimension de l'espace local $N_{\Omega_k}$ soit linéairement proportionnelle au
rayon $r_{\Omega_k}$ du cercle qui lui est circonscrit, ainsi qu'au nombre d'onde de l'équation $\kappa$ et à un paramètre utilisateur $\alpha$. Nous définissons donc dans la fonction de construction de l'espace de Trefftz la formule suivante : pour chaque cellule $\Omega_k$,

$$N_{\Omega_k} = \alpha \kappa r_{\Omega_k}.$$

### 9.5.2 Premières observations

Nous affichons en figures 50 et 49 les maillages raffinés successivement (l'écran est en rouge), ainsi que les formes des matrices C correspondantes (ou $D^{-1}C$, qui ont la même structure). On y observe nettement la différenciation progressive de la taille des espaces locaux : sur le maillage le plus grossier, les blocs de la matrice correspondante sont tous de taille identique, puisque, les cellules étant toutes de même taille, les espaces locaux ont tous le même nombre d'ondes planes dans leur base. A mesure que l'on raffine, les espaces locaux correspondant aux plus petites cellules voient la taille de leur base diminuer, et ainsi la hauteur (i.e. le nombre lignes) du bloc correspondant. Pour une cellule dont une de ses voisines est plus petite qu'elle (resp. plus grande), le bloc correspondant à leur interface est ainsi un rectangle allongé verticalement (resp. horizontalement). Nous avons également affiché en rouge sur ces représentations des matrices les blocs diagonaux correspondant aux cellules pour lesquelles la condition de bord liée à l'écran de la formulation variationnelle entre en jeu.

<figure><img src="image_55.png" /><figcaption>FIGURE 49 – Test 1, approximation de référence, à gauche, son maillage, à droite, la matrice correspondante.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>FIGURE 50 – Test 1, maillages raffinés et matrices correspondantes des simulations 1 à 4.</figcaption></figure>
9.5.3 Résultats

On observe en figure 51 les erreurs des approximations et de leurs gradients par rapport à notre référence. Pour ce qui est des erreurs ponctuelles des approximations (en haut sur la figure), si l'approximation la moins raffinée semble bien avoir globalement la plus importante et si la plus raffinée semble bien avoir la plus faible (le gain étant faible), la seconde semble plus proche de la référence que la troisième, ce devrait être l'inverse. De plus, si l'on prend le maximum de ces erreurs sur l'ensemble des points d'évaluation, qui se situe à proximité du bord de l'écran, c'est la seconde qui est la plus proche de la référence, ce que l'on peut observer notamment en haut à droite. Pour ce qui est des erreurs des gradients (en bas de la figure), le comportement semble également erratique : si l'on voit une amélioration globale entre la première et la dernière approximation (à gauche), le maximum des erreurs des gradients sur les points d'évaluation nous montre que les deux intermédiaires sont plus précises, ce que l'on voit notamment à approche de l'extrémité de l'écran (à droite). En résumé, les erreurs des approximations comme de leurs gradients restent globalement du même ordre malgré le raffinement.

<figure><img src="image_57.png" /><figcaption>FIGURE 51 – Test 1, Erreurs par rapport à l'approximation de référence des approximations (en haut) et de leurs gradients (en bas). Zoom progressif sur la singularité.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>FIGURE 52 – Test 1, tailles des espaces locaux des différentes approximations.</figcaption></figure>

La raison se trouve en figure 52, présentant la taille des espaces locaux en fonction de $r_{\Omega_k}$. Celle-ci nous montre que la dépendance est bien globalement linéaire comme espéré (les paliers sont le fait de la recette de construction des espaces de Trefftz). Mais elle nous renseigne aussi sur la taille minimale de ces espaces locaux, pour les plus petites cellules : du fait de cette dépendance, nous avons perdu un ordre et demi environ pour cette grandeur, passant de 196 pour le maillage grossier à 4 seulement pour les plus petites cellules du maillage de référence. Or ces cellules se situent précisément aux points de singularité du bord de l'écran. Ce peu de fonctions de bases à ces endroits nous paraît expliquer que l'erreur y soit particulièrement élevée, et ne diminue pas avec les raffinements. Il nous semble nécessaire de ne pas trop perdre de degrés de liberté dans nos approximations à ces endroits-là, nous décidons donc de réduire la dépendance de $N_{\Omega_k}$ envers la taille de cellules dans le test suivant. Ceci dit, quoi qu'il ne permette pas la convergence espérée, le fait que la précision de nos approximations reste à peu près stable avec cette diminution drastique de la taille des espaces locaux sur les plus petites cellules nous semble un résultat satisfaisant.

## 9.6 Test 2

### 9.6.1 Paramètres

Nous modifions la dépendance du nombre d'ondes $N_{\Omega_k}$ par cellule par rapport au rayon $r_{\Omega_k}$ du cercle circonscrit à celle-ci, de sorte qu'elle soit en racine carrée, et par homogénéité nous y incluons $\kappa$, pour obtenir la formule suivante :

$$N_{\Omega_k} = \alpha \sqrt{\kappa r_{\Omega_k}}.$$

Pour nous placer en contexte plus favorable, nous choisissons également une pente plus douce pour l'évolution de la longueur caractéristique minimale dans le raffinement des maillages, en définissant, pour $0 \le n \le 5$, $lc_{min} = \frac{0.2}{2^{n/2}}$. Ceci réduira encore la diminution de la taille des espaces locaux des plus petites cellules. On fixe $lc_{max} = 0.2$ et $dist_{max} = 0.2$. Nous présentons les grandeurs associées aux simulations dans le tableau ci-dessous. L'écart de taille du système de
référence par rapport aux autres s'explique par le choix de $\alpha = 10$ pour plus de précision, contre $\alpha = 5$ pour les approximations.

<table><thead><tr><th>simulation</th><th>lc<sub>min</sub></th><th>nb cellules</th><th>taille système</th></tr></thead><tbody><tr><td>1</td><td>0.2</td><td>70</td><td>4820 x 4820</td></tr><tr><td>2</td><td>0.1414</td><td>80</td><td>5460 x 5460</td></tr><tr><td>3</td><td>0.1</td><td>98</td><td>6564 x 6564</td></tr><tr><td>4</td><td>0.0707</td><td>116</td><td>7347 x 7347</td></tr><tr><td>5</td><td>0.05</td><td>132</td><td>7878 x 7878</td></tr><tr><td>réf</td><td>0.0354</td><td>162</td><td>15876 x 15876</td></tr></tbody></table>

<figure><img src="image_59.png" /><figcaption>FIGURE 53 – Test 2, maillages.</figcaption></figure>

Les maillages et matrices associés sont représentés en figures 53. On y voit que la différenciation des espaces locaux est moindre que dans le premier test, ce que la figure 54 présentant la taille des espaces locaux nous confirme : entre les premier et dernier maillages, la longueur caractéristique minimale des cellules est divisée par 5 à 6, tandis que le nombre de fonctions de base des cellules correspondantes ne l'est que par deux. Nous nous trouvons donc dans la situation favorable où, relativement à $r_{\Omega_k}$, le nombre de degrés de liberté est plus important pour les plus petites cellules qui sont celles proches de la singularité, nous pouvons donc espérer un gain de précision.

<figure><img src="image_60.png" /><figcaption>FIGURE 54 – Test 2, taille des espaces locaux d'approximation en fonction de $r_{\Omega_k}$.</figcaption></figure>
## 9.6.2 Résultats

Nous obtenons ce qui nous semble aller dans le sens d'une convergence : la figure 55 montre (colonne de gauche, zooms sur la singularité du bord de l'écran) qu'avec les raffinements successifs l'approximation tend à s'annuler en se rapprochant de cette extrémité, de même que son gradient tend à croître en cette zone. La même tendance se dégage des évolutions des erreurs ponctuelles (colonne du milieu), quoi que l'on n'y observe toujours pas de tendance nette entre les trois premières simulations. C'est surtout la visualisation de l'évolution des maximums de ces erreurs ponctuelles en fonction de $lc\_min$ (colonne de droite) qui nous permet de conclure à cette amélioration.

L'insatisfaction domine à l'issue de ces deux tests. Outre qu'il n'est pas aisé de conserver la maîtrise du raffinement de maillage comme désiré, étant tributaire des algorithmes de GMSH (ce que la pratique de cet outil nous permettra d'améliorer), les disparités sont trop grandes dans les tailles d'espaces locaux entre ces tests et ceux de fin de la section précédente où nous avions pu atteindre un haut degré de précision. De plus, le fait de ne pas disposer de solution analytique de référence pour mesurer l'erreur commise est une source supplémentaire d'incertitude dans la méthode comme dans les résultats, en plus de compliquer le calcul des erreurs (dans la gestion des deux maillages, de l'approximation et de la référence).

<figure><img src="image_61.png" /><figcaption>FIGURE 55 – Test 2, approximations (en haut), et leurs gradients (en bas). Zooms sur la singularité (à gauche), erreurs sur le segment [0.75, 1.] (au milieu), et évolutions des maximums des erreurs ponctuelles (à droite).</figcaption></figure>

## 9.7 Approximation d'une singularité de coin

Les considérations précédentes nous conduisent à deux choses :

- Premièrement, nous modifions la formule définissant $P$, le paramètre de troncation en Fou-rier utilisé pour la construction des espaces d'approximation (évoqué en 2.4.4), le reliant à $N_{\Omega_k}$ et au rayon de la cellule $\Omega_k$ considérée. Cette modification va nous permettre d'obtenir
bien plus de degrés de liberté dans les petites cellules que dans nos deux tests précédents. Nos essais nous conduisent alors à finalement écarter la dépendance de $N_{\Omega_k}$ en $\sqrt{r_{\Omega_k}}$, pour ne garder que la règle de dépendance linéaire d'origine, dont nous présentons les résultats ici.

- Deuxièmement, afin de pouvoir évaluer plus sereinement nos approximations, nous modifions quelque peu le protocole de notre expérience numérique. Nous abandonnons le problème de diffraction avec écran et nous concentrons sur l'évaluation de la capacité de notre méthode de Trefftz à approcher correctement une singularité. Pour cela nous utilisons à nouveau une fonction de Fourier-Bessel comme solution de référence, mais cette fois avec un ordre non entier de 0.5, créant ainsi une singularité sur son gradient, et en plaçant son origine en (0., 0.). Cela a pour effet de reproduire l'idée d'un point singulier du même type que sur l'extrémité d'un écran plat, mais sur le bord du maillage. Nous conservons par contre le principe d'effectuer nos mesures le long d'un segment qui approche ce point.

Nous effectuons des tests comparatifs sans et avec raffinement, ce qui nous permettra peut-être de mieux cerner la pertinence du raffinement. Les deux maillages considérés sont composés respectivement de 69 et 123 cellules, et nous effectuons pour chacun d'eux 6 simulations en augmentant $\alpha$, le paramètre utilisateur multiplicatif qui intervient dans la taille des espaces locaux : pour $1 \le n \le 6$, on choisit $\alpha_n = 2^{n-1}$ pour la n-ème approximation.

<figure><img src="image_62.png" /><figcaption>FIGURE 56 – Maillages avec points d'évaluation (à gauche) et tailles des espaces locaux en fonction de $r_{\Omega_k}$ (droite). Tests sans raffinement (en haut) et avec raffinement (en bas).</figcaption></figure>
<table><caption>TABLEAU 1 – Simulations en fonction de $\alpha$, versions sans et avec raffinements. Tailles des systèmes, nombre d'itérations de Gmres, maximums des erreurs ponctuelles sur le segment [0., 1.], et de leurs gradients sur le segment [$\epsilon$, 1.].</caption><thead><tr><th rowspan="2">α</th><th colspan="4">sans raffinement</th><th colspan="4">avec raffinement</th></tr><tr><th>taille</th><th>nb iters</th><th>||err(u)||<sub>∞</sub></th><th>||err(∇u)||<sub>∞</sub></th><th>taille</th><th>nb iters</th><th>||err(u)||<sub>∞</sub></th><th>||err(∇u)||<sub>∞</sub></th></tr></thead><tbody><tr><td>1</td><td>2733</td><td>173</td><td>0.189</td><td>511</td><td>3356</td><td>232</td><td>0.19</td><td>550</td></tr><tr><td>2</td><td>5200</td><td>173</td><td>0.183</td><td>497</td><td>6215</td><td>298</td><td>0.135</td><td>477</td></tr><tr><td>4</td><td>10208</td><td>193</td><td>0.154</td><td>488</td><td>11989</td><td>356</td><td>0.096</td><td>422</td></tr><tr><td>8</td><td>19669</td><td>211</td><td>0.1453</td><td>476</td><td>23440</td><td>351</td><td>0.098</td><td>423</td></tr><tr><td>16</td><td>38469</td><td>274</td><td>0.1</td><td>439</td><td>46003</td><td>340</td><td>0.099</td><td>409</td></tr><tr><td>32</td><td>76244</td><td>351</td><td>0.0915</td><td>371</td><td>90467</td><td>448</td><td>0.084</td><td>377</td></tr></tbody></table>

Nous voyons déjà en figure 56 que nous arrivons cette fois à obtenir des tailles d'espaces locaux globalement comparables à celles nous ayant permis d'atteindre un niveau élevé de précision auparavant. Celles des plus petites cellules de la version non raffinée sont presque un ordre de grandeur supérieures à celles des plus petites cellules de la version raffinée, ce qui suggère possiblement de meilleures approximations. Le tableau ci-dessous montre cependant que les tailles des systèmes sont supérieures dans les versions avec raffinement.

Ce tableau montre que les maximums des erreurs ponctuelles des approximations et de leurs gradients évoluent de façon sensiblement identique entre la première et la dernière simulation pour les deux maillages, avec un nombre d'itérations globalement à l'avantage du maillage non raffiné.

Mais la figure 57 présente un tout autre aspect des choses. On y voit les mesures des erreurs d'approximations et gradients le long du segment [0, 1] de la droite réelle, la singularité de la solution étant à l'origine. Et la différence entre les deux maillages est très nette : sur la version raffinée (à droite), la singularité est globalement contenue dans le premier quart du segment, et l'on voit une claire amélioration de l'erreur sur tout le reste de l'intervalle d'une simulation à la suivante, quoi que la précision reste du même ordre dès la troisième simulation ($\alpha = 8$), d'environ $10^{-7}$ pour les approximations et $10^{-6}$ pour leurs gradients. Au contraire, sur la version avec maillage plus grossier (à gauche), on voit l'impact bien plus important de la singularité, qui se propage sur plus de la moitié gauche de l'intervalle de mesure à travers les interfaces entre cellules voisines. Ces erreurs aux interfaces sont très visibles sur les courbes aussi bien des approximations que des gradients par les pics périodiques, et l'on a le sentiment que l'erreur décroît en nombre d'interfaces. Ceci entraîne un degré de précision globalement moins élevé, de même qu'une tendance bien moins nette à la convergence d'une simulation à l'autre.

Grâce à la modification de la formule de dépendance du paramètre $P$ envers $N$, nous avons clairement pu montrer dans ce test l'intérêt de procéder à un raffinement, avec la propriété d'un maillage raffiné de capturer une singularité et la contenir de façon bien plus locale qu'un maillage grossier.
<figure><img src="image_63.png" /><figcaption>FIGURE 57 – Comparaisons des évolutions des erreurs des approximations (en haut) et de leurs gradients (en bas) entre maillage non raffiné (colonne de gauche) et raffiné (colonne de droite).</figcaption></figure>
# 10 Déflation

## 10.1 Présentation

On considère un système linéaire $Ax = b$, avec $A \in \mathbb{C}^{n \times n}$, $x, b \in \mathbb{C}^n$, avec $n \in \mathbb{N}$.

Citons [10, p. 1] : "Dans leur formulation originale, les méthodes de déflation éliminent des valeurs propres du spectre de $A$. Si ce spectre est très éclaté ou si des valeurs propres sont proches de l'origine, cette technique conduit souvent à une convergence plus rapide de la méthode de Krylov [nda : de gmres dans notre cas].[...] Dans les méthodes de déflation, des systèmes singuliers doivent être résolus. En effet, les opérateurs de projection sont les ingrédients majeurs de cette méthode. Le système original est multiplié par une projection qui conduit à un système linéaire singulier, dont le spectre peut être plus favorable à la méthode de Krylov. Les projections sont souvent construites en connaissance du spectre de $A$ et de certains vecteurs propres pertinents. Les espaces de déflation sont constitués de ces vecteurs propres, ou d'approximations de ceux-ci ".

On se donne deux matrices $Y, Z \in \mathbb{C}^{n \times r}$, avec $r \in \mathbb{N}$ tel que $r << n$. Un choix classique, que nous ferons, est de prendre $Y := AZ$, pour un $Z$ à définir (voir plus bas).

L'hypothèse essentielle est que $\mathbb{C}^n = AZ \oplus Y^\perp = Z \oplus (A^{-1}Y)^\perp$.

On se donne ensuite deux projecteurs, typiquement

$$P := I - AZ(Y^*AZ)^{-1}Y^*$$

et

$$Q = I - Z(Y^*AZ)^{-1}Y^*A.$$

$P$ est donc un projecteur orthogonal sur $Y^\perp$ parallèlement à $AZ$, et $Q$ est un projecteur orthogonal sur $(A^{-1}Y)^\perp$ parallèlement à $Z$. On a la relation $PA = AQ$.

On peut écrire $x \in \mathbb{C}^n$ solution du système linéaire de la façon suivante :

$$ 
\begin{aligned}
x &= (I - Q)x + Qx \\
&= Z(Y^*AZ)^{-1}Y^*Ax + Qx \\
&= Z(Y^*AZ)^{-1}Y^*b + Qx.
\end{aligned}
$$

En pratique, $Y^*AZ \in \mathbb{C}^{r \times r}$ est factorisable (par LU, QR ou SVD).

Si l'on calcule $\tilde{x}$ solution de $PA\tilde{x} = Pb$, alors pour $x$ solution du système linéaire on peut écrire

$$ 
\begin{aligned}
Qx &= A^{-1}AQx \\
&= A^{-1}PAx \\
&= A^{-1}Pb \\
&= A^{-1}PA\tilde{x} \\
&= A^{-1}AQ\tilde{x} = Q\tilde{x},
\end{aligned}
$$

et l'on a donc

$$x = Z(Y^*AZ)^{-1}Y^*b + Q\tilde{x}.$$

Résoudre le système linéaire déflaté

$$PA\tilde{x} = Pb$$

permet donc de résoudre le système linéaire d'origine.
Une propriété importante de P est que, pour $x = Zz$ avec $z \in \mathbb{C}^r$,

$$ 
\begin{aligned}
PAx &= PAZz \\
    &= (AZ - AZ(Y^*AZ)^{-1}Y^*AZ)z \\
    &= 0,
\end{aligned}
$$

nous nous en servirons pour tester la précision des opérateurs que nous aurons défini dans le contexte numérique d'arithmétique finie.

Il nous faut encore définir $Z$. L'idée dans notre cas, assez différente du principe exposé en début de chapitre, est de tirer parti du fait que les ondes propagatives se concentrent plus à l'intérieur des cellules que les ondes évanescentes, concentrées sur les bordes (voir [13, section 2.1 p. 4]), en conservant des espaces locaux les ondes les plus proches d'ondes planes propagatives du domaine global (non du maillage), en réalisant des produits scalaires entre ces ondes globales et celles des espaces locaux.

Soit une base de Trefftz $\{\phi_{j1 \le j \le J_K}^k\}_{1 \le k \le K}$ d'un maillage à $K$ cellules, de tailles locales respectives $J_K$. On génère $r$ ondes planes propagatives $\{\psi_s\}_{1 \le s \le r}$ sur le domaine global, et l'on construit chaque colonne $Z_s$ de $Z$, décomposée en superposition de vecteurs locaux $Z_s^k$ de tailles respectives $J_K$,

$$ Z_s = \begin{pmatrix} Z_s^1 \\ \vdots \\ Z_s^k \\ \vdots \\ Z_s^K \end{pmatrix}, $$

en résolvant, par SVD régularisée, le système linéaire suivant

$$ \left( \cdots \begin{pmatrix} \phi_j^k, \phi_l^k \end{pmatrix}_{\mathbb{H}_\omega^1(\Omega_k)} \cdots \right) \left( \begin{pmatrix} \vdots \\ (Z_s^k)_j \\ \vdots \end{pmatrix} \right) = \left( \cdots \begin{pmatrix} \psi_s, \phi_l^k \end{pmatrix}_{\mathbb{H}_\omega^1(\Omega_k)} \cdots \right), $$

où, sur la cellule $\Omega_k$, le produit scalaire $\mathbb{H}_\omega^1(\Omega_k)$ est défini, pour deux ondes planes $u, v$, par

$$ (u, v)_{\mathbb{H}_{\omega}^{1}(\Omega_{k})} := (u, v)_{\mathbb{L}^{2}(\Omega_{k})} + \omega^{-2}(\nabla u, \nabla v)_{\mathbb{L}^{2}(\Omega_{k})} $$

(où l'on a pris $\omega$ le nombre de l'onde de l'équation de Helmholtz considérée pour ne pas confondre $\kappa$ et $k$).

Par abus de langage nous appellerons indifférement espace de déflation soit l'espace de Trefftz composé d'ondes propagatives cré sur le domaine global, soit sa projection (c'est-à-dire l'image de $Z$).

## 10.2 Expériences numériques

Nous avons commencé par implémenter la méthode décrite ci-dessus, en suivant les étapes suivantes :

- On construit $Z$ par projection $\mathbb{H}_\omega^1(\Omega_k)$,

- On stocke $Y = AZ$,

- On calcule la svd régularisée de $Y^*AZ$,

- On implémente les produits matrice-vecteur pour $P$, $PA$ et $Q$,

- On résoud $PA\tilde{x} = Pb$ (par Gmres),

- et enfin on calcule $x = Z(Y^*AZ)^{-1}Y^*b + Q\tilde{x}$.
### 10.2.1 Premiers essais

Nous expérimentons d'abord la déflation sur le problème de l'écran du chapitre précédent. Nous sommes alors directement confrontés à plusieurs problèmes : d'abord, le résidu de gmres stagne systématiquement quand son paramètre restart est plus petit que le nombre d'itérations nécessaire à la convergence. Ensuite, si gmres semble converger vers une solution correcte, le nombre d'itérations nécessaire est de loin supérieur à la version sans déflation. Par exemple, pour $\kappa = 40$, $N=200$ et un espace de déflation de taille $r = 100$, sur le maillage avec écran à 70 cellules de la section précédente, il faut 2039 itérations pour la version avec déflation, contre 217 pour la version sans (nous les comparons avec les mêmes éléments A et b du système linéaire, construit sur le même espace de Trefftz). Il semble bien que l'information que nous avons ajoutée au système linéaire de va pas dans le sens espéré.

Nous regardons alors la norme de $PAZx$, qui devrait être nulle comme nous l'avons remarqué, pour un $x \in \mathbb{C}^T$ aléatoire (mais le même pour les différentes mesures ci-dessous). Or pour l'exemple cité, elle est de 0.0136 environ...

Nous souhaitons investiger les seuils de régularisation $\epsilon_D$ et $\epsilon_M$ des svd de $D$ et $M = Y^*AZ$ respectivement, et montrons ci-dessous les valeurs obtenues pour différents seuils, en gardant les paramètres décrits plus haut.

<table><thead><tr><th>ε<sub>D</sub></th><th>ε<sub>M</sub></th><th>||PAZx||<sub>2</sub></th><th>nb itérs gmres</th></tr></thead><tbody><tr><td>1e-10</td><td>1e-10</td><td>1.362e-2</td><td>2039</td></tr><tr><td>1e-8</td><td>1e-10</td><td>1.054e-3</td><td>1795</td></tr><tr><td>1e-6</td><td>1e-10</td><td>3.35-4</td><td>1492</td></tr><tr><td>1e-4</td><td>1e-10</td><td>5.57-5</td><td>937</td></tr><tr><td>1e-2</td><td>1e-10</td><td>2.51e-5</td><td>547</td></tr><tr><td>1e-10</td><td>1e-8</td><td>0.104282</td><td>2030</td></tr><tr><td>1e-10</td><td>1e-12</td><td>0.000395</td><td>2048</td></tr><tr><td>1e-10</td><td>1e-14</td><td>0.000115</td><td>2056</td></tr></tbody></table>

Nous nous trouvons dans une situation très étrange où le fait de définir un seuil plus haut pour les valeurs singluières de la svd régularisee de $D$, et donc d'en réaliser une approximation plus grossière, améliore la situation, en tout cas au sens où le nombre d'itérations décroît. Mais on est loin d'améliorer la performance par rapport à la version sans déflation. Et c'est au contraire en régularisant moins $M$ que la norme étudiée baisse, mais cette fois au détriment du nombre d'itérations...

Nous remarquons par ailleurs que la décroissance du résidu de Gmres procède par paliers, parfois très longs en nombre d'itérations, et toujours situés au-dessus de la norme de $PAZx$ (sauf pour les plus fortes régularisations de $D$), allant dans le sens d'une corrélation entre les deux phénomènes : nous prenons un $x$ aléatoire, $||PAZx||_2$ est donc toujours différent, mais $||PA||_2$ semble être en cause.

<figure><img src="image_64.png" /><figcaption>FIGURE 58 – Déflation essai 1, palier du résidu de gmres et $||PAZx||_2$, $\epsilon_D = \epsilon_M = 1e-10$.</figcaption></figure>
### 10.2.2 Second test avec écran - en scindant l'espace grossier

Nous nous disons qu'il est possible qu'en présence de l'écran au milieu du domaine, créer un seul espace grossier global au domaine n'est peut-être pas optimal. Nous décidons donc de partitionner le domaine en deux domaines rectangles, l'un correspondant aux ordonnées positives, l'autre aux négatives. Nous créons pour chacun d'eux un espace de Trefftz de déflation, toujours composé d'ondes propagatives, tous deux étant définis avec des paramètres identiques hormis leur domaine. L'idée est ensuite de sélectionner, pour chacune des cellules du maillage sur lequel on procède aux approximations, lequel de ces deux espaces on inclut dans son espace local d'approximation. On calcule le barycentre de chaque cellule, et on inclut dans son espace local d'approximation l'espace grossier correspondant à son ordonnée. En conséquence la construction de Z s'en trouve modifiée puisqu'il faut prendre garde à n'évaluer, pour chacun de ses s sous-colonnes $Z_s^K$ correspondant à la cellule K, que les produits scalaires avec les ondes planes de l'espace grossier correspondant. Un exemple est montré en figure 59.

<figure><img src="image_65.png" /><figcaption>FIGURE 59 – Déflation test 3, maillage partitionné en deux zones correspondant aux deux espaces grossiers créés, et structure de la matrice Z correspondante.</figcaption></figure>

Nous testons avec $\kappa = 10$, $N = 120$ et $r = 20$, et comparons les approximations avec et sans déflation réalisées avec le même espace de Trefftz. Malheureusement, le résultat n'est toujours pas celui espéré : si l'approximation avec déflation obtenue semble similaire à celle sans déflation, il faut 1318 itérations à gmres pour converger, contre 293 pour la version sans déflation...

### 10.2.3 Inclusion de l'espace grossier dans l'espace d'approximation

Nous doutons de la pertinence comme de l'implémentation de la projection que nous opérons. Nous décidons d'essayer de simplifier la situation, en incluant directement l'espace de déflation dans notre espace d'approximation. Cela se réalise d'une part en incluant les $r$ ondes propagatives globales dans chacun des espaces d'approximation locaux, et d'autre part en construisant $Z$ non par projection $H_\omega^1(\Omega_k)$ mais en le définissant avec des 1 sur la diagonale de chaque bloc de lignes
correspondant à une cellule : pour chaque cellule $1 \le k \le K$,

$$ (Z_s)^k = \begin{pmatrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 0 \end{pmatrix}, $$

ce qui est bien équivalent à une projection idéale.

Nous nous plaçons cette fois sur le maillage de base à huit cellules pour encore simplifier le contexte, et pour $\kappa = 40$, $N=300$, $r = 100$ et $\epsilon_D=1e-10$, nous comparons les mesures de $||PAZx||_2$ selon le mode de construction de $Z$: nous obtenons une valeur de l'ordre de 2.0e-3 pour la version projetée, et de 2.2e-5 pour la version sans projection avec espace de déflation inclus dans celui d'approximation. C'est un premier pas vers le mieux, mais malheureusement nos tests avec gmres n'ont pas montré pour autant d'amélioration de vitesse de convergence avec cette version non projetée.

### 10.2.4 Seconde implémentation

Un enseignement de ce stage ayant été l'importance de l'ordre dans lequel on effectue les opérations (notamment les produits matriciels) en contexte d'arithmétique finie, nous suspectons que notre implémentation de $Y^*AZ$ n'est pas correcte de ce point de vue, et que sa pseudo-inverse svd régularisée génère de ce fait une instabilité numérique trop importante. nous développons donc totalement l'expression de $Y^*AZ$, afin d'implémenter toutes les fonctions réalisant les différents produits matrice-vecteur intermédiaires, pour tâcher de minimiser l'accumulation d'erreurs.

On commence par écrire $AZ = (I - D^{-1}C)Z = Z - D^{-1}(CZ)$, puis

$$ \begin{aligned} A^*(AZ) &= (I - D_{\epsilon}^{-1}C^*)(Z - D_{\epsilon}^{-1}(CZ)) = (I - C^*D_{\epsilon}^{-*})(Z - D_{\epsilon}^{-1}(CZ)) \\ &= Z - C^*(D_{\epsilon}^{-*}Z) - D_{\epsilon}^{-1}(CZ) + C^*(D_{\epsilon}^{-*}(D_{\epsilon}^{-1}(CZ))) \\ &= Z - C^*(U_{\epsilon}(\Sigma_{\epsilon}^{-1}(V_{\epsilon}^*Z))) - V_{\epsilon}(\Sigma_{\epsilon}^{-1}(U_{\epsilon}^*(CZ))) + C^*(U_{\epsilon}(\Sigma_{\epsilon}^{-1}(V_{\epsilon}^*(V_{\epsilon}(\Sigma_{\epsilon}^{-1}(U_{\epsilon}^*(CZ)))))) \end{aligned} $$

où l'on a utilisé la factorisation SVD pseudo-inverse régularisée de $D$, et enfin

$$ \begin{aligned} Z^*(A^*(AZ)) &= Z^*Z - Z^*(C^*(U_{\epsilon}(\Sigma_{\epsilon}^{-1}(V_{\epsilon}^*Z)))) \\ &\quad - Z^*(V_{\epsilon}(\Sigma_{\epsilon}^{-1}(U_{\epsilon}^*(CZ)))) \\ &\quad + Z^*(C^*(U_{\epsilon}(\Sigma_{\epsilon}^{-2}(U_{\epsilon}^*(CZ))))), \end{aligned} $$

en utilisant le fait que $D_{\epsilon}^{-*}D_{\epsilon}^{-1} = U_{\epsilon}\Sigma_{\epsilon}^{-1}V_{\epsilon}^*V_{\epsilon}\Sigma_{\epsilon}^{-1}U_{\epsilon}^* = U_{\epsilon}\Sigma_{\epsilon}^{-2}U_{\epsilon}^*$ (nos tests ont plutôt montré une meilleure précision en faisant cette simplification).

Nous développons de même l'expression de $P$

$$ \begin{aligned} P &= I - AZ(Y^*AZ)^{-1}Y^* \\ &= I - (I - D^{-1}C)ZM_{\epsilon}^{-1}Y^* \\ &= I - ZM_{\epsilon}^{-1}Z^*A^* + D^{-1}CZM_{\epsilon}^{-1}Z^*A^*, \end{aligned} $$
où $M_{\epsilon}^{-1}$ est la svd pseudo-inverse régularisée de $Y^*AZ$. nous faisons de même pour $Q$. Malheureusement, cette nouvelle version n'a pas amélioré la vitesse de convergence, et prend toujours plus d'itérations que la version sans déflation pour converger.

Nous avons enfin souhaité investiguer les propriétés spectrales des opérateurs construits. Pour $\kappa = 40$, $N = 300$, $r = 100$, les seuils de svd à 1e-10, et sur le maillage de référence, nous comparons, pour les deux méthodes de construction de $Z$, nos deux méthodes de construction de $PA$ (notées $PA_1$ et $PA_2$ respectivement).

<table><thead><tr><th>Z</th><th>||PA<sub>1</sub>Zx||<sub>2</sub></th><th>||PA<sub>2</sub>Zx||<sub>2</sub></th><th>cond(A)</th><th>cond(PA<sub>1</sub>)</th><th>cond(PA<sub>2</sub>)</th></tr></thead><tbody><tr><td>version 1</td><td>0.00225</td><td>0.00225</td><td>3.73e6</td><td>3.32e16</td><td>8.25e9</td></tr><tr><td>version 2</td><td>9.31e-8</td><td>4.16e-6</td><td>3.73e6</td><td>6.92e17</td><td>2.44e12</td></tr></tbody></table>

Pour la première version de $Z$, les normes de $PA_1Zx$ et $PA_2Zx$ sont moins bonnes que celles obtenues avec la seconde version (au sens de plus éloignées de 0). Les conditionnements de $PA_1$ et $PA_2$ sont supérieurs à celui de $A$, ce qui fait sens puisque par construction ces opérateurs deviennent singuliers. Le fait que ceux de la seconde version de $Z$ soit plus grand paraît donc cohérent avec le résultat sur les normes.

Nous affichons également les valeurs propres de $A$, $PA_1$ et $PA_2$ pour les deux méthodes de construction de $Z$ en figure 60.

<figure><img src="image_66.png" /><figcaption>FIGURE 60 – Déflation, comparaison des valeurs propres de $A$, $PA_1$ et $PA_2$. Première version de construction de $Z$ à gauche, seconde version à droite.</figcaption></figure>

Pour les deux modes de construction de $Z$, les spectres de $PA_1$ et $PA_2$ semblent rigoureusement identiques. On voit bien les valeurs propres ramenées à l'origine sous l'action de $P$, quoique le fait que cela produise un nuage autour de l'origine semble nous indiquer que cette action n'est pas correcte. On observe aussi malheureusement que pour les deux modes de construction, le rayon spectral excède à présent 1, ce qui ne devrait pas être le cas sous l'action sous d'un projecteur et ne peut qu'empêcher la convergence des solveurs itératifs.

Nous avons encore pu vérifier que les spectres de nos opérateurs $P$ étaient bien ceux de projecteurs, mais notre stage touchant à sa fin, nous n'avons malheureusement pas pu mener à bien cette partie du travail, et comprendre si cette situation défavorable était due à une instabilité numérique, à une implémentation défectueuse ou à un mauvais choix d'espace de déflation.
# 11 Conclusion

La seule partie de notre travail non documentée dans ce rapport est la version distribuée (MPI) des algorithmes de construction et de résolution du système, sur laquelle nous avons passé plusieurs semaines et qui n'a pu aboutir faute d'avoir trouvé une implémentation de gmres fonctionnant en MPI avec les structures de données COO-blocs du type de nos matrices creuses. Nous aurions pu implémenter notre version de l'algorithme mais cela nous aurait éloigné de notre sujet.

Hormis les pistes de l'opérateur T et de la déflation sur lesquelles nous nous sommes penchés, d'autres prolongements ont été évoqués, comme l'utilisation de l'oversampling pour essayer de capturer des modes supplémentaires de la solution, le fait de doubler la précision de toutes nos variables pour améliorer l'indépendance linéaire des ondes évanescentes des bases de Trefftz, ou enfin l'idée de conduire une analyse de perturbation pour comprendre ce qui ne fonctionne pas dans notre méthode de déflation.
Références

[1] URL : https://team.inria.fr/alpines/fr/.

[2] URL : https://who.rocq.inria.fr/Emile.Parolin/.

[3] Grégoire ALLAIRE. *Analyse numérique et optimisation*. Editions de l'Ecole Polytechnique, 2012. ISBN : 978-2-7302-1255-7.

[4] Hélène BARUCQ et al. « Local strategies for improving the conditioning of the plane-wave Ultra-Weak Variational Formulation ». In : *Journal of Computational Physics* 441 (2021), p. 110449. ISSN : 0021-9991. DOI : https://doi.org/10.1016/j.jcp.2021.110449. URL : https://www.sciencedirect.com/science/article/pii/S0021999121003442.

[5] Olivier CESSENAT et Bruno DESPRÉS. « Application of an Ultra Weak Variational Formulation of Elliptic PDEs to the Two-Dimensional Helmholtz Problem ». In : *SIAM Journal on Numerical Analysis* 35 (1998), p. 255-299. URL : https://api.semanticscholar.org/CorpusID:120068011.

[6] Francis COLLINO, Souad GHANEMI et Patrick JOLY. « Domain decomposition method for harmonic wave propagation : a general presentation ». In : *Computer Methods in Applied Mechanics and Engineering* 184.2 (2000), p. 171-211. ISSN : 0045-7825. DOI : https://doi.org/10.1016/S0045-7825(99)00228-5. URL : https://www.sciencedirect.com/science/article/pii/S0045782599002285.

[7] David L. COLTON et Rainer KRESS. *Inverse acoustic and electromagnetic scattering theory*. 2013. URL : http://link.springer.com/book/10.1007%2F978-1-4614-4942-3.

[8] Victorita DOLEAN, Pierre JOLIVET et Frédéric NATAF. *An introduction to domain decomposition methods*. Algorithms, theory, and parallel implementation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015. ISBN : 978-1-611974-05-8. URL : http://dx.doi.org/10.1137/1.9781611974065.ch1.

[9] Nicola GALANTE, Andrea MOIOLA et Emile PAROLIN. *Stable approximation of Helmholtz solutions in the ball using evanescent plane waves*. 2024. arXiv : 2401.04016 [math.NA]. URL : https://arxiv.org/abs/2401.04016.

[10] Luis GARCIA RAMOS, René KEHL et Reinhard NABBEN. « Projections, Deflation, and Multigrid for Nonsymmetric Matrices ». In : *SIAM Journal on Matrix Analysis and Applications* 41.1 (2020), p. 83-105. DOI : 10.1137/18M1180268. eprint : https://doi.org/10.1137/18M1180268. URL : https://doi.org/10.1137/18M1180268.

[11] Ralf HIPTMAIR, Andrea MOIOLA et Ilaria PERUGIA. « Trefftz discontinuous Galerkin methods for acoustic scattering on locally refined meshes ». In : *Applied Numerical Mathematics* 79 (2014). Workshop on Numerical Electromagnetics and Industrial Applications (NELIA 2011), p. 79-91. ISSN : 0168-9274. DOI : https://doi.org/10.1016/j.apnum.2012.12.004. URL : https://www.sciencedirect.com/science/article/pii/S0168927413000044.

[12] HPC Wiki. URL : https://hpc-wiki.info/hpc/Scaling.

[13] PAROLIN, EMILE, HUYBRECHS, DAAN et MOIOLA, ANDREA. « Stable approximation of Helmholtz solutions in the disk by evanescent plane waves ». In : *ESAIM : M2AN* 57.6 (2023), p. 3499-3536. DOI : 10.1051/m2an/2023081. URL : https://doi.org/10.1051/m2an/2023081.
[14] A. Moiola R. HIPTMAIR et I. PERUGIA. « A survey of Trefftz methods for the Helmholtz equation. In Building bridges : connections and challenges in modern approaches to numerical partial differential equations, » in : *Lect. Notes Comput. Sci. Eng.* 114 (2016), p. 237-278.

[15] Yousef SAAD. *Iterative methods for sparse linear systems*. SIAM, 2003, p. I-XVIII, 1-528.
ISBN: 978-0-89871-534-7.