Méthodes multigrilles

1. Introduction

La discrétisation d'une équation équation aux dérivées partielles linéaire conduit à la résolution d'une équation linéaire de la forme :

A U = F (1)

où U est une matrice colonne contenant les inconnues, c'est-à-dire les valeurs de la fonction inconnue aux points du maillage. Voir par exemple le cas de l'équation de Poisson bidimensionnelle en coordonnées cartésiennes.

La méthode de Gauss-Seidel permet d'obtenir la solution par itérations. Alors qu'un maillage très fin est nécessaire pour obtenir une solution détaillée, celui-ci conduit à une convergence très lente des basses fréquences spatiales. Lorsque la taille de la maille est augmentée, la vitesse de convergence des basses fréquences est considérablement augmentée. Les méthodes multigrilles ([1]) consistent à utiliser successivement des grilles (ou maillages) de différentes tailles, de manière à obtenir une solution détaillée dans les hautes fréquences, tout en assurant une relaxation rapide des basses fréquences.

On s'intéresse plus particulièrement aux équations à deux dimensions. Les inconnues sont donc représentées sur un maillage à deux dimensions. On note U_i,j l'inconnue au nœud (i,j) du maillage.

Les dimensions de la maille sont Δx et Δy. Pour alléger les notations, on considère des mailles carrées et on pose :

h = Δ x = Δ y (2)

2. Méthode à deux grilles

2.a. Principe

Soit h le pas de la grille sur laquelle on cherche la solution. Le système à résoudre s'écrit :

A^{h} U^{h} = F^{h} (3)

Dans cette notation, h est bien sûr un indicateur de la grille et non un exposant. La méthode à deux grilles consiste à utiliser aussi la grille de pas 2h pour réduire les erreurs basses fréquences lors du processus d'itération. Soit ${\tilde{U}}^{h}$ la solution approchée obtenue après des itérations sur la grille de pas h. L'erreur est par définition :

E^{h} = {\tilde{U}}^{h} - U^{h} (4)

où U^h désigne la solution exacte sur cette grille.

Si on définit le résidu par

R^{h} = F^{h} - A^{h} {\tilde{U}}^{h} (5)

alors l'erreur vérifie l'équation :

A^{h} E^{h} = R^{h} (6)

L'idée de la méthode à deux grilles est d'utiliser la grille 2h pour évaluer l'erreur $E^{h}$ . Il s'agit donc de résoudre l'équation

A^{2 h} E^{2 h} = R^{2 h} (7)

Pour cela il faut utiliser la matrice A^2h correspondant à la discrétisation de l'équation différentielle sur la grille 2h. Il faut aussi disposer du résidu sur cette grille, ce qui nécessite une méthode permettant d'obtenir le résidu sur la grille 2h à partir du résidu de la grille h. Cette opération est appelée la restriction et peut être notée de la manière suivante :

R^{2 h} = I_{h}^{2 h} R^{h} (8)

Lorsque la restriction est faite, on effectue des itérations de Gauss-Seidel sur la grille 2h afin d'obtenir une estimation ${\tilde{E}}^{2 h}$ de l'erreur sur cette grille. L'avantage d'utiliser pour cela la grille 2h est d'obtenir une relaxation plus rapide des erreurs de basses fréquences spatiales.

L'erreur obtenue sur la grille 2h doit être extrapolée à la grille h. Cette opération est appelée la prolongation :

{\tilde{E}}^{h} = I_{2 h}^{h} {\tilde{E}}^{2 h} (9)

Il reste à corriger la solution approchée avec cette erreur :

{\tilde{U}}^{h} \leftarrow {\tilde{U}}^{h} + {\tilde{E}}^{h} (10)

puis à continuer les itérations sur la grille h. Le passage à la grille 2h sera répété régulièrement jusqu'au niveau de convergence souhaité.

2.b. Échanges entre les grilles

Les grilles ont dans les deux directions un nombre de mailles égal à une puissance de 2. Pour simplifier les notations, on suppose que les deux dimensions sont égales (grille carrée). Le nombre de points sur chaque dimension pour la grille de pas h est N^h=2^p+1. Sur la grille de pas 2h il est N^2h=2^p-1+1. La correspondance entre les deux grilles apparait sur la figure suivante (p=2) :

Figure pleine page

Les bords et les coins du domaine sont définis par des nœuds des deux grilles. Il en est de même des objets définis à l'intérieur du domaine.

Le point (i,j) de la grille 2h coïncide avec le point (2i,2j) de la grille h. À chaque point (i,j) de la grille 2h, on peut associer 4 points de la grille h, par exemple les points (2i,2j), (2i+1,2j), (2i,2j+1) et (2i+1,2j+1). Pour les points situés sur les bords à droite et en haut, il y a seulement deux points de maille h pour chaque point de maille 2h.

L'opération de prolongation, c'est-à-dire le passage de l'erreur de la grille 2h à la grille h peut être réalisée par une interpolation linéaire :

\begin{matrix} {\tilde{E}}_{2 i, 2 j}^{h} = {\tilde{E}}_{i, j}^{2 h} & (11) \\ {\tilde{E}}_{2 i + 1, 2 j}^{h} = \frac{1}{2} ({\tilde{E}}_{i, j}^{2 h} + {\tilde{E}}_{i + 1, j}^{2 h}) & (12) \\ {\tilde{E}}_{2 i, 2 j + 1}^{h} = \frac{1}{2} ({\tilde{E}}_{i, j}^{2 h} + {\tilde{E}}_{i, j + 1}^{2 h}) & (13) \\ {\tilde{E}}_{2 i + 1, 2 j + 1}^{h} = \frac{1}{4} ({\tilde{E}}_{i, j}^{2 h} + {\tilde{E}}_{i + 1, j}^{2 h} + {\tilde{E}}_{i, j + 1}^{2 h} + {\tilde{E}}_{i + 1, j + 1}^{2 h}) & (14) \end{matrix}

Pour appliquer ce schéma, il suffit de parcourir la grille 2h, en faisant varier les indices i et j de 0 à N^2h-1. Sur une architecture parallèle, chaque paire (i,j) peut être traitée par une unité de calcul.

L'opération de restriction, c'est-à-dire le transfert du résidu de la grille h à la grille 2h peut être simplement réalisée en recopiant les résidus de chaque nœud (2i,2j) de la grille h sur le nœud (i,j) correspondant de la grille 2h. Cette opération de restriction est appelée injection. Il est toutefois préférable d'utiliser une moyenne pondérée des points voisins ([1]) :

\begin{matrix} R_{i, j}^{2 h} = & \frac{1}{16} (R_{2 i - 1, 2 j - 1}^{h} + R_{2 i - 1, 2 j + 1}^{h} + R_{2 i + 1, 2 j - 1}^{h} + R_{2 i + 1, 2 j + 1}^{h}) & (15) \\ + \frac{1}{8} (R_{2 i, 2 j - 1}^{h} + R_{2 i, 2 j + 1}^{h} + R_{2 i - 1, 2 j}^{h} + R_{2 i + 1, 2 j}^{h}) & (16) \\ + \frac{1}{4} R_{2 i, 2 j}^{h} & (17) \end{matrix}

2.c. Cas des grilles périodiques

Il peut arriver qu'une grille présente une périodicité dans une direction. C'est le cas avec les coordonnées polaires. Dans ce cas, l'itération de Gauss-Seidel se fait avec une condition limite périodique sur les deux bords correspondants. Dans la direction de la périodicité, la correspondance entre les grilles h et 2h est modifiée, comme le montre la figure suivante :

Figure pleine page

Le nombre de points est multiplié par deux lorsqu'on passe de 2h à h. On peut donc prendre par convention un nombre de points de la forme 2^p. En comparaison, le nombre de points dans une direction non périodique est nécessairement de la forme 2^p+1. Cela revient à enlever la dernière rangée dans la grille rectangulaire.

3. Schémas multigrilles

3.a. Définition des grilles

Dans les schémas multigrilles, il faut disposer d'une cascade de grilles avec des pas de plus en plus grand. Pour y parvenir, on définit le problème à discrétiser sur la grille la plus grossière. Soit 2^q le nombre de mailles dans chaque dimension de cette grille. L'entier q est choisi le plus petit possible compte tenu de la complexité du problème. Le nombre de nœuds est 2^q+1.

La grille deux fois plus fine comporte 2^q+1 mailles dans chaque dimension, soit 2^q+1+1 nœuds. On continue la division de grille jusqu'à la grille la plus fine, qui comporte 2^p mailles. On dispose ainsi de p-q grilles de pas h, 2h, 4h jusqu'à 2^p-qh.

Au moment de la discrétisation, les objets sont définis sur la grille d'exposant q. La matrice A^h et le second membre F^h sont définis en cascade pour toutes les grilles jusqu'à la plus fine.

3.b. Schéma multigrilles en V

On commence par effectuer quelques itérations (1 ou 2) de Gauss-Seidel sur la grille la plus fine, de pas h, avant de transférer le calcul de l'erreur sur la grille 2h, comme expliqué plus haut.

Pour calculer l'erreur sur la grille 2h, on effectue quelques itérations sur cette grille, avant de transférer sur la grille 4h le calcul de l'erreur (erreur de l'erreur de la première grille). On procède ainsi jusqu'à la grille la plus grossière, de pas 2^p-qh. Sur cette grille, on procède à suffisamment d'itérations pour faire converger l'erreur (on peut même faire une résolution directe si q est petit). L'erreur obtenue sur cette grille est transférée sur la grille de pas 2^p-q-1h, sur laquelle on effectue quelques itérations avant de transférer l'erreur sur la grille de pas 2^p-q-2h. On procède ainsi jusqu'à la grille de pas h, où l'erreur vient corriger U^h.

On a ainsi décrit le cycle de parcours de grilles suivant :

h \to 2 h \to 4 h \dots \to 2^{p - q} h \to \dots 4 h \to 2 h \to h (18)

Figure pleine page

Dans une phase descendante du cycle, on effectue quelques itérations de Gauss-Seidel sur la grille 2h avant de calculer le résidu, que l'on transfère par restriction sur le second membre de l'équation de la grille 4h. Dans la phase ascendante, l'erreur obtenue sur la grille 4h est prolongée vers la grille 2h pour venir corriger l'erreur qui était stockée sur la grille 2h. Quelques itérations sont faites avant le passage à la grille plus fine. Sur la grille la plus grosse, on effectue des itérations jusqu'à convergence. On peut aussi faire une résolution directe sur cette grille si elle est assez petite.

Une fois le cycle en V terminé, on le répète si le niveau de convergence souhaité n'est pas atteint. À chaque étape sur une grille, le nombre d'itérations de Gauss-Seidel effectué est très réduit (1 ou 2). Cela est suffisant pour relaxer les modes dont la période spatiale est de l'ordre du pas de la grille. Pour les modes de période plus grande, il est bien plus efficace de passer à la grille deux fois moins fine.

Les paramètres à choisir pour un cycle en V sont : le nombre d'itérations pour obtenir la convergence sur la grille la plus grosse, le nombre d'itérations n₁ pour chaque étape dans le sens descendant, et le nombre d'itérations n₂ pour chaque étape dans le sens montant.

3.c. Schéma multigrilles complet

Dans le schéma multigrilles complet, on se sert des différentes grilles pour obtenir la solution du problème sur chacune. On commence par effectuer des itérations de Gauss-Seidel sur la plus grosse grille (disons 8h), pour le système AU=F, jusqu'à convergence vers la solution pour cette grille. Cette solution est prolongée vers la grille 4h, ce qui donne une première approximation de la solution pour cette grille. Cette approximation est corrigée au moyen de quelques cycles en V (par restriction du résidu). La solution corrigée est ensuite transférée à la grille 2h, et ainsi de suite jusqu'à l'obtention de la solution corrigée sur la grille la plus fine h.

Figure pleine page

Les paramètres à choisir pour ce schéma sont : le nombre d'itérations pour obtenir la convergence sur la grille la plus grosse, le nombre de cycles en V pour chaque correction, le nombre d'itérations n₁ pour chaque étape dans le sens descendant, et le nombre d'itérations n₂ pour chaque étape dans le sens montant.

4. Exemple

from pylab import *
import numpy
import poisson.main

On considère l'équation de Laplace bidimensionnelle en coordonnées cartésiennes, avec une condition limite de dirichlet (valeur nulle) sur les bords du domaine. Un segment avec une condition limite de Dirichlet est placé dans le domaine.

La grille la plus grosse sert à définir le problème. Elle est carrée, avec q=4, ce qui fait 17x17 points. Le nombre de grilles est 4, ce qui signifie que la grille la plus fine comporte 257x257 points.

laplace = poisson.main.Poisson(4,4,5,1.0,1.0)
laplace.laplacien()
laplace.dirichlet_borders(0.0)
laplace.dirichlet_polygon(6,3,[4],[0],1.0)

On effectue un cycle en V avec 2 itérations sur chaque grille avant le passage à une grille voisine dans le sens descendant, et 1 itération dans le sens montant. Sur la grille la plus grosse, on effectue 50 itérations. Le cycle est répété 10 fois.

result1 = laplace.multigrid_Vcycles_norm(1,4,50,2,1,10)
V1 = laplace.get_array()

Pour comparaison, on effectue aussi un schéma multigrille complet :

laplace.init_array()
result2 = laplace.multigrid_full_norm(3,5,50,2,1)
V2 = laplace.get_array()

Pour comparaison, on effectue aussi le même calcul avec des itérations de Gauss-Seidel avec sur-relaxation optimale seulement sur la grille la plus fine :

laplace.init_array()
result3 = laplace.opencl_iterations_norm(5,50,omega=laplace.get_optimal_omega())
V3 = laplace.get_array()

On trace la norme en fonction du nombre d'itérations. Pour les méthodes multigrilles, le nombre d'itérations est le nombre équivalent aux itérations sur la grille la plus fine. Par exemple, une itération sur la grille 2h équivaut à un quart d'itérations sur la grille h (car il y a quatre fois moins de points).

figure(figsize=(10,5))
plot(result1[0],result1[1],label='V cycle')
plot(result2[0],result2[1],label='Full')
plot(result3[0],result3[1],label='SOR opt')
legend(loc="lower right")
xlabel('niter')
ylabel('norm')
grid()

plotA.pdf

figure(figsize=(6,6))
contour(V1,20,extent=laplace.extent)
imshow(~laplace.get_mask(),extent=laplace.get_extent(),alpha=0.5, cmap=cm.gray)

plotB.pdf

figure(figsize=(6,6))
contour(V3,20,extent=laplace.extent)
imshow(~laplace.get_mask(),extent=laplace.get_extent(),alpha=0.5, cmap=cm.gray)
laplace.close()

plotC.pdf

Faisons la même comparaison avec une grille plus fine (1025x1025 points). On utilise alors 7 grilles.

laplace = poisson.main.Poisson(4,4,7,1.0,1.0)
laplace.laplacien()
laplace.dirichlet_borders(0.0)
laplace.dirichlet_polygon(6,3,[4],[0],1.0)
result1 = laplace.multigrid_Vcycles_norm(1,7,50,2,1,10)
V1 = laplace.get_array()
laplace.init_array()
result2 = laplace.multigrid_full_norm(3,7,50,2,1)
V2 = laplace.get_array()
laplace.init_array()
result3 = laplace.opencl_iterations_norm(5,50,omega=laplace.get_optimal_omega())
V3 = laplace.get_array()
laplace.close()
figure(figsize=(10,5))
plot(result1[0],result1[1],label='V cycle')
plot(result2[0],result2[1],label='Full')
plot(result3[0],result3[1],label='SOR opt')
legend(loc="lower right")
xlabel('niter')
ylabel('norm')
grid()

plotD.pdf

Voici le résultat avec une grille de 2049x2049, avec utilisation de 8 grilles :

import time
laplace = poisson.main.Poisson(4,4,8,1.0,1.0)
laplace.laplacien()
laplace.dirichlet_borders(0.0)
laplace.dirichlet_polygon(6,3,[4],[0],1.0)
t = time.clock()
result1 = laplace.multigrid_Vcycles_norm(1,8,50,2,1,10)
t1=time.clock()-t
V1 = laplace.get_array()
laplace.init_array()
t = time.clock()
result2 = laplace.multigrid_full_norm(3,8,50,2,1)
t2=time.clock()-t
V2 = laplace.get_array()
laplace.init_array()
t = time.clock()
result3 = laplace.iterations_norm(5,50,omega=laplace.get_optimal_omega())
t3=time.clock()-t
V3 = laplace.get_array()
laplace.close()
figure(figsize=(10,5))
plot(result1[0],result1[1],label='V cycle')
plot(result2[0],result2[1],label='Full')
plot(result3[0],result3[1],label='SOR opt')
legend(loc="lower right")
xlabel('niter') 
ylabel('norm')
grid()

plotE.pdf

Temps de calcul :

print([t1,t2,t3])
--> [5.6, 1.5399999999999991, 23.78]

Pour comparaison, la méthode de Gauss-Seidel avec sur-relaxation est réalisée sur CPU. Pour cette méthode, la convergence n'est pas atteinte. On voit bien l'avantage considérable des méthodes multigrilles sur cette taille de maillage.

On remarque que les mêmes paramètres permettent d'obtenir la convergence quelque soit la taille de la grille. Le nombre d'itérations (équivalent grille fine) nécessaires est indépendant de la taille de la grille. Si Q est le nombre de points de la grille la plus fine, le nombre d'itérations ne dépend pas de Q donc la complexité globale est en Q. Par contraste, la méthode de Gauss-Seidel avec sur-relaxation optimale nécessite un nombre d'itérations en Q^1/2, ce qui fait une complexité globale Q^3/2.

Références

[1] W.L. Briggs, V.E. Henson, S.F. McCornick, A multigrid tutorial, (SIAM, 2000)