# La méthode des éléments finis et les ordinateurs parallèles SIMD

Olivier Pironneau

► To cite this version:

Olivier Pironneau. La méthode des éléments finis et les ordinateurs parallèles SIMD. [Rapport de recherche] IRIA-RR-361, IRIA. 1979, pp.28. hal-04716622

HAL Id: hal-04716622

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

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

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

LA MÉTHODE
DES ÉLÉMENTS FINIS
ET LES ORDINATEURS
PARALLÈLES SIMD

Olivier PIRONNEAU

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

Rapport de Recherche N° 361

Septembre 1979
LA MÉTHODE DES ÉLÉMENTS FINIS ET LES ORDINATEURS
PARALLELES SIMD
Olivier Pironneau

Résumé :

Cette note de travail, préparée dans le cadre d'un groupe de travail sur le calcul parallèle au Laboria, a pour objectif de démontrer que les méthodes d'éléments finis sont aussi bien parallélisables que les méthodes de différences finies. On se place dans le cas d'un grand nombre de processus élémen-taires et une architecture SIMD. On étudie trois équations (équation de la chaleur, problème de Dirichlet, équations de Navier-Stokes) on étudie l'implémentation et le nombre d'opérations néces-saires et on compare éventuellement aux méthodes de différences finies. On montre que si P est le nombre de processeurs les temps calculs sont divisés par P ou P/2 pour les méthodes explicites. Pour les méthodes implicites, le temps calcul est conditionné par la résolution des systèmes linéai-res.

Abstract :

This Laboratory report is connected with the workshop on parallel computation at Laboria.
Its purpose is to show that finite element methods are as easy to parallelize as the finite different
method.

We assume that the computer is made of a large number of elementary processors in a SIMD archi-
tecture. We study 3 equations : the heat equation, the Dirichlet problem and the Navier-Stokes
equations ; we investigate a possible implementation with super-elements, we count the number of
operations and we compare with the similar finite difference methods.

It is found that the computing time is divided by a factor proportional to the number of proces-
sors. However for implicit methods the computing time is mostly due to the solution of linear
systems thus it would be desirable to modify the architecture so as to be able to solve properly
such systems.
## Introduction :

Dans un groupe de travail sur la calcul parallèle, FRAILONG [1] a proposé une architecture adaptée aux résolutions d'équations aux dérivées partielles par différences finies explicites et le problème s'est posé de savoir si cette architecture pouvait être utilisée avec autant de satisfaction sur des discrétisations par éléments finis, seules méthodes permettant de prendre en compte des domaines compliqués pour des équations paraboliques ou elliptiques.

L'architecture considérée est un réseau de processus élémentaires (microprocesseur). Chaque processeur est relié à ses quatre voisins Nord, Sud, Est, Ouest (figure 2) et le temps de transfert d'une information d'un processeur à son voisin est à peu près du même ordre de grandeur que le temps de calcul d'une addition ou d'une multiplication. La grande restriction impo-sée par ce type d'architecture est donc de rendre le transfert des informations d'un point quelconque du réseau à un autre extrêmement couteux.

On envisage environ 10.000 processeurs élémentaires et le temps d'une opé-ration élémentaire serait de 10 µ seconde. Le problème est de savoir si on peut résoudre les équations de Navier-Stokes en un temps raisonnable avec 50.000 points de discrétisation sur cette machine. Cette étude montre qu'avec une telle machine il suffit de 1 mn 36 " pour 60 itérations de temps d'un schéma implicite en temps, en éléments finis en espace $P^2/P^1$ avec 90.000 noeuds. En plus cette étude montre que la partie la plus coûteuse des calculs est la résolution des systèmes linéaires (et non la construction des matrices). Malheureusement, il n'existe pas de bons schémas explicites (sans résolution de système linéaire) pour les équations de Navier-Stokes ; en revanche si on rajoute des connections supplémentaires aux réseaux il serait possible peut être d'utiliser la méthode d'élimination de Gauss en parallèle ce qui diviserait par 10 les temps calculs.

On se place dans le cadre des problèmes bidimensionnels, mais tout s'applique aux problèmes tridimensionnels à condition de rajouter 2 liaisons à chaque micro-processeur et d'avoir un réseau tridimensionnel (ce qui est possible).
Enfin les comptes d'opérations sont nécessairement approximatifs car
sur des méthodes aussi complexes il est très difficile de ne pas en oublier ;
cette étude sera complétée par des simulations réelles pour vérifier les chif-
fres, par exemple sur la machine DAP de ICL dont la structure est identique
au nombre de processeurs près et au temps de calcul près à celle étudiée dans
cette note.
# 1. Résolution de l'équation de la chaleur par des méthodes explicites.

## 1.1. Le problème

Nous considérons le problème d'évolution de la température $\varphi$ dans un ouvert $\Omega$ de $\mathbb{R}^2$ pendant un intervalle de temps $[0,T[$ :

$$ (1.1) \quad \left\{ \begin{array}{l} \frac{\partial \varphi}{\partial t} - \Delta \varphi = f \quad \text{dans } \Omega \times ]0,T[ \\ \varphi(t=0) = 0 \quad \text{dans } \Omega \ ; \quad \varphi = 0 \quad \text{sur le bord } \Gamma \quad \text{de } \Omega \text{ pour tout } t \text{ de } ]0,T[. \end{array} \right. $$

## 1.2. Discrétisation par différences finies

Si $\Omega$ est un rectangle $]0,a[ \times ]a,b[$ on peut discrétiser (1.1) en espace par différences finies centrées ce qui donne :

$$ (1.2) \quad \left\{ \begin{array}{ll} \displaystyle \frac{\partial \varphi_{i,j}}{\partial t} - \frac{1}{h_1^2} (\varphi_{i+1,j} - 2\varphi_{i,j} + \varphi_{i-1,j}) - \frac{1}{h_2^2} (\varphi_{i,j+1} - 2\varphi_{i,j} + \varphi_{i,j-1}) = f_{i,j} & 1 \le i,j \le N-1 \ ; \\ \\ \varphi_{i,j}(t=0) = 0 \ ; \quad \varphi_{0,j} = \varphi_{N,j} = \varphi_{i,0} = \varphi_{N,j} = 0 \quad 0 \le i,j \le N \ . & \end{array} \right. $$

où

$$ \varphi_{i,j} = \varphi\left(\frac{i}{N}a, \frac{j}{N}b\right), \quad h_1 = a/N, \quad h_2 = b/N. $$

Avec des notations matricielles ce système s'écrit :

$$ (1.3) \quad \frac{\partial \Phi}{\partial t} - A \Phi = F \quad \forall t \in ]0,T[ \ ; \quad \Phi(t=0) = 0 \ ; $$

où $\Phi = (\varphi_{1,1}, \varphi_{1,2}, \dots, \varphi_{1,N-1}, \varphi_{2,1}, \dots, \varphi_{2,N-1}, \dots, \varphi_{N-1,1}, \dots, \varphi_{N-1,N-1})^t$ où F est formé à partir des $f_{i,j}$ de la même manière et où A est une matrice symétrique $(N-1)^2 \times (N-1)^2$ représentée sur la figure 1.

<figure><img src="image_4.png" /><figcaption>**Figure 1** : la matrice du système (1.3) provenant d'une discrétisation par différences finies.</figcaption></figure>
Le schéma de discrétisation en temps explicite le plus simple est le suivant :

$$ (1.4) \qquad \bar{\Phi}^{n+1} = \Phi^n + \Delta t (A \Phi^n + F^n) $$

où

$$ (1.5) \qquad \Phi^n = \Phi(n\Delta t), \quad F^n = F(n\Delta t). $$

Supposons qu'on dispose d'un ordinateur avec une unité centrale de P unités de calcul (processeur) reliées entre elles par un réseau rectangulaire Nord-Sud - Est-Ouest (cf. figure 2).

<figure><img src="image_5.png" /><figcaption>**Figure 2** : Réseau de processeurs avec communication Nord, Sud, Est, Ouest. Chaque processeur a une mémoire individuelle ; en un cycle d'horloge on ne peut lire que les informations contenues dans sa mémoire et les registres de ses voisins.</figcaption></figure>

Si

$$ (1.6) \qquad P \ge (N-1)^2 $$

on peut associer à chaque processeur un noeud du maillage ce qui rend le calcul de $\Phi^{n+1}$ à partir de $\Phi^n$ selon (1.4) particulièrement facile; En effet le réseau de communications permet de calculer

$$ \frac{1}{h_1^2}(\varphi_{1+i,j} - 2\varphi_{i,j} + \varphi_{1-i,j}) + \frac{1}{h_2^2}(\varphi_{1,j+1} - 2\varphi_{i,j} + \varphi_{1,j-1}) + f_{i,j} $$

au noeud $(i,j)$ pourvu que chaque $\varphi_{i,j}$ et $f_{i,j}$ soient stockés dans les processeurs correspondants.
Coût de l'itération :

<table><tr><td>11 additions / soustractions</td><td rowspan="3">par processeur et indépendamment</td></tr><tr><td>3 multiplications</td></tr><tr><td>4 transferts d'information</td></tr><tr><td colspan="2">des autres.</td></tr></table>

Soit au total 18 opérations élémentaires par itération de temps. Pour résoudre le problème entièrement il faudra à peu près $18 \times (\frac{T}{\Delta t}) \times \mu$ secondes si $\mu$ est le temps moyen en seconde du cycle de calcul de chaque processeur.

S'il y a plus de noeuds dans le réseau que de processeurs alors il faudra associer à chaque processeur non plus un noeud mais plusieurs noeuds voisins. Soit $NI \times P$ le nombre de noeuds dans le réseau différences finies, il faudra alors $18 \times NI \times (\frac{T}{\Delta t}) \times \mu$ secondes pour résoudre le problème.

Pour que ce schéma explicite soit stable il faut que

$$ (1.7) \qquad \Delta t \le \frac{1}{2} \min \{h_1^2, h_2^2\} $$

Le temps minimum de résolution du problème est donc :

$$ (1.8) \qquad CPU_t = 18 \times NI \times N^2 \times T \times \mu / \min \{a^2, b^2\} $$

**Exemple :**

$$ \begin{array}{l} a = b = T = 1 \ ; \ P = 10.000 \ ; \ \mu = 10^{-5} \ ; \ NI = 9 \ ; \ N = 300 \\ CPU_t = 24 \text{ mn.} \end{array} $$

On notera que ce temps calcul est surtout dû à la condition de stabilité (1.7) qui implique $O(N^2)$ pas de temps.

### 1.3. Discrétisation par éléments finis

La même méthode peut s'appliquer avec des éléments finis. On part de la formulation variationnelle du problème (1.1) dans laquelle $\frac{\partial \varphi}{\partial t}$ a été discrétisé selon le même schéma explicite :

$$ (1.9) \qquad \frac{1}{\Delta t} (\varphi^{n+1} - \varphi^{n, w}) + (\nabla \varphi^{n}, \nabla w) = (f^{n}, w) \quad \forall w \in H_0^1(\Omega) $$

où ( , ) désigne le produit scalaire de $L^2(\Omega)^m$, $m = 1$ ou $2$ :

$$ (1.10) \qquad (u,v) = \int_{\Omega} u(x)v(x)dx $$
où

$$
(1.11) \quad H_0^1(\Omega) = \{ w \in L^2(\Omega) : \nabla w \in (L^2(\Omega))^2 , w = 0 \text{ sur } \Gamma \}
$$

et où comme ci-dessus :

$$
\varphi^n(x) = \varphi(n\Delta t, x)
$$

Dans les méthodes d'éléments finis on discrétise $H_0^1(\Omega)$.
La méthode
la plus simple consiste à recouvrir $\Omega$ avec des triangles $T \in \mathcal{J}_h$ avec les pro-
priétés suivantes :

. T ∩ T' = ∅ ou 1 côté ou 1 sommet ∀T, T' ∈ J_h .

. h = longueur du plus grand des côtés des T de J_h.

. distance entre la frontière Γ_h de Ω_h = UT et Γ de l'ordre de h² et
tout sommet de Γ_h est sur Γ.

Alors on démontre que

(1.12) $V_{oh} = \{w_h \in C^0(\bar{\Omega}_h) : w_h|_T \text{ affine } \forall T \in \mathcal{J}_h ; w_h|_{\Gamma_h} = 0\}$ est une approximation de $H_0^1(\Omega)$ ; $C^0$ désigne l'espace des fonctions continues.

$V_{oh}$ est un espace de _dimension finie_ qui a pour base par exemple les fonctions
w^i qui valent 1 au sommet numéro i et zéro sur les autres qui sont _continues_ et _affines par morceaux_ :

$$
\{w^i\}_{i \in I} ; w^i \in V_{oh} ; w_i(q^j) = \delta_{ij} \quad \forall i,j \in I
$$

où {$q^j$} désigne l'ensemble des sommets des triangles de $\mathcal{J}_h$, I l'ensemble des indices des sommets qui ne sont pas sur $\Gamma$ et $\delta_{ij}$ le symbole de Kronecker.

On remplace donc (1.9) par

$$
\left\{
\begin{array}{ll}
\frac{1}{\Delta t} (\varphi_h^{n+1} - \varphi_h^n, w_h) + (\nabla \varphi_h^n, \nabla w_h) = (f^n, w_h) & \forall w_h \in V_{oh} \\
\\
\varphi_h^{n+1} \in V_{oh} ; &
\end{array}
\right.
$$

ce qui est équivalent à
<figure><img src="image_6.png" /></figure>

**Figure 3 :** Eléments non nuls (en foncés) de la bande supérieure d'une matrice provenant d'une discrétisation quelconque par éléments finis. (fournie par J. PERIAUX)

On notera que la structure est quasi-régulière seulement.

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

**Figure 4 :** Un maillage régulier d'éléments finis.
$$
(1.13) \quad \sum_{i \in I} \frac{1}{\Delta t} \varphi_i^{n+1} (w^i, w^j) = (f, w^j) + \frac{1}{\Delta t} (\varphi_h^n, w^j) - (\nabla \varphi_h^n, w^j) \quad \forall j \in I
$$

où $\varphi_i = \varphi_h(q^i)$. On démontre qu'on ne détériore pas la précision de la méthode
si on utilise l'approximation dite du "mass limping" :

$$
(w^i, w^j) = \delta_{ij} \sigma_i ; \quad \sigma_i = \frac{1}{6} \sum_{\substack{T \in \mathcal{J} \\ T \supset \{q^i\}}} a(T) ; \quad a(T) = \text{aire de } T.
$$

Dans ces conditions (1.13) devient :

$$
(1.14) \quad \varphi_j^{n+1} = \Delta t f_j^n + \varphi_j^n - (\nabla \varphi_h^n, \nabla w^j) \Delta t / \sigma_j \quad \forall j \in I
$$

Soit encore

$$
\bar{\Phi}^{n+1} = \Delta t F + \Phi^n - \Delta t A\Phi^n
$$

avec

$$
a_{ij} = (\nabla w^i, \nabla w^j) / \sigma_j
$$

La figure 3 montre un exemple d'une telle matrice A construite à
partir d'une triangulation quelconque numérotée par un algorithme de minimi-
sation de largeur de bande.

Rappelons brièvement que sur un triangle T de sommets q¹, q², q³,
si φ et ψ sont deux fonctions affines valant φᵢ, ψᵢ sur les sommets qᵢ alors

$$
(1.15) \quad (\varphi, \psi)_T = \int_T \varphi \psi \, dx = \left[ (\varphi_1 + \varphi_2)(\psi_1 + \psi_2) + (\varphi_1 + \varphi_3)(\psi_1 + \psi_3) + (\varphi_2 + \psi_3)(\psi_2 + \psi_3) \right] \frac{a(T)}{12}
$$

$$
(1.16) \quad a(T) = \frac{1}{2} \left[ (q_1^1 - q_1^3)(q_2^2 - q_2^3) - (q_1^2 - q_1^3)(q_2^1 - q_2^3) \right]
$$

$$
(1.17) (\nabla \varphi, \nabla \psi)_T = \int_T \nabla \varphi \cdot \nabla \psi \, dx = \left( \begin{vmatrix} \varphi_1 q_2^1 & 1 \\ \varphi_2 q_2^2 & 1 \\ \varphi_3 q_2^3 & 1 \end{vmatrix} \times \begin{vmatrix} \psi_1 q_2^1 & 1 \\ \psi_2 q_2^2 & 1 \\ \psi_3 q_2^3 & 1 \end{vmatrix} + \begin{vmatrix} q_1 \varphi_1 & 1 \\ q_1 \varphi_2 & 1 \\ q_1 \varphi_3 & 1 \end{vmatrix} \times \begin{vmatrix} q_1 \psi_1 & 1 \\ q_1 \psi_2 & 1 \\ q_1 \psi_3 & 1 \end{vmatrix} \right) / 4 a(T)
$$
Pour évaluer la méthode plaçons nous d'abord dans les conditions optimales : on suppose qu'il y a un processeur par sommet et que $\Delta t, f_i, \varphi_i$ ($\nabla w^i; \nabla w^j$) $\forall j$ tel que $q^j$ soit voisin de $q^i$, ont été stockés dans la mémoire locale du processeur i. Nous verrons au paragraphe suivant que les coûts de fabrication des $a_{ij}$ sont négligeables.

On suppose enfin que le maillage est régulier c'est-à-dire engendré par 3 familles de courbes quasiparallèles (cf. figure 4) ; on peut alors associer à chaque processeur un noeud du maillage de telle sorte que si deux noeuds sont voisins les 2 processeurs correspondants sont aussi voisins. Les communications se font donc par :

1 transfert direct pour les noeuds voisins nord, sud, est, ouest.

1 transfert indirect représenté par 2 transferts directs pour les noeuds nord-ouest et sud-est.

Pour calculer $(\nabla \varphi_h^n, \nabla w^i)$ il faut :
7 multiplications, 6 additions et 8 transferts directs.

Pour calculer $\varphi_i^{n+1}$ il faut donc
25 cycles de calcul par itération de temps et le calcul de chaque $\varphi_i^{n+1}$ se fait en parallèle sur les i. Au total $CPU_t = 25 \times N^2 \times T \times \mu/\text{min}\{a^2, b^2\}$ s'il y a assez de processeurs.

Avec les chiffres de l'exemple ci-dessus on obtiendrait $CPU_t = 32 \text{ mn}$ chiffre tout à fait comparable à celui obtenu pour les différences finies.
Même pour les systèmes à matrices variables :

$$ (1.18) \quad \frac{\partial \varphi}{\partial t} - \nabla \cdot (\rho(x,t) \nabla \varphi) = f $$

les comparaisons restent honorables :

Différences finies : 24 cycles

Eléments finis : 60 cycles.

Les difficultés en éléments finis proviennent des maillages non réguliers comme nous allons le voir au paragraphe suivant.
### 1.4. Parallélisation par macro-éléments

Si le maillage est quelconque et qu'on met un processeur par noeud ou par élément on ne peut pas en général mettre en bijection le maillage et le réseau de processeur en conservent la propriété :

> "A tout noeud voisin correspond un processeur voisin"

Donc pour calculer $(∇φ_h, ∇ψ_h)$ par exemple, ou $φ_h$ et $ψ_h$ sont deux fonctions de $V_{oh}$ ayant leur valeurs $φ_i^j$ et $ψ_i^j$ aux noeuds $q_i^j$ stockés dans le ième processeur, il faut transférer dans le processeur i les valeurs $φ_j^k, ψ_j^k$ qui sont dans le processeur j, lequel peut très bien se trouver à l'autre bout du réseau physique. Le processus sera donc très cher : $3√P$ transferts par itération de temps en moyenne avec une programmation inextricable.

Il faut donc récupérer à un niveau plus élevé une certaine régularité dans le maillage.

**Une methode consiste à construire le maillage en 2 étapes :**

* Etape 1 : **Quadranguler le domaine**
* Etape 2 : Trianguler chaque trapèze.

Pour que ce procédé donne une triangulation admissible il faut que les triangulations des trapèzes se raccordent. Pour cela il suffit de choisir les noeuds qui seront sur les frontières des rectangles et de construire les triangulations à partir de ces noeuds (cf. figure 4). Cette méthode est connue sous le nom de méthode des super-éléments.

<figure><img src="image_8.png" /><figcaption>**Figure 4 :** Quadrangulation d'un domaine et subdivision d'un trapèze en triangles.</figcaption></figure>
Maintenant nous <u>associons à chaque trapèze un processeur.</u>

En procédant de cette manière on a conservé la propriété : trapèze voisin ⇏
processeur voisin ; et l'accès à une donnée correspondant à un trapèze voisin
se fait directement pour les trapèzes voisins nord, sud, est, ouest et à
2 niveaux pour les trapèzes voisins nord-est, sud-ouest, nord-ouest, sud-ouest.

Montrons comment organiser les calculs pour effectuer :

$$
(1.19) \quad \varphi_j^{n+1} = \Delta t f_j + \varphi_j^n - \Delta t (\nabla \varphi_h^n, \nabla_w^j)/\sigma_j \quad \forall j \in I
$$

où $(\nabla \varphi_h^n, \nabla w^j)$ est calculé par (1.17).

Posons :

`$\{R_k\}_{k=1}^K$ ensemble des trapèzes`

`$\{q_i^i\}_{i \in I_k}$ ensemble des noeuds <u>intérieurs</u> au trapèze k`

`$\{q_i^i\}_{i \in F_k}$ ensemble des noeuds frontières du trapèze k.`

Si $j \in I_k$ (1.19) se calcule entièrement dans le processeur k sans avoir besoin
de donnée extérieure à condition d'avoir stocké au préalable dans la mémoire
du processeur k: $f_j, \varphi_j^n, \Delta t$ et : $\varphi_i^n, (\nabla w^i, \nabla w^j), i \in I_k \cup F_k$ ;

Cout du calcul de (1.19) pour tous les $j \in I_k$ :
NSI \times (2 \times NV + 6)

où

NSI = nombre d'éléments de $I_k$

NV = nombre moyen de noeuds voisins des noeuds de $I_k$.

Si $j \in F_k$ alors il faut transférer les $\varphi_j^i$ pour $q_j^i$ voisin de $q_j^j$ ; ce qui
coute en moyenne environ 3xNV transferts. Si on a stocké dans le processeur
k les $(\nabla w^i, \nabla w^j)$ correspondants à $j \in F_k$ et $q_j^i$ voisin de $q_j^j$ alors le coût
total pour effectuer (1.19) pour tous les $j$ de $F_k$ est de l'ordre de
$NSF \times (2 \times NV + 6 + 3 \times NV)$

où

$NSF$ = nombre de noeuds de $F_k$.

En fait on associera à chaque processeur $k$ les noeuds de la moitié de la frontière de $R_k$.

Donc au total pour résoudre le problème il faudra :

$$CPU_t = \left[ NSF \times (2 \times NV + 6) + \left(\frac{NSF}{2}\right) \times (2 \times NV + 6 + 3 \times NV) \right] \times \frac{N^2 \times T}{\min(a,b)} \times \mu$$

**Exemple :**

Avec les mêmes paramètres que ci-dessus on prend

$$NSF = 4, NV = 6, NSF = 12$$

ce qui donne à peu près $10^6$ noeuds pour 10.000 processeurs. On a alors :

$$CPU_t = 48 \text{ mn}$$

## 2. Résolution du problème de Dirichlet par éléments finis

### 2.1. Présentation de la méthode

Considérons maintenant le problème suivant :

$$ (2.1) \quad \begin{cases} - \Delta \varphi = f \quad \text{dans } \Omega \\ \varphi|_{\Gamma_h} = 0 \end{cases} $$

c'est-à-dire, sous forme variationnelle, avec les notations (1.10), (1.11) :

$$ (2.2) \quad \begin{cases} (\nabla \varphi, \nabla w) = (f, w) & \forall w \in H_0^1(\Omega) \\ \varphi \in H_0^1(\Omega) \end{cases} $$

Construisons une triangulation comme au paragraphe 1.

Soit $R_h$ une quadrangulation de $\Omega_h$ avec les propriétés habituelles

$$ \cdot \Omega_h = \bigcup_{k=1}^{K} R_k, \quad R_k \in R_h. $$

$$ \cdot R_k \cap R_l = 1 \text{ côté, } 1 \text{ noeud ou le } \emptyset, \forall k,l = 1, \dots, K \quad k \neq l $$

• $h$ = longueur de plus grand côté ; distance entre $F_h$ et $\Gamma$ de l'ordre de $h^2$ ; tous les noeuds de $\Gamma_h$ sont sur $\Gamma$.
Soit maintenant $\mathcal{J}_{kh}$ des triangulations des $R_k$ avec les propriétés du § 1.3 et en plus :

• Si $R_k \cap R_\ell \neq \emptyset$ alors tout noeud de $\partial R_k \cap \partial R_\ell$ appartient à $\mathcal{J}_{kh}$ et à $\mathcal{J}_{\ell h}$.

On denote comme précédemment par $\{q^i\}_1^M$ les noeuds de $\mathcal{J}_h = \cup_k \mathcal{J}_{kh}$,
par $I_k$ l'ensemble des noeuds intérieurs à $R_k$, par $F_k$ les noeuds frontières
de $R_k$ et par

(2.3) $\mathcal{J}_i$ l'ensemble des indices des noeuds de $\mathcal{J}_h$ qui ont $q^i$ pour voisin.

Comme d'habitude (2.2) discrétisé par des éléments finis $P^1$ donne :

$$
(2.4) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
(\nabla \varphi_h, \nabla w_h) & = (f_h, w_h) \quad \forall w_h \in V_{oh} \\[1em]
\varphi_h \in V_{oh} & = \{w_h \in C^0(\Omega_h) : w_h|_T \text{ affine } \forall T \in \mathcal{J}_h ; w_h|_{\Gamma_h} = 0\}
\end{array}
\right.
$$

c'est-à-dire :

$$
(2.5) \quad A \Phi = F
$$

où

$$
(2.6) \quad \Phi = (\varphi_1, \dots, \varphi_N)^t ; \quad \varphi_i = \varphi_h (q^i)
$$

$$
(2.7) \quad N = \text{nombre de noeuds de } \mathcal{J}_h \text{ n'appartenant pas à } \Gamma_h.
$$

$$
(2.8) \quad a_{ij} = (\nabla w^i, \nabla w^j)
$$

$$
(2.9) \quad w^i \in V_{oh} ; \quad w^i(q^j) = \delta_{ij} \quad \forall i,j = 1, \dots, N.
$$

$$
(2.10) \quad F = (f_1, \dots, f_N)^t ; \quad f_i = (f, w^i).
$$

Pour résoudre (2.5) on va utiliser la technique des macro-éléments.
La technique consiste à éliminer les inconnues $\varphi_i$ pour i correspondant aux
noeuds internes des $R_k$.
Soit $i \in I_k$ l'indice d'un noeud interne de $R_k$. L'ensemble des équations dans (2.5) qui contiennent $\varphi_i$ sont :

$$ (2.11) \quad \sum_{\ell \in J_j} a_{j\ell} \varphi_\ell = f_j \quad j \in J. $$

On peut par exemple utiliser la i-ème équation pour calculer $\varphi_i$ :

$$ (2.12) \quad \varphi_i = a_{ii}^{-1} \left( f_i - \sum_{\substack{j \neq i \\ \ell \in J_i}} a_{ij} a_{ii}^{-1} \varphi_\ell \right) $$

et remplacer dans (2.11) pour obtenir :

$$ (2.13) \quad \sum_{\substack{l \in J_j \\ l \neq i}} a_{jl} \varphi_l - a_{ji} a_{ii}^{-1} \sum_{\substack{l \in J_i \\ l \neq i}} a_{il} \varphi_l = f_j - a_{ji} a_{ii}^{-1} f_i, \quad j \neq i, j \in J_i. $$

Plus généralement si (2.5) s'écrit :

$$ (2.14) \quad \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} \Phi_f \\ \Phi_i \end{pmatrix} = \begin{pmatrix} F_f \\ F_i \end{pmatrix} $$

où $\Phi_i$ est le vecteur des composantes correspondant aux noeuds internes.

On peut se ramener à :

$$ (2.15) \quad (A_{11} - A_{12} A_{22}^{-1} A_{21}) \Phi_f = F_f - A_{12} A_{22}^{-1} F_c $$

Puis une fois $\Phi_f$ connu retrouver $\Phi_c$ par :

$$ (2.16) \quad \Phi_c = A_{22}^{-1} (F_c - A_{21} \Phi_f) $$

Bien sûr dès que le nombre d'éléments de $I_k$ est grand il est préféré d'utiliser des techniques d'élimination par factorisation (cf. BOSSAVIT-MEYER [2] par exemple).
La résolution de (2.15) sur un réseau orthogonal de microprocesseur ne pose pas de problème de transfert de données puisqu'on est maintenant sur le réseau régulier {$R_k$} type différences finies. On renvoie à Hellen [3] pour une étude détaillée de ce problème.

### Remarque 2.1

Moyennant quelques transferts à 1 niveau on peut aussi éliminer les noeuds de $F_k$ qui ne sont pas des sommets de $R_k$. Il ne restera plus que les inconnues relatives aux sommets des trapèzes. Si on affecte le noeud sud-ouest par exemple de chaque $R_k$ au processeur k alors le système à résoudre est

$$ (2.17) \quad \sum_{\ell \in I_k} \tilde{a}_{k\ell} \varphi_{\ell} = \tilde{f}_k \quad k = 1, \dots, K $$

où

$$ (2.18) \quad L_k = \{ \ell : q^\ell \text{ et } q^k \in R_k, \ k' = 1, \dots, K \} $$

### Remarque 2.2

Une méthode frontale parallèle (qui est une méthode directe) ou plus généralement toute méthode d'élimination paire-impaire serait particulièrement séduisante si en plus des liaisons nord-sud-est-ouest on avait des liaisons reliant un processeur sur $\sigma$, à plusieurs niveaux comme l'indique la figure 5.

<figure><img src="image_9.png" /><figcaption>**Figure 5** : Liaisons supplémentaires pour une résolution par élimination de Gauss.</figcaption></figure>

Les liaisons sont représentées dans l'espace pour plus de clarté.

* = noeuds éliminés au premier tour

x = " " " second tour.
### Remarque 2.3

Il serait aussi intéressant d'étudier le parallélisme possible des méthodes multi-grilles. Celles-ci nécessitent aussi les liaisons représentées sur la figure 5.

## 2.2. Evaluation et programmation de la méthode.

La construction de la matrice du système (2.15) se fait en parallèle sur chaque $R_k$. Ayant associé un processeur pour chaque $R_k$ on envoie dans le processeur les informations suivantes :

• coordonnées des noeuds $q^i$ de $J_{kh}$ et des $q^j$ voisins des $q^k \in \partial R_k$ :
$Q(L,I), L = 1,2 ; I=1,...,NSV$

• valeurs de $f$ aux noeuds des formules de quadrature pour tous les triangles de $J_{kh}$ et ceux qui ont un côté ou un sommet appartenant à $\partial R_k$ (cf. figure 6).

Par exemple si les noeuds de quadrature sont les sommets des triangles on enverra :
$F(I), I=1,...,NSV$

• les numéros des noeuds de chaque triangle :
$ME(L,J), L=1,...,3 ; J=1,...,NTV$.

<figure><img src="image_10.png" /><figcaption>Figure 6 : $J_{kh}$ et les triangles voisins.</figcaption></figure>

• Un tableau d'état des noeuds valant 0 si le noeud est interne, 1 si le noeud est sur $\partial R_k$, 2 si le noeud est voisin de $\partial R_k$ à l'extérieur de $R_k$, NG(I), $I=1,...,NSV$.

Pour construire A et F on utilise les formules (1.15), (1.17) en bouclant sur les triangles ce qui demande environ 46 opérations par triangle, soit sur chaque $R_k$ :

$46 \times NTV$ cycles de calcul.
Pour éliminer les équations et inconnues correspondant aux noeuds internes, il faut à peu près

$(NSF/2)^2/2$

$$3 \text{ NSI} \sum_i \approx \text{NSI} \times \text{NSF}^2 \times \max\{5, \frac{\text{NSF}^2}{48}=5\} \text{ cycles de calcul.}$$
$$\text{NSR} \times \text{NV}/2$$

où NSI est le nombre de noeuds internes de $R_k$, NSR le nombre de noeuds de $R_k$, NSF le nombre de noédus frontière de $R_k$ et NV le nombre moyen de noeuds voisin de chaque noeud. Cette estimation ne serait pas valable pour des grands NSR car il faudrait utiliser l'algorithme de factorisation LV (cf. BOSSAVIT-MEYER [2] par exemple).

Pour effectuer cette élimination et pour pouvoir recalculer $\varphi_i$ aux
noeuds internes, il faut stocker :

• la matrice $A : A(I,J), I=1,...,NSR, J=1,...,NV$

• un tableau de voisinage de chaque noeud : $M(I,J), I=1,...,NSR, J=1,...,NV$ (ce n'est pas la seule méthode).

• la matrice de (2.15) $A_1(I,J), I,J=1,...,NSF$, le second membre
$F(I) \ I=1,...,NSF$, et une matrice de calcul des $\varphi_i$ internes $L(I,J) \ I=1,...,NSF,$
$J=1,...,NSI$.

Même remarque que ci-dessus : cette méthode est mauvaise si NSR est plus grand que quelques dizaines. Ces évaluations peuvent facilement se vérifier expérimentalement puisque les programmes d'élimination de noeud interne et de construction des matrices éléments finis existent.

Plaçons nous dans les conditions optimales et supposons qu'on a éliminé aussi les noeuds frontières de $R_k$ sauf les 4 coins. Il reste alors P inconnues. Pour résoudre le système (2.15), il faut de l'ordre de P/10 itérations pour une méthode itérative du type relaxation et chaque itération coute au moins une évaluation du système linéaire c'est-à-dire 18 opérations et 12 transferts élémentaires par processeur ; on obtient environ

$3P$ cycles de calcul

il reste le calcul des $\varphi_i$ internes soit : $2 \times \text{NSF} \times \text{NSI}$ opérations.

Donc au total pour toute la résolution, il faudra :
$$CPU_t = [3P + NSI \times NSF^2] \times \max\{5, \frac{NSF^2}{48} - 5\} + 46 \times NTV + 2 \times NSF \times NSI] \mu$$

et à peu près

$$2 \text{ NSR} \times \text{NV} + \text{NSF}^2 + (\text{NSI}+1)\text{NSF} + 3 \times \text{NTV} + 4 \times \text{NSV}$$

mots de mémoire par processeur.

### Exemple :

Avec

$$P = 10.000, \mu = 10^{-5}, \text{NSR} = 16, \text{NSF} = 12, \text{NV} = 6, \text{NSV} = 28, \text{NTV} = 50$$

on obtient

$$CPU_t = [310^4 + 2500 + 2300 + 96]10^{-5} = .35''$$

On voit que c'est la résolution du système linéaire (2.15) qui coute le plus cher. Si on a pas éliminé les noeuds frontières de $R_k$ alors il faut multiplier le temps calcul par $(NSF)/2$ ; on obtient $CPU_t = 2.1''$.

La capacité mémoire nécessaire par processeur est de l'ordre de 660 mots.

### Remarque 2.4

Si l'élimination de Gauss pair-impair était possible il faudrait à peu près 40 $\log_2 P$ cycles de calcul : avec les chiffres ci-dessus on obtient

$$CPU_t = [1000 + 2500 + 2300 + 96]10^{-5} \approx .06''$$

### Remarque 2.5

Pour que cette méthode soit facile à programmer il faut que, à part quelques exceptions éventuellement, tous les processeurs effectuent le même programme. Il faut donc que tous les $R_k$ aient le même nombre de triangles et le même nombre de noeuds. Il faudrait regarder ce point plus en détails mais à première vue il semble possible de rajouter des triangles et des noeuds fantomes (triangle d'aire nulle ou très petite et noeud d'indice NG = 3 qui ne comptent pas dans le calcul mais permettrait d'avoir les mêmes indices de boucle pour tous les processeurs).
### Remarque 2.6

Pour les calculs tridimensionnels les estimations sont à peu près identiques mais NV est de l'ordre de 24 et la disproportion entre NSF et NSF - NSR est encore plus grande.

## 3. Résolution des équations de Navier-Stokes

La maitrise des systèmes météorologiques, aérodynamiques et d'échanege thermique par circulation de fluide passe par la maitrise numérique des équations de Navier-Stokes :

$$ (3.1) \quad \left\{ \begin{array}{l} \frac{\partial u}{\partial t} - v \Delta u + u \nabla u + \nabla p = f \\ \nabla \cdot u = 0 \\ u(t=0) = u_0 \quad \text{dans } \Omega, u = u_0 \quad \text{sur } \Gamma \times ]0,T[. \end{array} \right. \quad \text{dans } \Omega \times ]0,T[ $$

Pour la plupart des écoulements réels (nombre de Reynolds très grand) et surtout si l'écoulement est turbulent, ces équations sont très difficiles et onéreuses à résoudre numériquement. On se propose donc d'étudier brièvement quel gain on peut espérer d'un calculateur parallèle pour ce genre de problème.

### 3.1. Rappel des méthodes classiques et comparaison

Il n'existe pas de bonne méthode explicite pour résoudre (3.1) (ce n'est pas vrai des écoulements compressibles) puisque la deuxième équation n'a pas de dérivées en t.

La plupart des méthodes semi-explicites à pas fractionnaire applicable aux écoulements tridimensionnels partent de la remarque suivante :

$$ (3.2) \quad \Delta p = \nabla \cdot (\nabla p) = \nabla \cdot f - \nabla \cdot (u \nabla u) \quad \text{dans } \Omega $$

En rajoutant une condition aux limites à (3.2), ce qui malheureusement n'est pas facile, on peut connaisant $u$, trouver $p$.

Un schéma simple sera par exemple (si $u^0$ est indépendant de $t$) :
$$
(3.3) \quad \left\{
\begin{array}{l}
\displaystyle \frac{1}{\Delta t} (u^{m+1} - u^m) - v \Delta u^{m+1} + u^m \nabla u^m + \nabla p^m = f^m , \\[1em]
u^{m+1} = u_0 \text{ sur } \Gamma , \\[1em]
\Delta p^{m+1} = \nabla \cdot f^{m+1} - \nabla \cdot (u^{m+1} \nabla u^{m+1}) \\[1em]
\displaystyle \frac{\partial p^{m+1}}{\partial n} = f^{m+1} n + v (\Delta u^{m+1}) n - u_0 \nabla u^{m+1} \text{ n on } \Gamma
\end{array}
\right.
$$

Une autre méthode se ramène à chaque itération à la résolution d'un
problème de Stokes : on ne découple pas $u$ et $p$ ;

$$
(3.4) \quad \left\{
\begin{array}{l}
\displaystyle \frac{1}{\Delta t} (u^{m+1} - u^m) - v \Delta u^{m+1} + \nabla p^{m+1} = f^m - u^m \nabla u^m \\[1em]
\qquad \nabla \cdot u^{m+1} = 0 \\[1em]
\qquad u^{m+1} = u^0 \quad \text{sur } \Gamma .
\end{array}
\right.
$$

Enfin la méthode suivante, quasi-implicite, se ramène à la résolu-
tion d'un problème d'Oseen :

$$
(3.5) \quad \left\{
\begin{array}{l}
\displaystyle \frac{1}{\Delta t} (u^{m+1} - u^m) - v \Delta u^{m+1} + \frac{1}{2} (u^{m+1} \nabla u^m + u^m \nabla u^{m+1}) + \nabla p^{m+1} = f^m \\[1em]
\qquad \nabla \cdot u^{m+1} = 0 \\[1em]
\qquad u^{m+1} = u^0 \quad \text{sur } \Gamma .
\end{array}
\right.
$$

Chacune de ces trois méthodes peut être modifiée par l'utilisation
d'un schéma du type Crank-Nicholson ou même "saute-mouton" au lieu du schéma
classique pour $\frac{\partial u}{\partial t}$.

Enfin il existe une méthode implicite par (GLOWINSKI-PERIAUX-PIRONNEAU [4])
par moindres carrés abstraits dont nous ne parlerons pas parce qu'elle est
trop longue à présenter et probablement difficile à programmer sur un ordi-
nateur parallèle.

La difficulté du schéma (3.3) vient de l'approximation de $\Delta u^{m+1}$ sur $\Gamma$.
En fait le schéma (3.3) peut être considéré comme un cas particulier du
schéma (3.4) où le système linéaire serait résolu par une itération d'une
méthode itérative.
Le schéma (3.5) est nettement plus stable que le schéma (3.4) mais il demande de résoudre à chaque itération au système linéaire non symétrique.

On a vue que c'était de loin les résolutions de système linéaire qui étaient les plus onéreuses sur un ordinateur étudié _sauf si celui-ci est équipé d'un réseau de communication_ de la figure 5. Même dans ce dernier cas la résolution d'un système non symétrique n'est pas plus cher que celle d'un système symétrique (sauf peut être pour les problèmes de pivotage, éventuellement). Du point de vue étudié au paragraphe précédent, le schéma (3.5) semble donc le meilleur.

## 3.2. Etude du schéma (3.5)

Etant donné que la construction des matrices des systèmes linéaires n'est pas un facteur déterminant pour le temps calcul on va discrétiser (3.5) par éléments finis car en mécanique des fluides les frontières sont souvent compliquées.

Le système (3.5) est donc discrétisé par : chercher $u_h \in V_h, P_h \in Q_h$ tels que :

$$ (3.6) \quad \left\{ \begin{array}{l} \displaystyle \frac{1}{\Delta t} (u_h^{m+1}, v_h) + v (\nabla u_h^{m+1}, \nabla v_h) + \frac{1}{2} (u^{m+1} \nabla u^m + u \nabla u^{m+1}) (v_h) + (\nabla p_h^{m+1}, v_h) = (f^m, v_h) \\[1em] \qquad + \frac{1}{\Delta t} (u_h^m, v_h) \quad \forall v_h \in V_{oh} \\[1em] \quad (\nabla \cdot u_h^{m+1}, q_h) = 0 \quad \forall q_h \in Q_h \end{array} \right. $$

avec $u_h^m(q^i) = u^0(q^i) \quad \forall m, \forall q^i \in J_h \cap \Gamma_h$ ; $u_h^0$ donné.

où

$$ V_{oh} = \{ v_h \in V_h : v_h = 0 \text{ sur } \Gamma_h \} $$

et où $Q_h$ et $V_h$ sont des espaces de dimensions finies destinés à approximer $L^2(\Omega)$ et $(H^1(\Omega))^2(\text{car } \Omega \subset \mathbb{R}^2)$.

Il y a plusieurs façons de discrétiser $L^2(\Omega)$ et $(H^1(\Omega))^2$ pour obtenir des schémas convergents pour (3.5). Nous n'examinons que les choix les plus simples avec des triangles.
### Choix 1 :

$$ 
\begin{aligned}
Q_h &= \{ q_h : q_h = \text{constant} \text{ sur } T \quad \forall T \in \mathcal{J}_h \} \\
V_h &= \{ v_h \in C^0(\bar{\Omega}_h)^2 : v_h|_T \in P^2 + \text{"bulle"} \quad \forall T \in \mathcal{J}_h \}
\end{aligned}
$$

où $P^k$ est l'ensemble des polynômes de degré $k$ en $(x_1, x_2)$ et la "bulle" (cf. RAVIART [5]) est le polynôme de degré 3 représenté en coordonnées barycentriques par $\lambda_1, \lambda_2, \lambda_3$.

### Choix 2 :

$$ 
\begin{aligned}
Q_h &= \{ q_h \in C^0(\bar{\Omega}_h) : q_h|_T \in P^1 \quad \forall T \in \mathcal{J}_h \} \\
V_h &= \{ v_h \in C^0(\bar{\Omega}_h)^2 : v_h|_T \in P^2 \quad \forall T \in \mathcal{J}_h \}
\end{aligned}
$$

### Choix 3 :

$$ 
\begin{aligned}
Q_h &= \{ q_h \in C^0(\bar{\Omega}_h) : q_h|_T \in P^1 \quad \forall T \in \mathcal{J}_h \} \\
V_h &= \{ v_h \in C^0(\bar{\Omega}_h/2)^2 : v_h|_T \in P^1 \quad \forall T \in \mathcal{J}_h/2 \}
\end{aligned}
$$

où $\mathcal{J}_{h/2}$ est la triangulation obtenue en divisant tout triangle de $\mathcal{J}_h$ en 4 par les milieux des côtés.

Le choix le plus simple à programmer est le choix 2 mais celui qui donne le minimum d'opérations arithmétiques est le choix 3. Ces choix ne sont pas déterminants pour le temps calcul, nous allons donc évaluer la méthode avec le choix 3 car les opérations sont plus faciles à compter.

On vérifiera que (3.6) devient :

$$ (3.7) \quad \left[ \begin{array}{c|c} A + B & C \\ \hline t_C & 0 \end{array} \right] \quad \left[ \begin{array}{c} v \\ \hline P \end{array} \right] = \left[ \begin{array}{c} F \\ \hline 0 \end{array} \right] $$

où

$$ v = (v_1, v_2, \dots, v_N)^t, \quad P = (p_1, \dots, p_M) $$

$$ u_h = \sum (u_i \cdot v_i)^t w_i^i \quad ; \quad p_h = \sum p_i w_i^i \dots (w_i^i \text{ fonction de base de } \mathcal{J}_{h/2}) $$
et où

$$
a'_{ij} = \frac{1}{\Delta t} (w_{\frac{1}{2}}^i, w_{\frac{1}{2}}^j) + v (\nabla w_{\frac{1}{2}}^i, \nabla w_{\frac{1}{2}}^j) ; \quad A = \begin{pmatrix} A' & 0 \\ 0 & A' \end{pmatrix}
$$

$b_{ij}$ est une matrice antisymétrique dont les termes sont de la forme

$$
w_{\frac{1}{2}}^{i} \frac{\partial u^{m}}{\partial x_{k}} + u^{m} \nabla w_{\frac{1}{2}}^{i} .
$$

$$
\vec{c}_{ij} = -(w_{\frac{1}{2}}^{i}, \nabla w^{j}) .
$$

Les calculs pour calculer la matrice de (3.7) sont du même style
que ceux du calcul d'un laplacien vectoriel en 3 dimensions avec couplage.

Par quadrilatère $R_k$, il faudra donc à peu près :

2 × NV2 × NSR2 + 46 × NTV2 opérations pour construire A et F
ou NV2, NSR2, NTV2 se rapportent à $\frac{\pi}{4h}/2$.

37 × NTV2 opérations pour construire B

9 × NTV2 opérations pour construire C

Pour éliminer les inconnues relatives aux noeuds internes il faut
à peu près

$$
(NSI + 2 \times NSI2) \times (3NSF)^2 \times \max\left\{5, \frac{(3NSF)^2}{48} - 5\right\} \text{ cycles de calcul.}
$$

Puis pour résoudre le système linéaire, après élimination des noeuds frontières
qui ne sont pas des sommets de $R_k$ :

14P cycles de calcul.

Soit au total :

$$
\text{CPU}_t = (14\text{P} + 9 \times (\text{NSI} + 2 \times \text{NSI2}) \times \text{NSF}^2 \times \max\{5, \frac{3}{16} \text{NSF}^2 - 5\} + 92 \times \text{NTV2} + 2 \times \text{NV2} \times \text{NSR2})\mu \\
\text{par itération de temps.}
$$

Exemple

On reprend les chiffres précédents :

$$
\begin{align*}
P &= 10.000, \mu = 10^{-5}, \text{NSR2} = 16, \text{NSF2} = 12, \text{NV2} = 6, \text{NSV2} = 28, \\
\text{NSR} &= 8, \text{NSF} = 6, \text{NV} = 6, \text{NSV} = 14, \text{NTV2} = 50
\end{align*}
$$
$$CPU_t = 1.4 + 16 + 046 + 0034 \approx 1.6''$$

On pourra donc faire en un temps raisonnable quelques centaines d'itérations de temps avec un tel système :

### Remarque 3.1

Même avec ce système complexe c'est toujours la résolution du système linéaire qui prend le plus de temps. Si l'élimination de Gauss pair-impair complète était possible alors $CPU_t \approx .3$".

C'est maintenant la construction des matrices qui est le plus onéreux ; on pourra alors peut être utiliser le choix 1 car on peut pratiquement calculer $p_h$ par

$$p_h \approx - \frac{1}{\epsilon} \nabla \cdot u_h \qquad \epsilon < 1$$

### Remarque 3.2

Le décentrage (FORTIN [6]) semble particulièrement adapté à ce type de système car il augmente la stabilité de la méthode sans pénaliser sensiblement le temps calcul.

### Remarque 3.3

Pour modéliser les écoulements à grand nombre de Reynolds on utilise des **modèles de turbulence**. Une méthode simple par exemple consiste à remplacer $- \nabla \cdot \Delta u$ par $- \nabla \cdot ((v+v_t) |\nabla \times u|^2) \nabla u$.

Si on traite ce terme par

$$- \nabla \cdot ((v+v_t) |\nabla \times u|^2) \nabla u^{m+1})$$

alors la méthode précédente est encore applicable sans grand changement.
<table><thead><tr><th></th><th>Equation de la chaleur<br>schéma explicite</th><th>Problème de Dirichlet</th><th>Equation de Navier-Stokes<br>(par itération de temps implicite)</th></tr></thead><tbody><tr><td>Eléments finis</td><td>48'</td><td>0.35" 0.06"</td><td>1.6" 0.3"</td></tr><tr><td>Même schéma en<br>différences<br>finies</td><td>24'</td><td>0.30" 0.03"</td><td>?</td></tr><tr><td></td><td>1 2</td><td>1 2</td><td>1 2</td></tr></tbody></table>

Tableau récapitulatif des estimations de temps calcul pour des équations aux dérivées partielles sur une grille avec $10^5$ noeuds.

La colonne 2 referre au cas où on a les liaisons décrites dans la figure 5.
4. Conclusion

Après avoir évalué assez grossièrement le nombre d'opérations nécessaires pour résoudre 3 équations aux dérivées partielles classiques nous avons dégagé les points suivants :

Les méthodes d'éléments finis se parallélisent aussi bien que les méthodes aux différences finies avec un léger avantage pour ces der nières s'il n'y a pas de systèmes linéaires à résoudre (schémas explicites).

La partie la plus onéreuse d'un tel calcul est la résolution des systèmes linéaires sauf si le réseau de processeurs possède des connections permettant l'élimination de Gauss pair-impair.

Il reste bien sur à vérifier ces chiffres, ce qui est possible en faisant des tests sur le D.A.P. d'ICL. Enfin il faut tout de même souligner que les méthodes d'éléments finis sont plus difficiles à programmer et demandent plus de mémoire.
References :

[1] FRAILONG J.M., Definition of an integrated processing element for a parallem numerical architecture (to appear).

[2] BOSSAVIT A., MEYER B., Construction d'un programme de factorisation sélective. Proc. Conf. GAMNI 3-1979.

[3] HELLER D., A survey of parallel algorithms in numerical linear algebra, SIAM Review, Oct. 78, Vol. 20 no 4.

[4] GLOWINSKI R., PERIAUX J., PIRONNEAU O., Abstract least squares. Proc. Conf. "finite element for mixte problems", Tel Aviv 76. N. GEFFEN Ed.

[5] RAVIART P.A., Lectures Notes for the Navier-Stokes equations, Univ. of Paris 6, LAN Report 189.

[6] FORTIN M., Résolution numérique des équations de Navier-Stokes par éléments finis de type mixte, Rapport Laboria N° 184.