Cyberlog-corp
Pour une Bretagne libre...

Accueil > Recherche et développement > Parallélisation d’algorithmes > Mécanique des fluides sur multiprocesseur (Bernard Couapel)

Mécanique des fluides sur multiprocesseur (Bernard Couapel)

samedi 2 février 2013, par bc

Toutes les versions de cet article : [Deutsch] [français]

Ce rapport décrit l’implémentation d’une méthode utilisant des volumes finis pour la résolution des équations de Navier-Stokes sur le système multiprocesseur DIRMU. L’exemple testé est le courant laminaire incompressible sur une marche descendante, afin de le comparer à sa résolution sur un ordinateur séquentiel.

La méthode de calcul utilise l’algorithme SIMPLE dans lequel on travaille sur des systèmes d’équations linéaires itératifs pour les deux composantes de la vitesse et une correction des pressions. L’algorithme s’arrête lorsque les équations de masse et de quantité de mouvement sont vérifiées simultanément.

Résolution de pb de mécanique des fluides sur DIRMU - Bernard Couapel
Résolution de problème à deux dimensions de mécanique des fluides sur le multiprocesseur DIRMU
de
B. Couapel*, G. Scheuerer**, J. Volkert***

 


Juillet1987
Version PDF

 

* Institut de Formation Supérieure en

Informatique et Communication

Université de Rennes I

Avenue du Général Leclerc

F—35042 Rennes

** Lehrstuhl für Strömungsmechanik

Universität Erlanqen-Nürnberg

Egerlandstr. 13

D—8520 Erlangen

*** Institut für Mathematische Maschinen

und Datenverarbeitung III

Universität Erlangen-NUrnberq

Martensstr. 3

D—8520 Erlanqen

 

Résumé

Ce rapport décrit l’implémentation d’une méthode utilisant des volumes finis pour la résolution des équations de Navier-Stokes sur le système multiprocesseur DIRMU. L’exemple testé est le courant laminaire incompressible sur une marche descendante, afin de le comparer à sa résolution sur un ordinateur séquentiel.

La méthode de calcul utilise l’algorithme SIMPLE dans lequel on travaille sur des systèmes d’équations linéaires itératifs pour les deux composantes de la vitesse et une correction des pressions. L’algorithme s’arrête lorsque les équations de masse et de quantité de mouvement sont vérifiées simultanément.

 

La configuration du calculateur est choisie de telle sorte que la zone calculée soit divisée en plusieurs bandes. Chaque processeur travaille sur une bande. Théoriquement, avec une telle organisation, on doit obtenir un speedup proportionnel au nombre de processeurs utilisés. Les mesures faites dans le cadre de ce travail montrent que des valeurs proches de l’optimal sont obtenues pour des maillages fins. Les pertes dues aux communications entre processeurs vont de 5 à 20 % selon le maillage et le nombre de processeurs. Les améliorations de la méthode de parallélisation doivent surtout porter sur l’algorithme de résolution des systèmes d’équations linéaires. La méthode ‘ligne par ligne’ est utilisée dans ce travail. Ultérieurement, la recherche sera plus axée sur la parallélisation d’algorithmes plus efficaces.

 

 

PREFACE

Ce travail a été réalisé dans le cadre du domaine de recherche défini par la fondation Volkswagen "Solutions d’équations de base de la mécanique des fluides adaptées aux architectures modernes des calculateurs".

 

Les auteurs remercient le Prof.Dr.W.Händler, Prof.Dr.A.Bode et Prof.Dr.F.Durst pour leur soutien et l’intérêt constant qu’ils ont porté & ce travail. Mme Dipl.Ing.M.Barcus, Mr Dr.I.1.Peric ainsi que Prof.Dr.J.H.Ferziger sont aussi remerciés pour leurs nombreuses discussions et leur aide à la finition de ce travail.

 

 

1. Introduction

 

1.1 Présentation du problème

Les méthodes modernes de calcul de mécanique des fluides s’appuient sur la résolution des équations de Navier-Stokes. Dans le cas laminaire, i.e. avec des nombres de Reynolds petits, celles-ci sont résolues de façon directe.

Pour les courants & turbulences intéressants sur le plan technique, on observe en général des formes de ces équations en relation avec un modèle de turbulence (Bradshaw 1981).

Pour les deux types de courants, il faut intégrer un système d’équations différentielles partielles du deuxième ordre, non linéaires et couplées entre elles.

Pour les courants à recirculation, cette solution consomme beaucoup de temps de calcul et de place mémoire, même pour des géométries à deux dimensions.

Pour les courants à trois dimensions qui sont intéressants à étudier sur le plan pratique, ces besoins sont encore beaucoup plus importants.

Les calculs de configurations techniques complètes telles que avions ou automobiles nécessitent des ordinateurs d’une capacité mémoire d’un milliard de mots et des performances de l’ordre de Giga-Flops.

 

Ces besoins importants entraînent que, notamment aux U.S.A., la mécanique des fluides ait un rôle de plus en plus important dans le développement de calculateurs & hautes performances (Kutler 1985).

Jusqu’à maintenant il s’agissait de calculateurs vectoriels à processeurs puissants, mais moins de parallélisme. Ces dernières années, on a découvert que des progrès importants ne pouvaient plus être réalisés avec de telles architectures sur le plan de la vitesse de calcul et de la capacité mémoire.

C’est pourquoi on s’oriente vers des calculateurs à haut degré de parallélisme i.e. multiprocesseurs dans lesquels plusieurs processeurs travaillent simultanément.

Théoriquement, on peut obtenir des gains en temps de calcul et en place mémoire proportionnels au nombre de processeurs utilisés.

 

Ce concept présente les avantages d’une grande flexibilité et de possibilités d’extension. Les multiprocesseurs offrent le moyen de résoudre d’importants problèmes de la mécanique des fluides qui sont impossibles & travailler sur les machines disponibles actuellement.

Le développement des multiprocesseurs n’ayant débuté que depuis quelques années, il n’y a encore que peu de machines comprenant un grand nombre de processeurs.

Cela entraîne qu’il n’y ait encore que peu d’utilisations dans le domaine de la mécanique des fluides.

Le paragraphe suivant traitera de l’état de la recherche dans ce domaine.

 

1.2 Etat de la recherche

 

Depuis quelques années, les multiprocesseurs comprenant un nombre restreint de processeurs sont techniquement réalisables et disponibles commercialement, Les systèmes ont actuellement en général 2 à 4 processeurs.

De tels multiprocesseurs sont utilisés dans le domaine commercial pour obtenir une amélioration du rendement (ex. I.B.M. 3801). La tendance est aussi aux multiprocesseurs dans le domaine des calculateurs scientifiques à hautes performances (ex CRAY-XMP, CRAY 2, ETA 10).

 

La plupart des machine pourvues d’un plus grand nombre de processeurs qui sont commercialisées ont des liaisons relativement faibles entre les différents éléments.

Un exemple est le calculateur ALLIANT qui comprend 8 processeurs sur un même bus.

Les différentes réalisations de 1’HYPERCUBE sont une exception. Elles sont disponibles depuis peu sur le marché et comprennent des systèmes de liaisons pour 2~ processeurs.

Les multiprocesseurs à système de liaison d’un autre ordre de grandeur (couplage par mémoires multiports, réseaux à plusieurs niveau) sont encore & l’état de recherche.

 

A l’Institut pour Machines Mathématiques et Informatique (IKMD> de l’Université d’Erlangen-Nürnberg on fait des recherches depuis le milieu des années 1970 sur des systèmes couplés par mémoires multiports.

Actuellement, on utilise le multiprocesseur modulaire DIRMU 25 (Händler 1985>. A partir d’éléments de base comprenant un processeur, un coprocesseur et de la mémoire, on peut facilement configurer différentes architectures jusqu’à 26 processeurs (ex anneaux, champs, hypercubes etc). Ce système sert & développer et mesurer les performances de structures et de programmes d’application parallèles.

 

Sur ce système quelques problèmes de mécanique des fluides ont déjà été travaillés avec succès.

Par exemple la résolution de l’équation de Poisson au moyen de différences finies et la méthode de Gau~3-Seidel comme algorithme de résolution, ainsi que la résolution des équations de Navier-Stokes avec une méthode comprenant plusieurs maillages (Multigrid) (Geus 1985).

Dans les deux cas, le gain en temps mesuré est presque linéaire par rapport au nombre de processeurs utilisés.

 

Dans le cadre du projet SUPRENUM du ministère fédéral pour la recherche et la technologie, Thole (1985) a travaillé l’équation de Poisson à trois dimensions sur un calculateur Hypercube. Les processeurs étaient organisés sous la forme d’un champ pouvant comprendre jusqu’à 32 éléments.

La résolution à été réalisée au moyen de multigrid avec des cycles V.

Le travail sur un maillage grossier et les temps de communication importants ont donné comme résultat un gain de temps de l’ordre de 50%.

 

 

 

 

1.3 But du travail

 

Dans ce travail, nous décrirons une application sur DIRMU 25 concernant la résolution des équations de Navier—Stokes avec une méthode à volumes finis. Comme exemple test, nous prendrons le cas d’un flot à deux dimensions, laminaire et incompressible sur une marche descendante. Les résultats de la version parallèle du programme doivent être comparés avec ceux de la version séquentielle. Le but de ce travail est de mesurer quels gains en temps de calcul sont obtenus lorsqu’on calcule le flot sur un multiprocesseur. En même temps on montrera les points faibles et problèmes dans l’implémentation parallèle ainsi que les propositions de solution.

 

 

1.4 Contenu du rapport

Le chapitre 2 décrira la géométrie du problème traitée et le système d’équations utilisé pour sa solution ainsi que les conditions de départ.

Le chapitre 3 consiste en une discussion sur les méthodes de calcul au moyen de volumes finis.

Le chapitre 4 est consacré au multiprocesseur DIRMU. Les principes de base ainsi que les composants matériels et logiciels de ce système sont expliqués.

Le chapitre 5 présente l’implémentation parallèle de l’application numérique. Les paragraphes 5.1 et 5.2 décrivent la configuration du sytème choisie et la différentiation des différents type de variables. Le déroulement du programme parallèle est décrit dans le paragraphe 5.3. Le paragraphe 5.4 traite des aspects importants de la synchronisation entre les processeurs du calculateur.

Les résultats se trouvent chapitre 6 sous forme de tableaux et diagrammes. Il sont discutés et critiqués.

Les conclusions se trouvent dans le chapitre 7 qui termine ce rapport.

 

2. Modèle mathématique

 

2.1 Géométrie et système de coordonnées

La figure 2.1 présente la géométrie et le système de coordonnées de l’exemple-test étudié: le courant à deux dimensions, laminaire sur une marche descendante. La zone d’intégration a une longueur L=0,25 ni; la hauteur du canal à la sortie est H=0,05 m. La hauteur de la marche est h=H/2.

Fig. 2.1 Géométrie de l’exemple test et système de coordonnées.

 

 

2.2 Equations de Navier-Stokes

 

Le courant présenté figure 2.1 peut se décrire par les conservations physiques de la masse et de la quantité de mouvement. Si l’on considère un courant stationnaire, à deux dimensions, laminaire et de densité constante, les équations suivantes pour un fluide Newtonien (Bird , 1960) sont valables:

 

Conservation de la masse:

 

 

Conservation de la quantité de mouvement en x

Conservation de la quantité de mouvement en y

 

Dans ces équations, U et V sont les vitesses dans les directions x et y, P est la pression statique. Les grandeurs p et ja sont la densité et la viscosité dynamique du fluide.

Les équations (2.2-1) à (2.2-3) sont appelées équations de Navier-Stokes. Elles constituent un système d’équations différentielles partielles elliptiques couplées. Pour obtenir une solution unique, il est nécessaire de définir des conditions de départ à la périphérie de la zone d’intégration. Celles-ci sont données dans le paragraphe suivant.

 

2.3 Conditions de départ

Les côtés de la zone d’intégration décrite fig. 2.1 sont formés de parois, d’une zone de flot d’entrée et d’une zone de flot de sortie. Les conditions de départ sont définies comme suit:

Au niveau des parois, les vitesses sont nulles soit:

U = V = 0 (2. 3-1)

A l’entrée du flot, on a un courant qui a un profil de vitesse

 

parabolique suivant la formule:

n est la valeur du point le plus haut de la marche. La vitesse normale V est nulle à l’entrée de la zone d’intégration. Avec les valeurs choisies Umax = 1,5 m/s, p = 1 kg/m3 et -¼ = 0,00375 Pas, nous obtenons un nombre de Reynolds égal à 100 par la formule:

 

Du côté du flot de sortie qui doit se trouver suffisamment éloigné de la zone de courant de retour, on supposera négligeables les variations des vitesses U et V dans la direction x.

 

 

A l’aide des relations données, les équations de Navier-Stokes présentée dans le paragraphe 2.2 sont définies clairement.

 

 

 

3. Méthode de résolution numérique

 

La résolution des équations différentielles de Navier-Stokes est faite au moyen de la méthode de conservation de volumes finis (Finite-Volumen-Methode CAST) décrite par Peric et Scheuerer (1987). Les éléments importants de cette méthode sont les suivants:

 

 

  • Discrétisation de la zone de calcul (Distribution des volumes de contrôle)
  • Discrétisation des équations de conservation
  • Algorithmes de résolution des systèmes correspondants d’équations différentielles linéaires.

 

Ceux-ci seront discutées dans le paragraphe suivant.

 

 

3.1 Distribution des volumes de contrôle

Dans la méthode des volumes finis, la zone d’intégration est divisée en un nombre fini de volumes de contrôle. Les points du maillage se trouvent tous au centre de gravité des volumes de contrôle (fig. 3.1). Jusqu’à maintenant, les équations de Navier-Stokes ont été résolues presque exclusivement (Patankar, 1980) au moyen de distribution décalée des volumes de contrôle (staggered grids). Cette méthode a aussi été utilisée pour ce travail.

Les vitesses U et V sont présentées f ig.3.1. Elles sont décalées par rapport à la pression P dans la zone d’intégration.

Les valeurs des vitesses U sont à la limite du volume de contrôle avec x=const, et les vitesses V à la limite correspondante avec y=const.

Cette distribution présente l’avantage d’éviter des problèmes d’oscillations du champ des pressions. Ce point est discuté en détail par Patankar (1980).

Fig. 3.1: Distribution des volumes de contrôle avec maillages décalés.

 

 

3.2 Discrétisation en volumes finis

 

La discrétisation. est faite selon l’équation de transport généralisée suivante:

 
 

Dans le cas de l’équation de U, -¦U et S-¦=--´P/-´x. L’-¹quation de V peut -ºtre construite de façon analogue. L’équation (3.2-1) sera formellement intégrée sur un volume de contrôle, dans lequel les intégrales du côté gauche peuvent être réécrites en intégrales de surfaces à l’aide du théorème de Gauß. De la figure 3.2, il suit:

 

Cette équation intégro-différentielle est équivalente à l’équation de sortie (3.2-1). Elle décrit l’équilibre entre les flux convectif et diffusifs entrant et sortant des surfaces du volume de contrôle, et les termes—sources qui se trouvent dans le volume de contrôle. Dans la méthode des volumes finis, l’équation (3.2-2) sera approximée à la place de, (3.2-1). Ceci donne une discrétisation conservative qui ne doit présenter aucune perte numérique de masse et de quantité de mouvement.

 

 

Fig. 3.2 : Géométrie d’un volume de contrôle.

 

Dans le premier temps de l’approximation, on exploite l’intégrale de l’équation (3.2-2) et l’équation est linéarisée. Pour cela, nous définissons les hypothèses suivantes:

 

  • Les courants le long des volumes de contrôle ainsi que
  • les termes-sources dans les volumes de contrôle sont constants.

 

Pour linéariser la relation, les courants de masse Cx=pU et Cy=pV seront chaque fois calculés avec les valeurs de 1’ itération précédente. Cela donne l’expression suivante:

 

 

Le deuxième temps de l’approximation consiste & exprimer les valeurs des fonctions -¦e, -¦w, -¦n, -¦s et leurs gradients correspondants sur les surfaces externes des volumes de contr-Äle par rapport aux points du maillage. Pour les premiers, nous utiliserons la méthode UPWIND (Patankar 1980). La valeur -¦ est approchée suivant une fonction en escalier comme fig.3.3, dans laquelle la direction du courant de masse est chaque fois prise en compte. Pour la valeur -¦w de la figure 3.3, on obtient par exemple:

où MAX(a,b) signifie la plus grande valeur entre a et b. Les courants de masse Cx,m avec m w, e, n, s sont déterminés par interpolation linéaire des valeurs des points du maillage.

 


 

Fig. 3.3 : Approximation des termes de correction avec la méthode Upwind.

 

Pour les gradients (-´-¦/-´x)e, (-´-¦/-´x)w, etc... , nous calculons la différence centrale, i.e. le partage de -¦ est approché par des droites. D’après la figure 3.4, nous obtenons:

 

avec les expressions correspondantes pour les autre gradients dans l’équation (3.2-3).

Fig 3.4 Approximation par différence centrale des termes de diffusion.

 

 

Si nous employons cette expression dans l’équation nous obtenons pour chaque volume de contrôle une différentielle implicite de la forme:

 

 

 

 

 

 

 

 

 

 

Il faut remarquer ici que les coefficients de la relation (3.2-7) sont différents du fait de l’utilisation de maillage décalé pour les équations de U et V, puisque par exemple, les courants de masse doivent toujours être utilisés à d’autres endroits.

 

 

3.3 Couplage pression - vitesse

 

Les équations discrétisées de quantité de mouvement sont les suivantes:

 

 

 

 

Celles-ci pourraient être résolues avec un algorithme approprié selon U et V, si la pression P était connue. Comme la pression n’a pas d’équation de transport mais qu’il n’y a qu’une condition de compatibilité entre les vitesse (2.2-1), celle-ci doit être recherchée par une méthode spéciale dans le cas de fluide incompressible qui nous intéresse.

 

Four cela, nous utilisons l’algorithme SIMPLE (Semi-Implicit-Method for Pressure-Linked Equations) développé par Patankar et Spalding (1972). .Les différentes étapes de cet algorithme sont:

 

 

1) Evaluation d’un champ de pression EP*], ex. EF*J0

 

2) Résolution des équations discrétisées de quantité de mouvement (3.3-1) et (3.3-2) avec le champ de pression évalué. Ceci nous donne les champs de vitesse [U*] et [V*].

 

  1. Comme [P*] est inexact, les vitesses calculées [U*] et [V*] ne remplissent pas l’équation de continuité.

 

4) Les correction de pression et de vitesse IF’ J , EU’ J seront faite de sorte que les valeurs corrigées

[U**] = [U*] + [U'] (3.3-3a)

[V**] = [V*] + [V'] (3.3-3b)

[P**] = [P*] + [P'] (3.3-3b)

correspondent à l’équation de continuité. Après ces corrections, les équations discrétisées de quantité de mouvement ne sont plus exactes.

5) Les étapes 2) - 4) sont répétées jusqu’à ce que les équations de quantité de mouvement et de continuité soient simultanément exactes.

Ensuite nous déduisons les relations correspondantes. La discrétisation de l’équation de continuité en un volume de contrôle de pression donne:

Dans un maillage décalé, les courants de masse peuvent être déterminés sans interpolation puisqu’ils sont mémorisés au niveau des surfaces externes des volumes de contrôle de pression. D’après les champs de vitesses EU*J et IV*J nous obtenons une 'source de masse' Sm suivant:

Celle-ci peut être annulée si l’on considère les vitesses corrigées par U’ et V’. L’équation définissant ces corrections est donc:

La soustraction de cette relation à l’équation (3.3-5) donne une expression des corrections des vitesses:

D’un point de vue purement physique, les corrections de vitesses ne peuvent avoir lieu que si les pressions sont modifiées de la valeur correspondant à P’ . La corrélation entre les corrections de pressions et de vitesses se déduit des équations discrétisées de la quantité de mouvement. Après la soustraction des équations définies par les champs [U*] et [P*] des équations de quantité de mouvement exactes (3.3-1) et (3.3-2) on obtient une équation de la forme:

 

Dans la méthode SIMPLE, l’influence implicite des corrections de vitesse sur les points voisins de Ui,j est négligée, ce qui donne l’équation semi-implicite

Des expressions analogues pour les autres corrections de vitesse peuvent être déduites dans l’équation (3.3—7). La substitution de ces expressions donne l’équation de correction de pression

L’équation de correction de pression à la même forme que les équations discrétisées de quantité de mouvement et peut être résolue suivant la même méthode. Avant cela, il faut spécifier les conditions correspondantes à la périphérie. Ceci est décrit dans la suite de ce rapport. Les détails sont donnés par Pakantar (1980).

 

Dans l’exemple de courant sur une marche descendante, les vitesses normales à la périphérie sont connues. Sur les parois elles sont nulles et sur les plans d’entrée et de sortie, elles sont connues. Ceci peut être utilisé pour la formulation des conditions de correction de pression à la périphérie.

Les corrections de vitesse sont nulles au niveau des parois, ce qui après reformulation conduit à des coefficients correspondants nuls pour l’équation de correction des pressions. Ceci correspond physiquement à une condition de Neumann pour la correction des pressions P’ à la périphérie, dans laquelle les gradients normaux aux parois de la zone d’intégration sont nuls.

 

3.4 Méthode de sous relaxation

 

D’après les hypothèses qui ont mené à l’équation (3.3-9), les corrections de pression qui en résultent sont un peu trop importantes. Ceci peut mener à une divergence de l’algorithme SIMPLE si l’on ne fait pas une sous relaxation. C’est pourquoi dans le programme CAST, la pression n’est pas corrigée au moyen de l’équation (3.3-3c) mais par:

[P**] = [P*] + -±P [P‘] (3.4-1)

P est le facteur de sous relaxation pour la pression.

Les vitesses U et V seront sous relaxées comme suit. A gauche de l’équation algébrique (3.2-7) pour la variable générale -¦=(U,V), la valeur -¦ de l’it-¹ration précédente sera ajoutée et soustraite de la façon suivante:

L'index m court sur les point voisins E, W, K, et S. Le terme qui se trouve entre parenthèses correspond à la variation de -¦ pour une itération. Celui-ci est sous relaxé par le facteur -±-¦ dans le programme, soit:

Nous obtenons une bonne convergence de la résolution des équations de Navier-Stokes lorsque (Peric 1987):

Dans ce travail les facteurs de sous relaxation suivants seront utilisés:

 

 

3.5 Algorithme de résolution

 

Les équations (3.3-1), (3.3-2) et (3.3-10) présentent des systèmes linéaires et algébriques pour le calcul des champs [U], [V] et [P']. Elles peuvent s’écrire sous forme de matrices. La matrice de coefficients [A] comporte, dans le cas de la discrétisation UPWIND, cinq éléments différents de zéro, positifs par ligne et à dominante diagonale. Les méthodes de résolution de tels systèmes étant très coûteuses, nous utilisons des méthodes itératives habituelles. Nous allons présenter deux de ces méthodes qui ont été utilisées dans ce travail, la méthode Gau-²-Seidel paragraphe(3.5.2) et la méthode itérative ligne par ligne paragraphe (3.5.1).

 

 

3.5.1 Méthode itérative ligne par ligne

La matrice de coefficients [A] peut, comme présenté fig. 3.5a, être écrite sous forme d’un système tridiagonal. Les sous matrices [Ap] sont tridiagonales. [As] et [An] sont des matrices diagonales. Chacun des vecteurs (-¦j) comporte une ligne o-É y=const dans le maillage num-¹rique.

La transformation du système d’équations présenté fig. 3.5b constitue la base du processus itératif de la méthode ligne par ligne. A gauche, il n’y a que des matrices tridiagonales qui peuvent être résolues efficacement par l’algorithme de Thomas. Le processus itératif fonctionne comme suit:

La solution est d’abord calculée sur la deuxième ligne de volumes de contrôle, i.e. le vecteur [-¦2]. Pour cela, le produit [An][-¦3] de la partie droite du syst-¸me (fig. 3.5) est nécessaire. Comme [-¦3] n’est pas encore disponible, on utilise au départ des valeur évaluées. Après que le vecteur [-¦2] ait été calculé avec l’algorithme de Thomas (i.e. avec une méthode directe), on passe au vecteur [-¦4] puis plus généralement au vecteur [-¦j ].

Ici, les produits [As][-¦j-1] et [An][-¦j+1] sont n-¹cessaires. Des valeurs de d-¹part évaluées seront aussi utilisées.

L’ensemble de la zone de résolution est ainsi balayé avec à chaque fois des valeurs évaluées ou de la dernière itération (fig. 3.5b) pour la partie droite de l’équation.

Au deuxième passage sur le champ d’intégration, le calcul commence sur la troisième ligne.

Les produits [As][-¦j-1] et [An][-¦j+1] n-¹cessaires pour le calcul seront ceux calcul-¹s au premier passage. puis on passe deux lignes plus loin etc, jusqu’à ce que la zone de travail soit balayée une deuxième fois.

Cela nous donne un nouveau vecteur de solution complet [-¦].

La procédure entière est répétée jusqu’à ce que le critère de convergence ou nombre d’itérations (sweeps) donné soit atteint.

Dans l’algorithme SIMPLE, il n’est pas nécessaire à chaque fois de résoudre chaque équation jusqu’à la convergence, puisque après la résolution, les valeurs de la fonction obtenues sont corrigées et les matrices de coefficients sont recalculées.

Il faut cependant être sûr que les équations U, V et notamment l’équation de. .correction de pression, convergent, sinon l’ensemble de l'algorithme divergera après quelques itérations.

Dans ce travail, le nombre d’itérations (Sweeps), i.e. le nombre de passages sur les matrices de coefficients, est prédéfini. A chaque itération, il y a 5 Sweeps pour les équations U et V, et 40 pour l’équation de correction des pressions.

Les équations U et V sont plus simples à résoudre, puisque à chaque fois, on peut utiliser comme valeurs évaluées celles qui se trouvent dans les champs après 1’itération précédente. D’autre part, il y a des conditions Diriclet sur trois côtés de la zone de calcul.

Ceci n’est pas valable pour la correction des pressions; le calcul est à chaque fois démarré avec [P’]=0 à défaut de meilleures valeurs. D’autre part, pour l’équation de pression, les côtés de la zone remplissent la condition de Neumann, ce qui défavorise la convergence.

 

3.5.2 Méthode Gau-²-Seidel

La méthode Gau-²-Seidel itère sur les points. Le calcul se fait selon

Dans ce travail, nous utilisons la méthode Red-Black. La zone de calcul est donc balayée deux fois, comme pour la méthode ligne par ligne.

Au premier passage, les nouvelles valeurs de -¦ seront calcul-¹es

La zone de travail est représentée comme un échiquier. A chaque fois on saute le point voisin (on ne calcule que les points 'rouges' )

Du côté droit de l’équation (3.5-1) les valeurs sont évaluées pour le premier sweep. Pour les passages suivants, les valeurs sont celles calculées au sweep précédent.

Au deuxième passage, ce sont les points ‘noirs’ qui sont calculés. Four cela on utilise dans la partie droite de l’équation (3.5-1) les valeurs -¦'ij calculées au premier passage.

L’algorithme continue comme décrit dans le paragraphe précédent jusqu’à ce que le critère de convergence ou le nombre d’itérations prédéfini soit atteint.

 

 

3.5.3 Evaluation des algorithmes de résolution

Les algorithmes de résolution sont une partie importante du programme car ce sont eux qui demandent le plus de temps de calcul.

Il est donc indispensable d’utiliser des méthodes les plus efficaces possibles. Il faut cependant trouver un compromis pour les calculateurs parallèles.

La méthode Gau-²-Seidel décrite dans le paragraphe 3.5.2 est très efficacement parallélisable. Cependant, elle devient excessivement lente avec l’augmentation du nombre de points de maillage ( surtout, pour des problèmes à conditions périphérique de Neumann comme pour l’équation de correction des pressions ). Le gain de temps obtenu par la parallélisation peut même dans le cas extrème être annulé par la mauvaise efficacité de cet algorithme.

Ceci explique qu’après son utilisation au début de ce travail il ait été abandonné.

L’efficacité de l’algorithme ligne par ligne est améliorée lorsque la résolution directe du système tridiagonal ne se fait pas seulement suivant la direction y=const. mais de façon alternée y=const et x=const. Mais dans ce cas, son degré de parallélisation diminue. Cette possibilité n’a donc pas été utilisée dans ce travail.

Dans la première version du programme CAST (Peric et Scheuerer 1987), la méthode de Stone (1968) a été utilisée pour l’équation de correction des pressions. Celle-ci se base sur une décomposition incomplète. L’efficacité de cette méthode est très supérieure pour la résolution de l’équation de correction des pressions aux deux algorithmes présentés dans ce rapport (Peric 1987 et Barcus 1987). Mais l’algorithme de Stone (1968) n’est pas suffisamment parallélisable. Dans un travail ultérieur, on calculera si une telle méthode n’a pas une efficacité globale meilleure malgré un degré de parallélisation plus faible.

 

3.6 Critère de convergence

Le critère de convergence est le suivant:

Après que les coefficients de l’équation (3.2-7) soient calculés, les résidus pour les volumes de contrôle U et V sont établis à partir de:

Dans le cas de l’équation de continuité, le résidu est la source de masse Sm. Pour toutes les variables, la somme des résidus absolus est calculée et mémorisée dans la grandeur F-¦. Dans ce travail, on utilise pour F-¦ les masses entrantes et les flots de quantit-¹ de mouvement. Il est alors attendu que ces grandeurs soient inf-¹rieures à un critère de convergence prédéfini -µ i.e.:

Pour ce travail, on a choisi -µ=0.005.

 

4. Multiprocesseur DIRMU

(DIstributed Reconfigurable MUltiprocessor kit)

 

4.1 Matériel

Le multiprocesseur DIRMU 25 (Händler 1~85) a été développé à l’université d’Erlangen—Nürnberg pour tester les configurations de multiprocesseurs et la programmation parallèle.

L’idée de départ de DIRMU est de proposer un élément de base permettant de construire un large éventail de configurations de multiprocesseurs adaptées à diverses tâches.

Un module de base est présenté figure 4. 1, Il comprend un module processeur (P-module) et un module mémoire multiport (M-module). Chacun comprend 8 connexions (ports). Les P-ports 1 - 7 peuvent être connectés aux M-ports d’autres modules DIRMU. Les M-modules sont connectés à plusieurs processeurs qui peuvent ainsi avoir directement accès aux données de la mémoire multiport.


Fig. 4.1 : module de base DIRMU.

 

A partir de ces éléments, on peut construire différentes configurations de calculateurs. La figure 4.2 montre plusieurs exemples.

En théorie, on peut combiner autant de processeurs que l’on veut. Cependant, il n’y au plus que huit processeurs directement reliés entre eux, il s’agit donc d’un structure à voisinage limité.

 


Fig. 4.2 Exemples de configurations DIRMU.

 

Un module de base est présenté en détail figure 4.3.

Le P-module contient un microprocesseur 8086 et un coprocesseur 8087 de la firme Intel, une mémoire RAM/ROM privée (actuellement de 320K/16K octets), des interfaces d'Entrées/Sorties (terminal, disque, imprimante, etc...) et 8 P-ports.

Le M-module est une mémoire multiport (64 K octets) comprenant 8 M-ports.

Par le P-port 0 et le M-port 0, chaque processeur peut accéder son propre M-module qui fait partie de son espace mémoire. Les P-ports 1,2,...,7 peuvent être connectés au moyen de câbles aux M-ports des autres unités DIRMU. Dans ce cas, les M-modules sont partagés entre plusieurs processeurs ce qui permet de partager les données se trouvant dans la mémoire multiport.

La figure 4 montre un exemple simple.

La zone de données A est vue par le processeur F0 à partir de l’adresse 1:6000H ( l=segment, 6000H=offset). La même zone de données peut être accédée par P2 à partir de l’adresse 2:6000H et par P1 à partir de l’adresse 7:6000H.

Fig. 4.3 : Structure d'un module de base DIRMU

 


Fig. 4.4 : Adressage dans les M-modules

 

Les adresses du 8086 sont décodées par le contrôleur du P-port. Lorsque le contrôleur du P-port décode le segment i, i=1,2,...,8, il sélecte le M-port i-1 et accède au module correspondant. L’offset reste inchangé.

Les mots de données sont de 16 bits. Le temps d’accès est de 800ns pour la mémoire locale, et de 1,2-¼s pour la mémoire multiports (sans conflits d'accès).

Le temps d’accès plus important pour la mémoire multiports est dû aux mémoires RAM plus lentes et aux retards des contrôles des P- ports et M-ports. Il faut souligner que toutes les mémoires multiports connectées à un processeur ont le même temps d’accès.

Chaque processeur a sa propre horloge. Si plus d’un processeur essaie d’accéder en même temps à la mémoire multiport, il y a un conflit d’accès qui est résolu par le contrôleur de M-port.

Un processeur attend jusqu’à l’obtention de la donnée. La règle est telle qu’un processeur ne peut pas accéder 2 fois à la mémoire de suite si un autre processeur attend d’accéder à la même mémoire.

 

Des mesure ont montré (Maehle 1985) que le temps d’accès plus important pour le M-module, ainsi que les retards dus aux conflits d’accès n’ont qu’une très faible influence sur le temps total d’exécution d’un programme d’application, car seules les données partagées sont stockées dans les M-modules. La majorité des accès mémoire se font sur la zone de code et les données dans la mémoire privée.

Un processeur peut bloquer la mémoire multiports pendant plus d’un cycle mémoire. Cela est nécessaire pour certaines primitives de synchronisation (ex: instruction EXCHANGE) et la consistance des données qui sont plus larges qu’un mot mémoire (ex. variables réelles).

La communication entre processeurs se fait par deux registres dans l’espace d’adressage de la mémoire multiports. Deux processeurs qui partagent le même M-module peuvent s’envoyer des interruptions en positionnant les bits de ces registres.

Ce mécanisme peut être utilisé pour le trafic de messages entre processeurs.

Un autre registre status supporte le diagnostic d’erreur des processeurs voisins. Si une erreur est détectée dans une unité (ex problème d’alimentation, instruction HALT illégale etc.), un bit correspondant est positionné dans le registre status et une interruption est envoyée à tous les voisins.

Le diagnostic du voisinage est un prérequis important pour les configurations DIRMU résistant aux pannes.

Chaque unité DIRMU a sa propre alimentation et se trouve dans un rack 3/4 19" ( voir fig. 4.5). On peut contrôler la connexion des P-modules et M-modules ainsi que la charge des ports sur la face avant par des diodes électro-luminescentes.

 

 

 

 

 


Fig. 4.5 Faces avant et arrière d’un module DIRMU.

 

 

4.2 Système d’exploitation et langage de programmation

 

Nous décrirons dans ce chapitre le sytème fonctionnant sur DIRMU, le DIRMOS (Diagnosis Integrated Reconfiguring Multiprocessor Operating System). Il se greffe sur le système Concurent CP/M.

Le développement d’applications parallèles sur DIRMU est basé sur le concept de rnaître/esclave.

Chaque unité DIRMU équipée d’un disque et d’un terminal (MASTER) peut être utilisée comme unité de développement en mode monoprocesseur. L’exécution d’une application parallèle est supportée par DIRMOS.

Dans les applications parallèles, le MASTER fait toutes les E/S disques et le dialogue avec l’utilisateur. Les autres processeurs coopérant à l’application parallèle sont appelés SLAVES.

Différents utilisateurs peuvent se partager la machine de telle sorte que chaque SLAVE soit assigné à un MASTER à la fois (space sharing).

 

Une des caractéristiques de DIRMOS est l’intégration de la gestion du diagnostic et de la configuration.

Pendant les moments libres, le diagnostic calcule la configuration de processeurs en bon état. Cela comprend les processeurs en bon état et (par un chemin d’accès) les processeurs atteignables ainsi que leurs connexions.

La configuration utilisée par le programme d’application (‘application configuration’) est mappée par DIRMOS (Configuration Finder) sur la machine en bon état.

DIRMOS fournit le numéro de port du voisin logique pour chaque processeur ainsi que sa position dans la configuration appelée par l’application.

D’autre part, DIRMOS comporte un chargeur parallèle qui charge l’application parallèle dans la mémoire privée de chaque processeur participant à celle-ci.

 

 

 

DIRMOS consiste en trois parties principales:

  • La partie résidente COMINT.
  • Le FINDER.
  • Le MANAGER.

COMINT est la partie résidente de DIRMOS. Il est chargé dans la mémoire privée de tous les modules à l’initialisation du système.

En tapant le caractère de commande à la console, l’utilisateur change l’état du processeur correspondant de FREE SLAVE en MASTER et entre dans 1’interpréteur de commandes.

L’application est démarrée en tapant le nom du programme parallèle. Avant que celui-ci ne soit chargé dans les SLAVES, l’utilisateur doit spécifier la configuration pour laquelle le programme a été écrit. Le FINDER approprié est chargé dans le MASTER, cherche la configuration la plus grande de ce type et fournit le nombre maximal de processeurs disponibles. L’utilisateur peut alors choisir la nombre désiré de processeurs pour son application.

Les processeurs réservés pour l’application sont mis en état de SLAVE par le configurateur et le programme est chargé en pipeline à partir du MASTER (fig. 4.6).


Fig. 4.6 : Chargement du programme parallèle.

L’application parallèle est chargée au dessus du COMINT en overlay.

Le MANAGER est un module en librairie qui est linké au programme d’application. Le MANAGER trouve les informations concernant le mapping fait par le FINDER dans la partie résidente (COMINT). Le MANAGER peut ainsi calculer la position logique du processeur dans la configuration de l’application. Le nombre de ports connectés aux voisins logiques est aussi calculé.

Après l’initialisation du MANAGER, l’APPLICATION a accès au numéro de port de ses voisins logiques et peut créer des structures de données partagées en appelant le LOCATOR.

Le LOCATOR est un module en librairie qui calcule les adresses des données partagées.

Le système DIRMOS supporte trois principaux types de configuration:

  • anneau (RING).. y
  • tableau <ARRAY)
  • cube à n dimensions(N-CUBE>

Si l’application demande une autre configuration (ex. une pyramide), il faut utiliser le configurateur spécial. L’utilisateur décrit alors sa propre configuration.

 

Le langage de programmation utilisé est le MODULA2 (Wirth 1985) qui s’apparente au PASCAL.

 

5. Parallélisation de la méthode de résolution

 

Dans ce chapitre, nous décrirons la parallélisation de la méthode de résolution présentée dans le chapitre 3.

Le premier paragraphe présentera la configuration du calculateur utilisée.

Le partage des variables du programme d’application et le déroulement de la méthode parallèle seront décrits dans les paragraphes 5.2 et 5.3.

La synchronisation nécessaire au déroulement parallèle du programme sera présentée paragraphe 5.4.

 

5.1 Configuration du calculateur

Le premier temps de la parallélisation d’une méthode de résolution est le choix d’une configuration de calculateur adaptée.

Pour le problème ici traité d’une zone d’intégration à deux dimensions rectangulaire, nous avons le choix entre deux structures (voir fig. 4.2)

  • anneau
  • champ

 

Du fait de la facilité des problèmes de synchronisation et de la mise en oeuvre aisée de la méthode de résolution ligne par ligne sans alternance d’orientation des lignes, la structure choisie est l’anneau d’une longueur maximale de 26 processeurs (voir fig. 5.1).

La zone de calcul est divisée en bande comme figure 5.2. Chaque bande est travaillée par un processeur. Sur l’anneau de la figure 5.2, le MASTER est noté M et les SLAVES Sj.

 

 

 

 

 

 

Fig. 5.1 : Configuration du multiprocesseur DIRMU

Processeurs Zone de calcul

Fig. 5.2 : Partage des variables

 

5.2 Partage des variables

Le deuxième temps de la parallélisation consiste à partager les variables entre les processeurs.

Les variables dépendantes du problème traité sont les champs de vitesses [U] et [V] ainsi que le champ de pressions [P]. Ce sont des matrices de longueur NX*NY où NX est le nombre de points de maillage dans la direction X et NY le nombre de points dans la direction Y. Leur calcul nécessite l’utilisation de douze autres matrice du même type pour

  • les coefficients des volumes finis (équation (3.2-8)) 5
  • les courants de masse Cx et Cy aux surfaces des volumes de contrôle 2
  • le champ de correction des pressions 1
  • le terme source S 2
  • des champs accessoires pour le calcul des coefficients de l’équation de correction des pressions. 2

L’utilisation de vecteurs de longueur NY pour la mémorisation des courants aux bords des volumes de contrôle par exemple n’entre presque pas en ligne de compte pour le calcul des besoins en mémoire.

Pour le partage des variables, les concepts de variables locales, régionales et globales seront utilisés. Il sont définis comme suit

  • Les variables locales ne sont pas partagées, elles sont en mémoire privée et initialisées dans le programme par le processeur Sj (ou M)
  • Les variables régionales sont partagées par des processeurs voisins. Elles se trouvent en mémoire multiport.
  • Les variables globales sont partagées par tous les processeurs utilisés dans l’application. Elles se trouvent aussi en mémoire multiport.

 

Les matrices de longueur NX*NY seront traitées comme des variables régionales, comme le montre la figure 5.2, chaque processeur mémorise une bande comprenant les données de NI et NJ/M volumes de contrôle. M est le nombre de processeurs, le nombre de volumes de contrôle correspond à NI=NX-2 et NJ=NY-2. Les processeurs 1 et M situés aux extrémité de la zone de travail ont une ligne de maillage en plus.

Les sommes de résidus absolus utilisées pour le critère de convergence, le courant de masse vers l’extérieur de la zone d’intégration et la valeur de référence pour la pression sont traitées comme variables globales.

Pour ce travail, l’ensemble des sous-matrices est mémorisé dans les mémoires multiports. Il serait plus efficace de ne mémoriser en mémoire multiport que les rangées de volumes de contrôles situées aux bords des bandes, le reste des données étant situé en mémoire privée (zone hachurée fig. 5.2). Ainsi, lors du calcul notamment de maillages fins, où contrairement à la figure 5.2 le rapport entre zone locale et zone régionale partagée entre processeurs est élevé, on pourrait remplacer la plupart des accès à la mémoire multiport par des accès à la mémoire locale plus rapide.

Cependant, ceci aurait pour effet d’augmenter sensiblement le travail de programmation et la transparence d’un programme parallèle (la ‘sous-matrice’ devrait être séparée en deux objets local et régional). C’est pourquoi cette possibilité n’a pas été retenue.

La démarche sera la même pour les prochains systèmes multiprocesseurs d’Erlangen. La mémoire multiport sera augmentée par rapport à la mémoire privée, afin que des objets de ce type puissent être mémorisés. La mémoire privée ne contiendra plus que le système d’exploitation et les données locales.

Ce travail comporte une limitation qui est due au caractère de test du programme. Le maître contient l’ensemble des champs de variables et les rassemble à la fin du calcul.

Ceci entraîne que la grandeur maximale du problème traité est limitée à la place mémoire du MASTER. Cette limitation a pour avantage de faciliter les mesure d’efficacité et du speedup.

Lors des mesures, le même problème sera mesuré d’abord sur un seul processeur, puis sur plusieurs, et les temps de calcul seront ensuite comparés.

La limitation peut être évitée par une modification de l’initialisation et de la collecte des résultats afin que pour des applications pratiques on puisse travailler sur des maillages beaucoup plus fins.

 

5.3 Déroulement du programme

 

Le premier temps de la méthode est l’initialisation. Les champs de vitesses et de pression sont initialisés avec des valeurs évaluées sur chaque processeur ainsi que les variables régionales et locales.

Ensuite la méthode de résolution se déroule sur chaque processeur.

D’abord, les coefficients des volumes de contrôle et les termes—source de l’équation de quantité de mouvement en X (3.3-1) sont calculés. Pour cela, on utilise les flots convectifs et diffusifs qui passent par les surfaces des volumes de contrôle. Une première boucle calcule les flots à travers les surfaces ouest des volumes de contrôle pour chaque processeur.

Pour les interpolations au niveau des volumes du bord de la zone de calcul, on prend à chaque fois des valeurs des processeurs voisins.

Dans une deuxième boucle, les processeurs calculent les flots à travers les surfaces Sud, Nord et Est des volumes de contrôle. Les flots à travers les surfaces est sont mémorisés dans un vecteur à une dimension de longueur NY et utilisés comme flot Ouest de la colonne suivante.

Aux frontières des sous-matrices, dans la méthode de volumes finis, il faut seulement faire attention à ce que le flot passant par la surface nord du volume de contrôle le plus haut de la zone Sj soit identique à celui passant par la surface sud du volume de contrôle le plus bas de la zone Sj+1.

Le déroulement décrit est le même pour toutes les autres colonnes de volumes de contrôle.

Le calcul des coefficients est facilement parallélisable puisqu’aucune synchronisation n’est nécessaire.

Après le calcul des coefficients de l’équation (3.2-8), les valeurs de [U] en mémoire sont utilisées dans l’équation de volumes finis (3.3-1) et la somme des résidus absolus est calculée dans toutes les sous-matrices selon l’équation (3.6-1). Les valeurs résultantes pour chaque sous-zone seront mémorisées jusqu’au tour de synchronisation qui suit le calcul de la quantité de mouvement en Y et la correction des pressions.

Au cours de cette méthode, l’algorithme de résolution décrit paragraphe 3.5.1 est appelé.

Les étapes de sa parallélisation seront 4écrites à la fin de ce paragraphe.

Après la résolution du champ [U], les valeurs du champ [V] sont calculées de façon analogue.

Ensuite, on fait une correction du flot de masse à la sortie de la zone de résolution, car puisque cette méthode est conservative, la quantité de masse à la sortie doit être égale à la quantité de masse à l’entrée.

Pour la réalisation de cette condition, on calcule le flot de masse sortant en faisant la somme globale des flux de masse sortant de la dernière rangée de volumes de contrôles soit:

Chaque processeur calcule d’abord la masse sortant de sa zone de travail, pùis après une synchronisation de tous les processeurs (voir paragraphe 5.4) pour le calcul de la somme globale, on obtient la vitesse à la sortie de la zone de travail par

 

Les courants de masse Cx et Cy sont recalculés avec les valeurs des champs de vitesses ainsi que le terme source pour l’équation de correction des pressions (3.3-5). En même temps que l’on calcule le courant de masse et le terme source, on détermine les coefficients de l’équation de correction des pressions.

Les coefficients d’équation de correction des pressions sont symétriques, i.e. à chaque fois, le coefficient sud de la ligne j de volumes de contrôle est le même que le coefficient nord de la ligne j-1.

Lors du calcul des coefficients des volumes de contrôle aux frontières des zones de chaque processeur, il faut soit une synchronisation régionale, soit un calcul spécifique pour ces valeurs (comme réalisé dans ce travail).

Après la détermination des coefficients, le résidu de masse est formé de la somme des valeurs absolues des termes source Sm (équ. (3.3-5)) et distribué à tous les processeurs.

Ensuite, l’algorithme de résolution est appelé.

Il faut ici faire attention à ce que l’équation de correction des pressions soit résolue de façon suffisamment précise, sinon on risque d’avoir des corrections de pression irréelles et de faire diverger l’algorithme.

Après la résolution de l’équation de correction des pressions, les pressions et les vitesses sont corrigées avec la relation (3.3-3). Cela est exécuté parallèlement sur toutes les sous-zones.

Ainsi se termine une itération complète de la méthode. Ensuite, une synchronisation globale permet d’ajouter les valeurs résiduelles des équations de correction de U, V et P pour obtenir des valeurs globales. Ces valeurs sont comparées à des conditions de fin données.

Si la valeur maximale de ces trois résidus est plus grande que la condition de fin, 1’itération décrite est relancée avec les nouvelles valeurs en mémoire.

Sinon le programme se termine.

Le maître envoie aux esclaves un ordre de récupération des données qui retournent leurs sous champs [U], [V] et [P]. Les résultats sont rassemblés et affichés.

Pour terminer ce paragraphe, nous décrirons la parallélisation de la méthode ligne par ligne.

Comme nous l’avons montré dans le paragraphe 3;5.1, cet algorithme se déroule en deux passages. Lors du premier, on utilise des valeurs évaluées pour la première itération, sinon les valeurs de l'itération précédente.

C’est pourquoi le premier passage peut se faire sur tous les processeurs simultanément sans synchronisation.

Par contre une synchronisation régionale est nécessaire lors du deuxième passage puisque le processeur Si a besoin de certains résultats des processeurs Sj-1 et Sj+1 et doit les attendre.

Les résolutions directes se faisant uniquement pour des lignes où y=const, une structure en anneau de processeurs permet un degré relativement important de parallélisation de la méthode ligne par ligne.

 

 

5.3 Synchronisation

 

La synchronisation est un élément primordial dans le développement d’algorithmes parallèles.

Nous distinguons deux catégories:

  • synchronisation régionale
  • synchronisation globale

 

La synchronisation régionale ne joue qu’entre processeurs voisins sur des variables communes. Chaque processeur interrompt le déroulement de son programme jusqu’à ce que les variables de son voisin aient pris les valeurs désirées.

La synchronisation globale est nécessaire par exemple lors du calcul de la somme des résidus absolus.

Elle se déroule comme suit:

  1. Le maître écrit sa somme partielle dans la mémoire de son voisin Si et positionne une variable booléenne a VRAI.
  2. Cela signale & Si que la somme partielle est valable.
  3. Si ajoute sa somme partielle à cette valeur, repositionne la variable booléenne du maître à FAUX et envoie la nouvelle somme partielle à S2 de la même façon etc...
  4. De cette manière, tous les processeurs sont interrogés.
  5. Le maître reçoit un signal du dernier processeur SM et lit la somme globale des résidus absolus dans sa mémoire multiports.

De la même façon, la somme globale des résidus est envoyée à tous les processeurs de l’anneau.

Si la somme des résidus est inférieure à une valeur donnée, le programme s’arrête.

Le type de synchronisation globale choisi présente le désavantage d’une perte de temps importante pour des anneaux comprenant de nombreux processeurs. Pour le cas de la machine DIRMU à 26 processeurs, ces pertes de temps sont encore acceptables

 

 

 

6. RESULTATS

 

6.1 Courants

La figure 6. 1 montre le vecteurs de vitesse calculés et les lignes de courant.

Les vitesses et pressions calculées par ce programme parallèle ont été comparées aux résultats d’un programme séquentiel fonctionnant sur le calculateur CDC Cyber 855 du centre de calcul régional d’Erlangen.

Les valeurs sont les mêmes à la précision de calcul près.

 

 

Figure 6.1 Vecteurs de vitesse et lignes de courants.

 

 

6.2 Accélération grâce à la paralléllisation

 Les mesures de temps de calcul ont été effectués pour deux maillages différents.

Le premier comprend 20x10 volumes de contrôle dans les direction respectives x et y.

Le deuxième comprend 20x46 volumes de contrôle.

Les temps de calcul, accélération grâce à la parallélisation (speedup) et l'efficacité sont regroupés dans les tableaux 6-1 et 6-2.

D’autre part, les speedup des deux tests sont présentés sous forme de fonction sur les figures 6.2 et 6.3.

Matrice 22x12 points.

RESOR: U=2,3E-6 V=2,24E-7 P=4,19E-3

Nombre de Processeurs

Nombre d'iterations

Durée

(s)

Speedup

Efficacité

(%)

1

29

446.700

1.00

100

2

29

265.523

1.68

84

3

29

187.245

2.39

80

4

29

179.876

2.48

62

6

29

98.981

4.51

75

Tableau 6-1 Resultats pour 20x10 volumes de contrôle.

 

Le speedup est défini comme le temps de calcul avec un processeur divisé par le temps de calcul pour M processeurs. L’efficacité est la relation entre le speedup réel par rapport au speedup théorique. Le speedup idéal est égal au nombre de processeurs.

Les deux tableaux montrent que pour des maillages plus fins, on obtient un speedup et une efficacité nettement meilleurs que pour un maillage grossier.

Matrice 22x48 points.

RESOR: U8,27E-4 V=3,35E-5 P=4,79E-3

Nombre de Processeurs

Nombre d'iterations

Durée

(s)

Speedup

Efficacité

(%)

2

235

8619.997

2.00

100.0

3

235

5828.656

2.96

98.6

4

235

4399.256

3.92

98.0

6

235

2962.094

5.82

97.0

8

235

2242.709

7.69

96.1

12

235

1523.944

11.31

94.3

16

235

1461.889

11.79

73.7

24

235

807.000

21.36

89.0

Tableau 6-2 Résultats pour 20x46 volumes de contrôle.

 

Dans le cas de maillage plus fin, on obtient une efficacité supérieure à 90% jusqu’à 12 processeurs. On a donc un plus grand degré de parallélisme.

Lorsqu’on dépasse 12 processeurs, l’efficacité tombe au dessous de 90%. La raison est que les bandes travaillées par chaque processeur sont de plus en plus faibles avec l’augmentation du nombre de processeurs, et donc le temps de Synchronisation devient plus important.

Cependant, même avec 24 processeurs, on gagne un facteur 21 sur le temps de calcul avec un seul processeur.

Pour le maillage grossier, l’efficacité reste au dessous de 90%. Ceci est aussi dû au temps de synchronisation important par rapport au calcul.

Le décrochage de la courbe du speedup ( ex. lorsqu’on passe de 12 à 16 processeurs fig. 6.3) est dû aux inégalités de charge entre les différents processeurs. Tous les processeurs n’ont pas le même nombre de volumes de contrôle à calculer, ce qui entraîne que les processeurs les moins chargés doivent attendre les processeurs les plus chargés.

La meilleure efficacité serait obtenue si le maximum de place mémoire de chaque processeur était utilisé. Le temps de synchronisation diminuerait par rapport au temps de calcul ce qui ferait augmenter l’efficacité.

 

 

Fig. 6.2 : Speedup obetnu pour 20x10 volumes de contrôle

 

 

 

 

Fig. 6.3 : Speedup obetnu pour 20x10 volumes de contrôle

 

 

6.3 Discussion et propositions d’améliorations

Dans ce paragraphe nous ferons quelques remarques sur l’optimisation de la méthode.

 

6.3.1 Algorithmique

 

L’algorithme de résolution est la partie qui demande le plus de temps de calcul. Dans la résolution séquentielle, on utilise 1’algorithme de Stone (1968) qui est efficace et basé sur un découpage. incomplet. Cette méthode présente le désavantage de ne pas être :parallélisable.

Dans le ‘programme parallèle, nous avons utilisé l’algorithme point par point Gau-²-Seidel et la m-¹thode ligne par ligne qui sont très parallélisables mais peu efficaces pour la résolution de grandes matrices. Ceci est un problème, surtout pour la résolution de l’équation de correction des pressions qui remplit les conditions de Neumann pour les bords et dont la convergence avec l’algorithme de Gau-²-Seidel est tr-¸s mauvaise. L’utilisation de la m-¹thode ligne par ligne n’apporte pas d’amélioration déterminante.

Il sera donc important de développer et tester ultérieurement des algorithmes parallélisables performants. Il faudra aussi vérifier dans quelle mesure l’utilisation de méthodes travaillant sur des points couplés est valable. Ce sont des méthodes itératives sur des points qui résolvent pour chaque point un système pour les composantes de la vitesse et la pression (Vanka 1986). Ces méthodes de relaxation peuvent être accélérées en faisant des cycles avec plusieurs maillages (multigrid).

Le calcul des coefficients et la méthode ligne par ligne sont très parallélisables sur une configuration en anneau. Une synchronisation n’est nécessaire qu’aux limites Nord-Sud entre les sous-zones.

Dans le cas d’une configuration en champ, il faut ajouter à la synchronisation Nord-Sud une synchronisation Est-Ouest aux limites des sous-zones. Ceci entraîne une augmentation du temps de synchronisation. D’autre part, le nombre de volumes de contrôle calculés par processeur diminue quand on passe de l’anneau au champ.

Le choix entre l’anneau et le champ dépend de la grosseur du problème et du nombre de processeurs disponibles.

Pour un grand nombre de processeurs, on peut choisir une structure de champ plus intéressante que l’anneau sur le plan de la synchronisation (bandes étroites, chaque synchronisation demande n étapes au lieu de sqrt(n) pour le champ).

Les pertes de temps sont dues à la synchronisation. Les synchronisations globales demandent deux tours d’anneau le premier pour le calcul de la variable et le deuxième pour la diffusion de sa valeur globale à tous les processeurs. Les synchronisations doivent donc être réduites au minimum. Ceci a été fait pour ce travail.

Une synchronisation peut parfois être évitée grâce à un calcul d’une variable par deux processeurs voisins. Ceci est plus rapide même au prix d’un calcul un peu plus long, qu’une synchronisation.

 

6.3.2 Matériel

La machine DIRMU a été conçue dans les années 70 et travaille avec un processeur 8086/8087 à la vitesse de 5 MHz.

Il est clair que l’utilisation de processeurs 32 bits de la nouvelle génération augmenterait sensiblement la vitesse de calcul.

D’autre part, le développement de coprocesseurs microprogrammables spécialisés dans le calcul d’une partie souvent utilisée de l’algorithme (ex. algorithme de résolution Gau-²-Seidel) permettrait un gain important en temps de calcul.

 

 

7. Conclusion

 

 

Ce travail a consisté en la parallélisation d’une méthode à volumes finis des équations de Navier-Stokes à deux dimensions et incompressibles sur le multiprocesseur DIRMU.

La configuration choisie est l’anneau de processeurs. La zone de calcul est divisée en bandes travaillées chacune par un processeur.

L’exemple test est le courant laminaire sur une marche descendante. Les résultats ont été comparés à ceux obtenus sur une machine séquentielle.

On peut tirer les conclusions suivantes du travail réalisé:

 

  • La méthode à volumes finis est très efficacement paralléllisable. Le travail sur une zone de 20x46 volumes de contrôle avec 24 processeurs permet de gagner un facteur 21,4 ce qui représente une efficacité de 89%.Une diminution du nombre de processeurs fait augmenter l’efficacité (97% pour 6 processeurs) puisque le temps de synchronisation diminue.Le caractère de test du programme n’a pas permis le travail sur des maillages fins.
  • L’efficacité importante obtenue ne doit pas cacher le fait que l’utilisation des méthodes Gau-²-Seidel et ligne par ligne demandent beaucoup plus d’itérations (sweeps) que la méthode séquentielle de Stone (1968) qui présente le désavantage d’être difficilement parallélisable. L’utilisation d’algorithmes parallèles diminue l’efficacité de la méthode. Il faut calculer si le gain obtenu grâce à la parallélisation compense cette perte d’efficacité. Sinon, la parallélisation de la méthode ne présente aucun intérêt.

 

A partir de ces résultats, les travaux suivants devront être réalisés:

  1. Test de différentes méthodes de résolution sur le plan de la parallélisation mais aussi de leurs performances.
  2. Ajout aux algorithmes de solution d’une méthode par plusieurs niveaux de maillage (multigrid) pour la résolution des équations de Navier-Stokes. Cette méthode qui accé1ère la convergence doit être absolument utilisée pour avoir dans l’avenir des temps de calcul acceptables sur des machines parallèles pour des maillages fins.
  3. Comparaison entre des méthodes couplées ou non pour la résolution des équations de Navier—Stokes sur des calculateurs parallèles. L’algorithme SIMPLE utilisé est une méthode non couplée: trois matrices sont successivement résolues pour les deux vitesses et la correction des pressions. Dans la méthode couplée, on résout une matrice par volume d~ contrôle, i.e. pour chaque volume de contrôle, on résout un système d’équations. Il y a dans cette méthode une part plus importante de temps de calcul par rapport à la synchronisation, ce qui pourrait améliorer l’efficacité de la parallélisation. Ces arguments doivent cependant être confirmés par des essais.

 

 

 

8 Bibliographie

 

 

Barcus, M., 1987, "Berechnung zweidimensionaler Strömungsprobleme mit Mehrgitterverfahren", Diplomarbeit, Lehrstuhl für Strömungsmechanik, Universität Erlangen-Nürnberg.

Bird, R.B., Stewart, W.E., Lightfoot, E.N., 1960, ‘Transport Phenomena", J. Wiley & Sons Inc., New York.

Bradshaw, P., Cebeci, T., Whitelaw, J.H., 1981, "Engineering Calculation Methods for Turbulent Flow", Academic Press, London.

Geus, L., 1985, "Parallelisierung eines Mehrgitterverfahrens für die Navier-Stokes-Gleichungen auf EGPA-Systemen", Arbeitsberichte des Inst. f. Mathem. Maschinen und Datenverarbeitung, Universität Erlangen-Nürnberg, Band 18, Nummer 3.

Händler, W., Maehle, E., Wirl, K., 1985, "The DIRI4U Testbed for High-Performance Multiprocessor Configurations", Proc. of the First Intern. Conf. on Supercomputing Systems, St. Petersburg, FL, pp. 468—475.

Kutler, P., 1985, "A Perspective of Theoretical and Applied Computational Fluid Dynamics", AIAA-J., Vol. 23, pp. 328-341.

Maehle, E., Wirl, K., Japel, D., 1985, "Experiments with Parallel Programs on the DIRMU Multiprocessor", Parallel Computing, Berlin, pp. 515-520.

Patankar, S.V., 1980, "Numerical Heat Transfer and Fluid Flow", Hemisphere Publ. Corp., Washington.

Patankar, S.V., Spalding, D.B., 1972, "A Calculation Procedure for Heat, Mass and Momentum Transfer in Three—Dimensional Parabolic Flows‘, Int. J. Heat Mass Transfer, Vol. 15, pp. 1787—1806.

Peric, M., 1987, "Efficient Semi-Implicit Solving Algorithm for Nine-Diagonal Coefficient Matrix", Num. Heat Transfer, Vol. 11, pp. 251—279.

Peric, M., Scheuerer, G., 1987, "CAST - A Finite Volume Method for the Computer Simulation of Two-Dimensional Flows", Lehrstuhl f. Strömungsmechanik, Universität Erlangen-Nürbnerg, LSTM-Bericht 165/T/87.

Stone, H.L., 1968, "Iterative Solution of Implicit Approximations of Multidimensional Partial Differential Equations‘, SIAM J. Num. Anal., Vol. 5, pp. 530-558.

Thole, C.-A., 1985, "Experiments with Multigrid Methods on the CalTech-Hypercube" Gesellschaft f. Mathem. und Datenverarbeitung, St. Augustin, GMD-Studien Nr. 103.

Vanka, S.P., 1986, "Block-Implicit Multigrid Solution of Navier-Stokes Equations in Primitive Variables", J. Comp. Physics, Vol. 65, pp. 138—158.

Wirth, N., 1985, "Programming in MODULA-2", 3. ed., Springer Verlag, Berlin.

 

SPIP | | Plan du site | Suivre la vie du site RSS 2.0
Habillage visuel © freelayouts sous Licence Creative Commons Attribution 2.5 License
{id_article}