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

IMPLEMENTATION D'UNE METHODE MULTIPOLE
RAPIDE DANS LE CODE DE DIFFRACTION/RADIATION
AQUAPLUS POUR LA SIMULATION DE FERMES DE
RECUPERATEURS D'ENERGIE DES VAGUES

IMPLEMENTING A FAST MULTIPOLE ALGORITHM
INTO THE DIFFRACTION/RADIATION CODE
AQUAPLUS FOR SIMULATING ARRAYS OF WAVE
ENERGY CONVERTERS

B. BORGARINO*, A. BABARIT*, P. FERRANT*,

* Equipe Hydrodynqmique et Génie Océanique, Laboratoire de Mécanique des Fluides École Centrale de Nantes, BP 92101, 44321 NANTES Cedex 3
bruno.borgarino@ec-nantes.fr

Résumé

Cet article décrit l'implémentation d'un algorithme multipole rapide dans le code Aquaplus. L'objectif est d'accélérer la résolution des problèmes de diffraction/radiation, pour simuler par la méthode des singularités une ferme composée de nombreux récupérateurs d'énergie des vagues. Les fondements mathématiques nécessaires (expansion multipole de la fonction de Green), et les principes de l'algorithme sont détaillés ici. Les résultats obtenus par la méthode multipole sont confrontés à ceux de la méthode classique. L'accélération apportée par la méthode multipole est étudiée, en fonction des différents paramètres de l'algorithme.

Summary

This paper describes the implementation of a fast multipole algorithm (FMA) into the
software Aquaplus. The objective is to speed up the resolution of the diffraction/radiation
problems, to simulate arrays of wave energy converters by the boundary elements method
(BEM). The necessary mathematical formulations (multipole expansion of the free-surface
Green's function), are briefly detailed here, along with the working principle of the FMA.
The results obtained by the FMA are compared to the classical method results. The
acceleration brought by the FMA is investigated, depending on its parameters.
# I – Introduction

Les récupérateurs d'énergie des vagues (REV) doivent être à terme déployés en fermes, regroupant de 10 à 100 systèmes. De fortes interactions entre les systèmes peuvent avoir lieu, provoquées par les champs diffractés et rayonnés par les corps flottants en mouvement. Des interactions constructives sont possibles, mais le principal défi est de réduire les interactions destructrices, lorsque les systèmes se masquent mutuellement la houle incidente. Les contraintes opérationnelles nécessitant de regrouper les REV dans un espace limité, il est indispensable d'évaluer ces interactions, afin de déterminer la configuration optimale de la ferme.

Les méthodes des singularités, (Boundary Element Methods, BEM) permettent d'évaluer rapidement les mouvements d'un corps flottant. Elles nécessitent la construction et la résolution d'un système linéaire dont la matrice est pleine. Le temps de calcul évolue donc avec $N_{panneaux}^3$ (nombre de panneaux pour la discrétisation de la surface immergée des corps). Cette tendance constitue une goulet numérique qui empêche pour le moment l'étude d'une ferme de récupérateurs d'énergie des vagues.

La méthode des singularités peut être rendue plus rapide par un solveur itératif GMRes, lui-même accéléré par l'algorithme multipole rapide (Fast Multipole Algorithm, FMA). Cet algorithme différencie les interactions entre panneaux relativement proches et relativement lointains. Il permet d'éviter de construire une grande partie du système linéaire, d'où une économie importante de temps de calcul. Le FMA présente une complexité en $O(N_{panneaux})$.

Le FMA a été développé par Greengard et Rokhlin [6] pour le calcul des champs potentiels dans les systèmes à grands nombres de particules. Cet algorithme a été utilisé en hydrodynamique dans un code BEM temporel par Fochesato [5] pour étudier la propagation des vagues au-dessus d'un bathymétrie complexe. Utsunomyia et Okafuji ont développé en 2007 [11] l'expansion multipole de la fonction de Green à surface libre et en profondeur infinie pour résoudre le problème de radiation autour d'une très grande structure flottante. Utsunomyia et Wantanabe ont fait de même pour la fonction de Green en profondeur constante [12].

Cet article présente le rôle du FMA dans l'accélération du GMRes ainsi que ses principes de fonctionnement. Les résultats du code sont exposés; on détaillera l'influence de certains paramètres du FMA sur le temps de calcul.

# II – Méthodes

## II – 1 Méthode des singularités

Aquaplus résout le problème de diffraction/radiation dans le domaine fréquentiel [3]. L'eau est modélisée par un fluide non-visqueux et incompressible. La vitesse du fluide dérive d'un potentiel $\phi$. Le problème aux limites à résoudre est le suivant :

$$ \Delta\phi = 0 \text{ dans tout le domaine fluide} \quad (1) $$

$$ \frac{\partial\phi}{\partial n} = 0 \text{ sur le fond} \quad (2) $$

$$ \frac{\partial\phi}{\partial n} = \vec{V}_i \cdot \vec{n} \quad \text{sur la surface } S_i \text{ du corps } i \quad (3) $$

$$ \frac{\partial^2\phi}{\partial t^2} + g\frac{\partial\phi}{\partial z} = 0 \text{ sur la surface libre} \quad (4) $$
Avec g l'accélération de la pesanteur, $\vec{V}_i$ la vitesse sur la surface du corps i. Soit :

$$ \Phi = \mathrm{Re}[\phi e^{(-i\omega t)}] \qquad (5) $$

Le système à résoudre pour obtenir $\Phi$ est le suivant :

$$ \frac{\sigma(M)}{2} - \frac{1}{4\pi} \iint_{S_i} \sigma(M') \frac{\partial G(M, M')}{\partial n} dS = \vec{V}_i \cdot \vec{n}_{xyz} \qquad (6) $$

$$ \Phi(M) = \iint_{S_i} \sigma(M') \cdot G(M, M') dS \qquad (7) $$

Avec $\omega$ la pulsation d'onde, $\sigma(M)$ la densité de source, $G(M, M')$ la fonction de Green entre point influencé $M$ et point influençant $M'$. Si la surface totale des corps est discrétisée en $N_{panneaux}$ panneaux, l' équation 6 s'exprime sous la forme d'un système linéaire

$$ K \cdot \sigma = b \quad (K_{ij} = \frac{\partial G(M_i, M_j)}{\partial n}, \sigma_i = \sigma(M_i)) \qquad (8) $$

de taille $N_{panneaux} \times N_{panneaux}$. Nous décrivons ici les méthodes employées pour accélérer la résolution de ce système lors de simulations à grands $N_{panneaux}$ (à noter que l'équation 7 se résout elle aussi par l'utilisation du FMA, voir V – 3).

## II – 2 Rôle du FMA dans la méthode itérative

La résolution de l'équation 8 par le solveur itératif GMRes [10] fait intervenir des produits matrice-vecteur de type $K \cdot r^{(i)}$ ($r^{(i)}$ est le résidu à l'itération i). Ces opérations sont coûteuses en temps de calcul, principalement à cause de la construction de K lors du premier produit. K reste inchangée par la suite, car elle ne depend que de la longueur d'onde de la houle et de la position des panneaux (dans l'hypothèse des petites déplacements, cette position ne varie pas). La méthode multipole scinde le produit $K \cdot r^{(i)}$ :

$$ K \cdot r^{(i)} = (K_{lointain} \cdot r^{(i)})_{\mathrm{FMA}} + K_{proche} \cdot r^{(i)} \qquad (9) $$

Dans l'équation 9, $(K_{lointain} \cdot r^{(i)})_{\mathrm{FMA}}$ est directement évalué par le FMA sans construire la matrice $K_{lointain}$. Ce produit traite des interactions entre groupes de panneaux éloignés les uns des autres (il concerne donc la majeure partie des interactions). Cette évaluation est très rapide, même lorsque $N_{panneaux}$ est grand. En effet, supposons que l'on ajoute un panneau dans le problème. Ce panneau interagira avec le groupe de panneaux lointains (seulement une opération supplémentaire) et non avec l'ensemble des panneaux de ce groupe (autant d'opérations supplémentaires que de panneaux dans le groupe), d'où une complexité réduite à $O(N_{panneaux})$.

La matrice $K_{proche}$ qui concerne le reste des interactions est construite par la méthode classique.

## II – 3 Principe du FMA

[8] est une bonne introduction à l'emploi du FMA dans un code BEM. Les principes de cet algorithme sont décrits dans [6]. Le FMA décompose l'espace physique en cellules, par la construction d'un arbre hiérarchique : la cellule élémentaire qui englobe l'espace de simulation est scindée en huit cellules, (enfants), elles-mêmes décomposées en huit enfants, et ainsi de suite. Le nombre des décompositions $Nl$ correspond au nombre de niveaux de l'arbre hiérarchique. La différentiation des interactions entre panneaux proches
<figure><img src="image_2.png" /><figcaption>Figure 1 – Principe du FMA, illustré en 2D. $Nl = 4$ niveaux de decomposition</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2 – Moments de l'expansion selon la translation (repérée par le centre de l'expansion $G_{...}$ vers laquelle elle pointe)</figcaption></figure>

et lointains dépend des positions relatives des cellules qui contiennent ces panneaux. La figure 1 illustre en 2D les principes du FMA qui permettent d'exprimer l'influence d'un groupe de panneaux sur un autre. Les cellules Niv. 1 sont les enfants des cellules Niv. 0, et ainsi de suite.

* On somme l'influence des panneaux contenus dans chaque cellule (Mexp sur la figure 1). Cette étape fait intervenir le moment de l'expansion multipole (équation 16), multiplié par le résidu sur chacun des panneaux.

* On translate cette influence dans l'espace, par des opérateurs permettant de se déplacer le long de l'arbre hiérarchique (de centre de cellule à centre de cellule) : M2M et L2L pour les petites distances, M2L pour les grandes distances.

* Arrivée dans la cellule cible, l'influence du groupe de panneaux lointains est répartie sur les panneaux contenus dans cette cellule (voir l'équation 21).

# III – Expansion multipole de la fonction de Green

## III – 1 Fonction de Green en profondeur infinie

On considère le cas d'une ferme de REV en grande profondeur. On utilisera la fonction de Green à surface libre en profondeur infinie, dans le domaine fréquentiel [9] :

$$G(S, F) = \underbrace{\frac{1}{r} + \frac{1}{r_1}}_{G_1} + \underbrace{\int_0^\infty \frac{2\nu}{k-\nu} e^{k(z+\zeta)} J_0(kR) dk}_{G_2} \quad (10)$$

Avec :

* $(\xi, \eta, \zeta)$ les coordonnées du point source S

* $(x, y, z)$ les coordonnées du point champ F

* $\nu$ le nombre d'onde

* $J_n$ la fonction de Bessel de première espèce à l'ordre $n$

* $r = \sqrt{(x - \xi)^2 + (y - \eta)^2 + (z - \zeta)^2}$

* $r_1 = \sqrt{(x - \xi)^2 + (y - \eta)^2 + (z + \zeta)^2}$
$$ \bullet R = \sqrt{(x - \xi)^2 + (y - \eta)^2} $$

Dans le cas d'une profondeur finie ou d'une bathymétrie complexe, on considère que l'accélération substantielle apportée par le FMA permettra de représenter le fond marin par un corps supplémentaire [5], [11].

### III – 2 Expansion multipole

Les formulations présentées dans [11] ont été développées pour le cas d'un corps très allongé horizontalement. Le centre de l'expansion multipole appartenant à la surface libre, c'est probablement un algorithme FMA 2D qui a été mis en place. Nous présentons ici les formulations étendues, dans lesquelles ce centre peut être en tout point du domaine fluide. Les formules suivantes concernent la fonction de Green elle-même et non ses dérivées (voir pour cela [2]).

La fonction de Green est décomposée en un champ proche $G_1$ (qui décroît selon $\frac{1}{r}$) et un champ lointain $G_2$ (qui décroît selon $\frac{1}{\sqrt{r}}$). Les expansions multipoles de ces deux fonctions, en un centre situé à une profondeur $z_G$ quelconque, sont les suivantes :

$$ G_1 = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} r_{\xi}^{n} \frac{P_{n}^{\vert m \vert} (\cos \theta_{\xi})}{(n + \vert m \vert)!} e^{-im\alpha_{\xi}} e^{im\alpha'_{x}} (n - \vert m \vert)! \left( \frac{P_{n}^{\vert m \vert} \cos \theta_{x}}{r_{x}^{n+1}} (-1)^{n+m} + \frac{P_{n}^{\vert m \vert} \cos \theta_{x'}}{r_{x'}^{n+1}} \right) \quad (11) $$

$$ G_2 = 2 \sum_{n=0}^{\infty} \sum_{m=-n}^{n} r_{\xi}^{n} \frac{P_{n}^{\vert m \vert} (\cos \theta_{\xi})}{(n + \vert m \vert)!} e^{-im\alpha_{\xi}} G_3(z'_{x}, R'_{x}) e^{im\alpha'_{x}} \quad (12) $$

$G_3$ est défini ainsi [11] :

$$ G_3(z_x, R_x) = \int_0^\infty \frac{1}{k-\nu} e^{kz} J_{\vert m \vert} (k R_x) k^n dk \quad (13) $$

$$ G_3 = \begin{cases} \sum_{p=1}^{n-|m|} \nu^{p-1} (n-p-|m|)! \frac{P_{n-p}^{\vert m \vert} (\cos \theta_x)}{r_x^{n-p+1}} & + \nu^{n-|m|} \int_0^\infty \frac{1}{k-\nu} e^{kz} J_{\vert m \vert} (k R_x) k^{|m|} dk & \text{if } n > |m| \\ \int_0^\infty \frac{1}{k-\nu} e^{kz} J_{\vert m \vert} (k R_x) k^{|m|} dk & \text{if } n = |m| \end{cases} \quad (14) $$

La dernière intégrale s'évalue par l'équation (15), $\forall m > 0$ :

$$ \int_0^\infty \frac{1}{k-\nu} e^{kz} J_{\vert m \vert} (k R_x) k^{|m|} dk = \frac{2}{\pi} \left( \frac{-\nu^m \pi^2 e^{kz}}{4} [(-1)^m \mathbf{H}_{-\mathbf{m}}(\nu R_x) - Y_m(\nu R_x)] \right. \\ \left. + \frac{\sqrt{\pi} (2R_x)^m \Gamma(m+\frac{1}{2})}{2} \int_0^z \frac{e^{\nu(z-s)}}{(s^2 + R_x^2)^{m+\frac{1}{2}}} ds \right) - \pi \nu^m i e^{\nu z} H_m^{(2)}(\nu R_x) \quad (15) $$

avec **H** la fonction de Struve, $H^{(2)}$ la fonction de Hankel de deuxième espèce, $Y$ la fonction de Bessel de deuxième espèce, et $\Gamma$ la fonction Gamma. Les coordonnées de $F$ sont $(r'_x, \theta'_x, \alpha'_x)$ si on les exprime dans le système centrée en $G'(x_G, y_G, -z_G)$ (voir figure 3). Le terme relatif au point source, présent dans $G_1$ (équation 11) et $G_2$ (équation 12) est le moment de l'expansion multipole :

$$ M_n^m = r_{\xi}^n \frac{P_n^{\vert m \vert} (\cos \theta_{\xi})}{(n + \vert m \vert)!} e^{-im\alpha_{\xi}} \quad (16) $$
<figure><img src="image_4.png" /><figcaption>Figure 3 – Système de coordonnées (cas simplifié, lorsque S, G et F sont alignés)</figcaption></figure>

### III – 3 Opérateurs de translation

Les opérateurs de translation permettent de déplacer le centre de l'expansion, en calculant le moment au nouveau centre en fonction du moment à l'ancien centre. La figure 2 résume les notations des moments résultant des différentes translations.

* La M2M (Moment To Moment) déplace le moment de l'expansion multipole d'une cellule enfant vers son parent. ($r_{x_M}$, $\alpha_{x_M}$, $\theta_{x_M}$) est la position du centre d'expansion d'origine vu depuis le nouveaux centre. $\epsilon_m = (-1)^m$ si $m < 0$, 0 sinon.

$$ \tilde{M}_n^m = \sum_{t=0}^{n} \sum_{s=-t}^{t} M_{n-t}^{m-s} \epsilon_{m-s} \epsilon_s \epsilon_m r_{x_M}^t \frac{P_t^{|s|} (\cos \theta_{x_M})}{(t + |s|)!} e^{-is\alpha_{x_M}} \quad (17) $$

* La M2L (Moment To Local) transforme l'expansion multipole au centre de la cellule influençante en expansion de Taylor au centre de la cellule influencée. Les coordonnées du centre de l'expansion multipole et de son symétrique vus depuis le centre de l'expansion de Taylor sont indexées par $x_{ML}$ and $x'_{ML}$.

$$ L_t^{s(1)} = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \epsilon_{m-s} \epsilon_s \epsilon_m M_n^m [(n+t) - |s-m|]! \left( \frac{P_{n+t}^{|s-m|} (\cos \theta_{x_{ML}})}{r_{x_{ML}}^{n+t+1}} (-1)^{n-m} - \frac{P_{n+t}^{|s-m|} (\cos \theta_{x'_{ML}})}{r_{x'_{ML}}^{n+t+1}} \right) e^{-i(s-m)\alpha_{x_{ML}}} \quad (18) $$

$$ L_t^{s(2)} = \sum_{n=0}^{\infty} \sum_{m=-n}^{n} \epsilon_{m-s} \epsilon_s \epsilon_m M_n^m [(n+t) - |s-m|]! e^{-i(s-m)\alpha_{x_{ML}}} \int_0^{\infty} \frac{k}{k-\nu} e^{-k|z_{x_{ML}}|} J_{|s-m|}(kR_{x_{ML}}) k^{n+t} dk \quad (19) $$

* La L2L (Local To Local) déplace le moment de l'expansion de Taylor d'un parent vers un de ses enfants. ($r_{x_L}$, $\alpha_{x_L}$, $\theta_{x_L}$) est la position du centre d'origine de l'expansion vu depuis le nouveau centre.

$$ \tilde{L}_u^{v(1)(2)} = \sum_{t=v}^{\infty} \sum_{s=-t}^{t} \frac{\epsilon_{u-s} \epsilon_u}{\epsilon_s} L_t^s (-1)^{(t-v)+(u-s)} r_{x_L}^{t-v} \frac{P_{t-v}^{u-s} (\cos \theta_{x_L})}{[(t-v) - |u-s|]!} e^{-i(u-s)\alpha_{x_L}} \quad (20) $$

On obtient finalement l'évaluation de $G_1$ et $G_2$ par :

$$ G_{(1)(2)} = \sum_{v=0}^{\infty} \sum_{u=-v}^{v} \tilde{L}_{v}^{u(1)(2)}(r_{x})^{v} \frac{P_{v}^{|u|} (\cos \theta_{x})}{(v + |u|)!} e^{iu\alpha_{x}} \quad (21) $$
avec $(r_x, \alpha_x, \theta_x)$ les coordonnées du point champ vu depuis le centre de l'expansion de Taylor. Le second terme de l'équation 21 est similaire à l'équation 16 car il exprime le transfert de l'influence des panneaux distants (contenue dans l'expansion de Taylor) entre deux points rapprochés (centre de la cellule - panneau qui lui appartient).

# IV – Implementation d'un algorithme multipole

Le code Distributed Parallel Multipole Tree Algorithm (DPMTA) a été développé jusqu'en 2002 à l'université de Duke (USA). Cet outil open-source est conçu pour être implanté dans des codes BEM. Les formulations précédentes y ont été intégrées. Les communications entre DPMTA et Aquaplus à chaque produit matrice-vecteur dans GMRes sont les suivantes :

* le DPMTA recoit les positions des panneaux, la longueur d'onde de la houle et le résidu $r^{(i)}$

* le DPMTA évalue le produit $K_{lointain} \cdot r^{(i)}$

* le DPMTA génère la liste des interactions proches à calculer directement et l'envoie à Aquaplus, qui construit $K_{proche}$ (seulement à la première itération)

* les contributions des interactions proches et lointaines sont sommées dans le GMRes pour obtenir le produit complet à l'itération $i$

# V – Validation du logiciel

Le bon accord entre les résultats d'Aquaplus version classique et du code modifié a été vérifié à chaque étape du calcul : précision de l'expansion multipole, des produits matrice-vecteur, et de la résolution complète. On résume ici ces résultats en soulignant l'influence des différents paramètres de l'algorithme multipole.

## V – 1 Convergence de l'expansion multipole

La fonction de Green a été calculée par la méthode classique et par l'expansion multipole en faisant varier les distances $R$ et $Z$ entre point source et point champ, ainsi que le nombre d'onde $\nu$. Cette fonction a été tracée dans le domaine non-dimensionnel ($\nu.R$, $\nu.Z$). Afin d'étudier la convergence de l'expansion multipole, les trois paramètres suivant ont été balayés : $Np$ (ordre de troncature de l'expansion multipole), $\mathbf{r}_{G1} = \frac{r_\xi}{r_x}$ (joue sur la convergence du champ proche) et $\mathbf{r}_{G2} = \nu \cdot r_\xi$ (convergence du champ lointain). On relève les différences entre les valeurs calculées par la méthode classique et par l'expansion multipole sur l'ensemble du domaine (figure 4, figure 5). Les différences entre les méthodes sont petites lorsque :

* l'ordre de troncature $Np$ est grand

* le point source est relativement proche du centre de l'expansion($\mathbf{r}_{G1}$ petit) pour $G1$

* la distance point source - centre de l'expansion est petite devant la longueur d'onde pour $G2$

Sur $G2$ la convergence tends à devenir mauvaise de manière brutale quand $\mathbf{r}_{G2}$ augmente. Dans l'algorithme complet, les valeurs de $\mathbf{r}_{G1}$ dépendent directement du "multipole acceptance criterion" (mac, [1]) dans le DPMTA. Ce critère permet de différencier interactions proches et interactions lointaines entre les cellules (100% interactions directes - $0 < mac < 1$ - 100% interactions FMA). Diminuer le mac améliore donc la concordance entre les deux méthodes :
<figure><img src="image_5.png" /><figcaption>Figure 4 – Moyenne et maximum de l’er-
reur relative sur $G_1$, selon $N_p$ et $\mathbf{r}_{G_1}$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5 – Moyenne et maximum de l’er-
reur relative sur $Im(G_2)$, selon $N_p$ et $\mathbf{r}_{G_2}$</figcaption></figure>

• en améliorant la convergence de l'expansion multipole ($\mathbf{r}_{G_1}$ sera statistiquement réduit)

• en augmentant la proportion des interactions calculées par la méthode directe.

$\mathbf{r}_{G_2}$ dépend de la taille des cellules dans lequel le moment de l'expansion est calculé. L'arête $a$ d'une cellule se calcule en fonction de son niveau $n$ et de l'arête $a_{max}$ de la cellule de niveau 0 : $a = a_{max}/2^n$. L'utilisateur peut choisir manuellement les dimensions de cette cellule pour mieux maîtriser le critère $\mathbf{r}_{G_2}$.

## V-2 Produit matrice-vecteur

On évalue les différences sur un produit $K.r^{(i)}$ dans GMRes selon la méthode employée. Le corps étudié est un hémisphère décrit par 400 panneaux. Le *mac* est tel que le FMA prenne en compte autant d'interactions lointaines pour $Nl = 3$ et 4. Les résultats apparaissent sur la figure 6. La précision sur $G_2$ est plus sensible à $Nl$ que sur $G_1$ ; en effet un petit $Nl$ (grands $r_\xi$) nécessite un ordre $Np$ plus élevé pour que l'expansion de $G_2$ converge. $\mathbf{r}_{G_1}$ n'est pas affecté par le nombre de niveaux, d'où une bonne convergence de $G_1$ dès un faible $Np$.

## V-3 Coefficients de radiation d'un hémisphère

Les problèmes de radiation en pilonnement et cavalement d'un hémisphère sont résolus en utilisant le FMA. Les résultats sont comparés à la ceux de la méthode de Gauss (solveur d'origine dans Aquaplus) et aux solutions analytiques de [7]. La période de houle (de 3 à 20 s) a été testée sur des sphères de rayons $a = 5$ m, 10 m et 30 m. Le calcul implique la résolution de $K.\sigma = \frac{\partial\phi}{\partial n}$ pour obtenir la distribution de sources, puis le produit $S.\sigma = \Phi$ pour obtenir le potentiel (équation 7), lui aussi calculé par le FMA. La figure 7 montre un bon accord entre les différentes solutions. Cependant, des calculs sur des plus grandes distances présentent des difficultés de convergence de $G_2$ (quand la distance panneau - centre de l'expansion multipole est grande devant la longueur d'onde). Une
<figure><img src="image_7.png" /><figcaption>**Figure 6** – Errur relative due aux FMA selon Np et Nl (ligne continue : maximum ; pointillés : moyenne ; carrés : Im(G₂), Nl = 3, mac = 0.9 ; triangles : Im(G₂), Nl = 4, mac = 0.45 ; ∇ : G₁, Nl = 3, mac = 0.9 ; losange : G₁, Nl = 4, mac = 0.45)</figcaption></figure>

solution consistant à placer les centres des cellules au barycentre des panneaux qu'elles
contiennent est à l'étude.

VI – Performance du code

L'accélération apportée par la méthode multipole a été évaluée pour la résolution du problème de radiation en pilonnement d'un hémisphère, représenté par un nombre de panneaux variables (de 100 à 9800). Deux résultats de temps CPU sont mis en valeur ici :

* le temps nécessaire pour effectuer le premier produit matrice-vecteur du GMRes

* le temps nécessaire pour obtenir la valeur du potentiel en chacun des panneaux (résolution complète)

Le nombre de niveaux Nl et le mac ont été fixés de manière à ce qu'environ 30% des interactions soient calculées par la méthode directe.

VI – 1 Premier produit matrice-vecteur

La figure 8 compare les temps CPU du 1er produit matrice-vecteur, calculés dans le GMRes de manière classique et par la méthode multipole. Le temps CPU du FMA est tracé pour Nl = 3, 4, 5. L'évolution du temps CPU en fonction de Npanneaux montre une rupture de pente, qui s'explique ainsi :

* Pour N<sub>panneaux</sub> faible, peu de cellules de l'arbre hiérarchique sont occupées. Augmenter N<sub>panneaux</sub> revient à peupler cet arbre, d'où une croissance rapide du nombre de M2L, étape la plus consommatrice de temps CPU [4]. Le temps CPU croît avec une forte pente.

* Une fois la majorité des cellules peuplée d'au moins un panneau, l'évolution du temps CPU suit une pente plus douce. L'ajout d'un panneau dans une cellule occupée n'entraîne qu'une expansion multipole supplémentaire.

On remarque cependant que la courbe Nl = 3 pour un grand Npanneaux suit une tendance assez similaire à celle du GMRes. Le faible nombre de cellules se traduit par une
<figure><img src="image_8.png" /><figcaption>Figure 7 – Amortissement *B* et masse ajoutée *A*. Pointillés : solutions analytiques ; triangles : Gauss ; carrés : GMRes + FMM. *Np* = 15, *Nl* = 4, *mac* = 1.0 pour le cavalement (1) ; *Np* = 10, *Nl* = 4, *mac* = 0.75 pour le pilonnement (0)</figcaption></figure>

importante proportion d'interactions directes, d'où une allure en $O(N_{panneaux}^2)$.

## VI – 2 Résolution complète

La figure 10 présente les temps CPU pour une résolution complète. Si l'on retrouve les mêmes tendances qu'auparavant, en revanche le nombre de panneaux à partir duquel le FMA accélère le GMRes est considérablement plus grand. En effet le temps employé par le FMA pour évaluer $(K_{lointain} \cdot r^{(i)})_{FMA}$ est supérieur au temps consacré au produit $K \cdot r^{(i)}$ (du moment que $K$ a été construite auparavant). Par conséquent :

* au premier produit $K \cdot r^{(i)}$, le FMA accélère nettement le GMRes en évitant la construction de $K$

* les produits suivants sont par contre ralentis par l'emploi du FMA

Faire passer $N_p$ de 10 à 5 améliore considérablement l'efficacité du FMA ([1] recommande un ordre de troncature compris entre 4 et 8). Selon les paramètres, il est intéressant d'utiliser le FMA à partir de 3000 - 6000 panneaux. La tendance du temps de calcul suggère fortement que l'intérêt du FMA deviendra évident au delà de 10000 panneaux.

La figure 9 montre le temps CPU cumulé à chaque étape du FMA : upward pass (multipole expansion + M2M), downward pass (M2L+ L2L), calcul direct. Le downward pass est dans tout les cas l'étape la plus consommatrice en temps de calcul, notamment quand $N_l$ est grand. Pour notre cas d'étude ($N_{panneaux} \approx 10^4$) il n'est pas pertinent de travailler avec $N_l > 4$.

## VII – Conclusion

Cet article présente l'implémentation d'une méthode multipole rapide dans Aquaplus. Cette méthode accélère les produits matrice-vecteur dans le solveur itératif GMRes. L'expansion multipole de la fonction de Green en profondeur infinie a été étendue pour être utilisée dans un algorithme 3D. Ces formulations présentent des difficultés de convergence lorsque le domaine de calcul est grand devant la longueur d'onde de la houle. Les trois solveurs disponibles doivent pour l'instant s'utiliser de manière complémentaire :
<figure><img src="image_9.png" /><figcaption>Figure 8 – Temps CPU pour le premier produit matrice-vecteur du GMRes. Trait continu : $N_p = 5$, pointillés : $N_p = 10$</figcaption></figure>

• Gauss : faible nombre de panneaux, toutes longueurs d'onde, pour permettre la validation et le calibrage des paramètres du FMA.

• Gmres + FMA : grand nombre de panneaux, grandes longueurs d'ondes

• GMRes seul : en complément du FMA (grand nombre de panneaux, petites longueurs d'ondes)

La combinaison de ces trois outils rend possible l'étude des fermes de récupérateurs d'énergie des vagues, sans négliger les interactions entre les systèmes. La manipulation du FMA reste délicate et seuls quelques jeux de paramètres ont prouvé leur efficacité. Les résultats présentées ici se placent à la frontière inférieure des possibilités de l'algorithme multipole. Par conséquent l'accélération qu'il apporte sera plus évidente lors de simulations plus importantes.

## Remerciements

Les auteurs remercient l' Agence de l'Environnement et de la Maîtrise de l'Energie (ADEME) et la Région des Pays de la Loire pour le financement de ces travaux de thèse.

## Références

[1] J. Board, Z. Hakura, W. Elliot, and W. Rankin. Scalable variants of multipole-accelerated algorithms for molecular dynamics applications. Technical report, Duke University, Department of Electrical Engineering, 1995.

[2] B. Borgarino, A. Babarit, and P. Ferrant. Extension of the free-surface green's function multipole expansion for the infinite water depth case. In *Proc of the 20th ISOPE Conference, Beijing (China)*, 2010.

[3] G. Delhommeau. Seakeeping codes aquadyn and aquaplus. In *19th WEGMENT School, Numerical Simulation of Hydrodynamics : Ship and Offshore Structures*, 1993.
<figure><img src="image_10.png" /><figcaption>Figure 9 – Temps CPU cumulé, résolution du pilonnement. Carrés : upward pass, cercles : downward pass, triangles : direct. $N_p = 5$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10 – Temps CPU pour la résolution du pilonnement. Trait continu : $N_p = 5$, pointillés : $N_p = 10$</figcaption></figure>

[4] W. Elliot and J. Board. Fast fourier transform accelerated fast multipole algorithm. *SIAM Journal on Scientific and Statistical Computing*, 17 :398–415, 1996.

[5] C. Fochesato. Modèles numériques pour les vagues et les ondes internes. PhD thesis, Ecole Normale Supérieure de Cachan, 2004.

[6] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. *Journal of Computational Physics*, 135 :280–292, 1997.

[7] A. Hulme. The wave forces acting on a floating hemisphere undergoing forced periodic oscillations. *Journal of Fluid Mechanics*, 121 :443–463, 1982.

[8] Y. J. Liu and N. Nishimura. The fast multipole boundary element method for potential problems : A tutorial. *Engineering Analysis with Boundary Elements*, 30 :371–381, 2006.

[9] J. N. Newman. Algorithms for free-surface green's function. *Journal of engineering mathematics*, 19 :57–67, 1985.

[10] Youcef Saad and Martin H. Schultz. Gmres : a generalized minimal residual algorithm for solving nonsynthetic linear systems. *SIAM Journal on Scientific and Statistical Computing*, 7, 1986.

[11] T. Utsunomiya and T. Okafuji. Wave response of a vlfs by accelerated green's function method in infinite water depth. *International Journal of Offshore and Polar Engineers*, 17 :30–38, 2007.

[12] T Utsunomiya and E Watanabe. Accelerated higher order boundary element method for wave diffraction/radiation problems and its applications. In *Proceedings of the 12th ISOPE Conference, Kitakyushu, Japan*, 2002.