FR2963692A1 - Dispositif informatique de calcul polyvalent - Google Patents

Abstract

Un dispositif informatique de calcul polyvalent du type comprend : - un calculateur-solveur (12), agencé pour recevoir une représentation matricielle de travail correspondant à un système d'équations, ainsi que des données de résidus, et pour fournir une solution du système d'équations à partir des données de résidus, - un adaptateur (10), agencé pour recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter, ainsi qu'une représentation matricielle filtrante pour ce système d'équations, et agencé pour calculer une représentation matricielle de travail correspondant à un système d'équations soluble par le calculateur-solveur , la représentation matricielle de travail étant contrainte à vérifier avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation matricielle de travail, L'adaptateur est agencé pour calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation numérique de représentation matricielle filtrante, tandis que le calculateur-solveur est agencé pour travailler sur cette matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci, tandis que ledit calcul itératif de l'adaptateur obéit à une règle de calcul où un bloc courant (i,j) de la matrice intermédiaire est défini par la différence entre le bloc (i,j) correspondant de la représentation matricielle initiale et une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliaire de la matrice d'approche.

Description

INRIAl21.FRD Dispositif informatique de calcul polyvalent

L'invention concerne la modélisation et la simulation des systèmes physiques complexes.

Dans de nombreux domaines de la physique moderne, les équations qui régissent un phénomène physique ne peuvent pas être résolues de manière théorique. C'est notamment le cas de tous les problèmes qui ont trait à la mécanique des fluides, par exemple dans la modélisation de l'exploitation d'un champ pétrolifère.

Dans ces situations, les équations différentielles sont résolues de manière numérique, c'est-à-dire par discrétisation des équations générales, en fonction des paramètres particuliers de la simulation.

15 Ces systèmes discrets sont résolus par l'utilisation de matrices de très grande taille, dans lesquelles les équations discrétisées forment les bases des systèmes. Mais ces matrices de base sont difficilement inversibles.

Pour résoudre ce problème, des méthodes itératives sont largement utilisées 20 aujourd'hui. Ces méthodes, par exemple celle dite « GMRES », sont basées sur des sous-espaces de Krylov. Dans le but d'accélérer la convergence des méthodes itératives, des "préconditionneurs" ont été créés. Ce sont des éléments qui calculent une matrice proche de la matrice de base, et dont l'inverse peut être appliqué de manière efficace à un vecteur arbitraire. 25 Les préconditionneurs sont intéressants, mais posent des problèmes avec certains vecteurs pour lesquels ils ne reproduisent pas fidèlement la matrice de base. Pour répondre à ce problème, des préconditionneurs "satisfaisant une propriété de filtrage" ont été développés, qui ont la particularité d'être fidèles à la matrice de base pour un 30 vecteur choisi particulier.

À ce jour, les méthodes pour produire des préconditionneurs satisfaisant une condition de filtrage nécessitent des matrices de base présentant une forme très spécifique, ce qui 10 limite fortement l'usage et l'utilité des préconditionneurs satisfaisant une propriété de filtrage.

L'invention vient améliorer la situation.

A cet effet, l'invention propose un dispositif informatique de calcul polyvalent du type comprenant : - un calculateur-solveur, agencé pour recevoir une représentation matricielle de travail et une représentation matricielle initiale correspondant à un système d'équations, ainsi que des données de résidus, et pour fournir une solution du système d'équations à partir des données de résidus, - un adaptateur, agencé pour recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter, ainsi qu'une représentation matricielle filtrante formant représentation matricielle filtrante pour ce système d'équations, et agencé pour calculer une représentation matricielle de travail correspondant à un système d'équations soluble par le calculateur-solveur, la représentation matricielle de travail étant contrainte à vérifier avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation matricielle de travail,

L'adaptateur est agencé pour calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation matricielle filtrante, tandis que le calculateur-solveur est agencé pour travailler sur cette matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci, tandis que ledit calcul itératif de l'adaptateur obéit à une règle de calcul où un bloc courant (i,j) de la matrice intermédiaire est défini par la différence entre le bloc (i,j) correspondant de la représentation matricielle initiale et une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliaire de la matrice d'approche.

Un tel dispositif est particulièrement avantageux, car il permet de mettre en oeuvre un préconditionneur satisfaisant une condition de filtrage, et ce quelque soit la matrice de base qui définit le système que l'on cherche à résoudre.

L'invention concerne également un procédé comprenant : a) recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter et une représentation matricielle filtrante, b) calculer une représentation matricielle de travail vérifiant avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation matricielle de travail, c) recevoir des données de résidus, et résoudre le système d'équations défini par la représentation matricielle initiale, à partir des données de résidus, de la représentation matricielle de travail et de la représentation matricielle initiale.

L'opération b) comprend calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation matricielle filtrante en répétant itérativement, pour chaque bloc d'indice courant (i,j) de la matrice intermédiaire : b1) calculer une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliaire de la matrice d'approche b12) calculer la différence entre le bloc (i,j) de la matrice de la représentation matricielle initiale et la somme issue de l'opération bl).

L'opération c) comprend travailler sur la matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci.

D'autres caractéristiques et avantages de l'invention apparaîtront mieux à la lecture de la description qui suit, tirée d'exemples donnés à titre illustratif et non limitatif, tirés des dessins sur lesquels : - la figure 1 représente une vue schématique d'un système de modélisation et de 10 simulation selon l'invention, - la figure 2 représente un diagramme de flux simplifié d'une opération de modélisation et de simulation au moyen du système de la figure 1, - la figure 3 représente un diagramme de flux simplifié d'une opération de la figure 2, - la figure 4 représente un exemple d'une fonction de calcul d'un préconditionneur selon 15 une opération de la figure 3, - la figure 5 représente un exemple de matrice obtenue après une première opération de réordonnancement selon une opération optionnelle de la figure 3, - la figure 6 représente un exemple de matrice obtenue après une deuxième opération de réordonnancement selon une opération optionnelle de la figure 3, et 20 - la figure 7 représente un exemple de modification de la fonction de la figure 4 pour mettre en place une parallélisation des calculs à partir d'une matrice réordonnée du type de celle représentée sur les figures 5 et 6.

Les dessins et la description ci-après contiennent, pour l'essentiel, des éléments de 25 caractère certain. Ils pourront donc non seulement servir à mieux faire comprendre la présente invention, mais aussi contribuer à sa définition, le cas échéant.

En outre, la description détaillée est augmentée de l'annexe A, qui donne la formulation de certaines formules mathématiques mises en oeuvre dans le cadre de l'invention. Cette 30 Annexe est mise à part dans un but de clarification, et pour faciliter les renvois. Elle est partie intégrante de la description, et pourra donc non seulement servir à mieux faire comprendre la présente invention, mais aussi contribuer à sa définition, le cas échéant.

La modélisation et la simulation des systèmes physiques sont devenues des enjeux de taille. Par exemple, dans l'exploitation d'un forage d'hydrocarbure, il y a une première phase pendant laquelle le pétrole sort naturellement. Puis, au fur et à mesure que la pression baisse, il devient nécessaire d'agir pour récupérer le pétrole.

Pour cela, il est possible par exemple d'utiliser un flux d'eau, qui est introduit dans le puits pour faire remonter la pression et faire rejaillir le pétrole. Mais ces opérations périlleuses nécessitent une connaissance poussée du puits et des réactions de celui-ci dans ces circonstances.

Les équations qui déterminent ce problème physique sont très complexes, et pour la plupart n'admettent que des solutions par discrétisation et méthode numérique de type différences finies ou volumes finis.

15 Les problèmes ainsi discrétisés peuvent alors être résumés à la formule (10) de l'Annexe A, dans laquelle A est la matrice de base qui définit le système d'équations discrétisé, x est le vecteur qui est recherché, et y est le vecteur résultat connu.

Ce type de problème est bien connu en algèbre, et il s'agit de trouver la matrice inverse 20 de A pour calculer x. Mais l'inversion des matrices est un problème complexe, qui monopolise des puissances de calcul qui croissent de manière exponentielle avec la taille de la matrice à inverser.

Pour cela, des méthodes itératives basées sur des sous-espaces de Krylov, comme 25 GMRES, sont largement utilisées aujourd'hui. Pour accélérer la convergence de ces méthodes, des "préconditionneurs" ont été proposés. Les préconditionneurs sont des matrices qui permettent d'approcher rapidement la matrice inverse de A. En utilisant le préconditionneur M, la méthode itérative résoud le système linéaire M-1 A x = M-1 b. 30 Dans ce mode de résolution, on calcule des opérations de type M-1 v et A v, où v est un vecteur, sans calculer de manière explicite l'inverse de M. 10

Comme cela a été expliqué plus haut, il existe une classe particulière de préconditionneurs, les préconditionneurs satisfaisant une propriété de filtrage.

Les préconditionneurs satisfaisant une propriété de filtrage présentent l'avantage supplémentaire de se comporter de manière identique à la matrice A pour un vecteur choisi, comme cela est explicité dans la formule (20) de l'annexe A, dans laquelle M est le préconditionneur, et t le vecteur choisi.

À ce jour, seules les matrices A qui sont tridiagonales par blocs sont utilisées pour 10 produire un préconditionneur satisfaisant une propriété de filtrage. Cela restreint considérablement leur champ d'application.

En outre, les méthodes de calcul de ces préconditionneurs sont pour la plupart séquentielles, ce qui les rend rapidement prohibitives en termes de coût de calcul, et 15 donc peu utilisables en pratiques. En effet, seules les matrices issues d'un maillage structuré peuvent être traitées en parallèle, ce qui restreint considérablement leur champ d'application.

La figure 1 représente un dispositif informatique de calcul polyvalent 2 selon 20 l'invention. Le dispositif 2 comprend un ensemble de capteurs 4, un numériseur 6, un discrétiseur 8, un adaptateur 10, un calculateur-solveur 12, et un pilote 14 qui les commande.

Dans l'exemple décrit ici, l'ensemble de capteurs 4 servent à obtenir les données qui 25 contraignent le système physique à modéliser, et le numériseur 6 sert à transformer ces données analogiques pour les injecter dans des équations théoriques.

Ces éléments sont pour ainsi dire indifférents au problème résolu par l'invention : ils servent à définir le cadre pour son application pratique. Aussi, leur réalisation pourra 30 être très variée.

Le discrétiseur 8 est appelé par le pilote 14 pour discrétiser les équations théoriques particularisées avec les données réelles, et pour en tirer un système d'équations linéaires.

Ce système présente en général une taille très importante, et ses lignes forment la matrice A. Là encore, cet élément peut être réalisé de nombreuses manières différentes.

Enfin, l'adaptateur 10 et le calculateur 12 sont appelés par le pilote 14 pour calculer le préconditionneur et pour tirer une solution correspondant à une situation particulière que l'on cherche à modéliser. Dans ce cas, les données de "second membre" de l'équation faisant intervenir la matrice A sont également appelées données de résidus, en référence aux méthodes de Newton.

Le pilote 14 peut appeler l'adaptateur 10 et le calculateur 12 pour faire évoluer cette solution particulière par pas de temps successifs, et pour donner ainsi une simulation de l'évolution du système physique modélisé.

Selon une première variante, l'adaptateur 10 est appelé une unique fois par le pilote 14 pour calculer un préconditionneur qui est utilisé pour toute la durée de la simulation, et le résultat calculé par le calculateur 12 pour un pas de temps donné est utilisé comme entrée au pas de temps suivant.

Selon une deuxième variante, l'adaptateur 10 peut être sélectivement appelé par le pilote 14, en fonction de l'évolution de la simulation, notamment si celle-ci tend à modifier le système à l'origine de la matrice A.

Ici encore, les techniques de simulation sont variées.

La figure 2 montre un diagramme de flux simplifié montrant les opérations résumées plus haut: - dans une opération 200, l'ensemble de capteurs 4 est appelé pour mesurer tous les paramètres nécessaires à la simulation, - dans une opération 220, le numériseur 6 et le discrétiseur 8 sont appelés pour modéliser le système de manière numérique, avec les mesures tirées de l'opération 200, - dans une opération 240, l'adaptateur 10 et le calculateur 12 sont appelés pour effectuer la simulation en tant que telle.

La figure 3 montre un diagramme de flux simplifié de l'opération 240.

Dans une opération 300, le pilote 14 transmet la matrice A tirée de l'opération 220 à l'adaptateur 10.

Dans une opération 320 optionnelle, l'adaptateur 10 réordonne les éléments de la matrice A pour permettre un traitement ultérieur en parallèle. Cette opération peut être réalisée de plusieurs manières, par exemple par une dissection emboîtée ou par une partition en plusieurs domaines indépendants qui peuvent également se recouvrir et avoir une sous-partition récursive.

Un tel recouvrement augmente légèrement les coûts de calcul, mais offre de meilleurs taux de convergence et une robustesse supérieure, comme cela est réalisé dans la méthode de Schwarz. L'adaptateur 10 permet ainsi d'obtenir une matrice B réordonnée qui comprend des blocs nuls. Ensuite, dans une opération 340, l'adaptateur 10 traite la matrice A, réordonnée ou pas, pour en tirer une représentation d'une préconditionneur M satisfaisant une propriété de filtrage. Enfin, dans une opération 360, le pilote 14 appelle le calculateur 12 avec la représentation du préconditionneur M pour réaliser la simulation. 20 Le dispositif formé par l'adaptateur 10, le calculateur 12 et le pilote 14, permet donc : - de calculer une représentation d'un préconditionneur vérifiant une propriété de filtrage pour une quelconque matrice A en entrée, et - paralléliser de manière massive les calculs liés au préconditionneur lorsque 25 l'adaptateur 10 est appelé pour réordonner la matrice A.

La figure 4 représente un diagramme de flux de calcul du préconditionneur selon l'invention.

30 Ce calcul repose sur la décomposition du préconditionneur M sous la forme de la formule (30) de l'annexe A. Cette décomposition en matrice LDU est connue en principe, mais le calcul des éléments est différent. Les formules (40), (50) et (60) de l'annexe A donnent la composition de ces éléments respectifs.15 La décomposition du préconditionneur sous la forme LDU est très avantageuse, car elle permet de résoudre le système réel sans avoir à inverser la matrice M. Plus précisément, la technique contient de nombreux algorithmes qui permettent une résolution simplifiée d'une équation matricielle lorsque la décomposition LDU est utilisée. Pour cette raison, le calcul du préconditionneur M en tant que tel n'est jamais réalisé, et seules ses composantes LDU sont calculées et stockées. Ensuite, le calculateursolveur 12 les appelle sélectivement pour résoudre le système.

10 Il serait néanmoins possible de calculer le préconditionneur M, en appliquant la formule (30).

Pour calculer les éléments de la décomposition LDU, les Demandeurs ont découvert une formulation basée sur le calcul d'une matrice C qui correspond à la somme des 15 matrices L, D et U (formule (70) de l'Annexe A).

Du fait de la forme respective des matrices L, D et U, il apparaît que chaque élément de la matrice C correspond à un unique élément de chacune des matrices L, D, U.

20 Ainsi : - Di, = C;,, et Du = 0 pour i différent de j, - L1 = Cu pour i>1 et i strictement supérieur à j, et L1 = 0 pour i inférieur ou égal à j, - Uu = Cu pour j> l et j strictement supérieur à i, et Uu = 0 pour j inférieur ou égal à i.

25 Les Demandeurs ont découvert que la matrice C peut être établie conformément à la formule (80) de l'annexe A, dans laquelle le terme Fkj satisfait la formule (90) de l'annexe A.

Dans la formule (80), la première ligne ne représente qu'une initialisation. 30 Conceptuellement, cette première ligne est équivalente à la deuxième ligne. En effet, pour i=1 ou j=1, alors min(i,j)-1 vaut zéro, ce qui signifie que la somme de cette deuxième ligne ne comprend aucun terme, et donne un résultat identique à la première ligne.5 La formule (90) exprime le fait que la matrice F qui regroupe les termes Fkk est une matrice qui approche le bloc diagonal Dkk pour la condition portant sur l'indice j. Ainsi, le fait que Fkk satisfait la formule (90) peut être vu comme une condition d'équivalence avec l'inverse du bloc diagonal Dkk pour la condition portant sur l'indice j.

Lorsque le vecteur Ukktj n'a aucune composante nulle, le calcul de FkJ est aisé. Cependant, pour tenir compte des cas où ce vecteur présente des composantes nulles, il est possible de modifier les approches de la littérature conformément à la formule (100) de l'annexe A.

Les Demandeurs ont également découvert un autre calcul pour la matrice C, dans laquelle, dans la formule (80), le terme Fkk est remplacé par le terme GkJ, qui est défini avec la formule (110) de l'annexe A. Pour calculer le terme Gki, on calcule d'abord le terme Fkj correspondant grâce à la formule (100), puis on applique la formule (110). 15 Il est aussi possible de définir le terme Fki satisfaisant la condition de filtrage (90) par la formule (120) combinée à la formule (140) de l'Annexe A, qui est issue des méthodes de déflation de l'algèbre linéaire. Dans le cas où les blocs Dkk sont symétriques, on peut utiliser la formule (130) de l'Annexe A, qui est une version simplifiée de la 20 formule (120).

Dans l'état actuel des recherches des Demandeurs, la formule (100) est préférée à la formule (110), pour des raisons de stabilité.

25 La figure 4 représente un exemple de fonction permettant de calculer les matrices L, D et U. Pour cela, chaque terme de la matrice est calculé, et affecté à la matrice L, D et U comme il convient.

En variante, les termes ne sont pas affectés à chaque matrice L, D, U, mais à la seule 30 matrice C, qui est directement utilisée par le calculateur-solveur 12. En fait, on a vu plus haut que la matrice C est équivalente aux matrices L, D, et U, et ces deux variantes représentent seulement des manières différentes d'exprimer le préconditionneur M. 10 Au vu de ce qui précède, le premier élément de la matrice D peut être calculé directement, puisqu'il correspond au premier terme de la diagonale de la matrice A (ou B si l'opération 320 est exécutée).

De même tous les termes non nuls de respectivement la première colonne de la matrice L et la première ligne de la matrice U sont initialisés avec le terme correspondant de la matrice A. Un indice i est initialisé à 1. Cela est réalisé dans une opération 400.

Par premier élément ou premier terme, on entend bloc. En effet, la matrice A (ou B si l'opération 320 est exécutée) est prédécoupée en blocs rectangulaires dont les côtés sont des paramètres qui peuvent être choisis librement. Seuls les blocs diagonaux de A doivent être carrés. À ce jour, les Demandeurs utilisent des valeurs de côtés telles que le produit de ces valeurs est égal à la taille de la mémoire tampon, c'est-à-dire qu'un bloc donné de la matrice A peut être stocké dans le tampon mémoire.

En variante, les Demandeurs utilisent également de valeurs de côtés telles que leur 20 produit est inférieur à la taille de la mémoire tampon. Si l'opération 320 est exécutée, la taille des blocs de la matrice B est déterminée par cette opération.

Ensuite, on exécute une boucle dite globale qui va calculer tous les autres termes de la matrice C, et donc tous les termes des matrices L, D, et U, ce qui permet de définir le 25 préconditionneur M.

La boucle globale consiste, à chaque itération, à calculer d'abord le terme diagonal, puis à calculer au moyen d'une boucle locale les autres termes, par indice croissant de ligne et de colonne. D'abord, l'indice i de la boucle globale est incrémenté dans une opération 402, puis dans une opération 404, le terme diagonal D;; est calculé conformément à la formule (80). 30 Le terme D. correspond comme on l'a vu précédemment au terme C,;. L'opération 404 inclut le calcul de la matrice F tel que la formule (90) soit satisfaite, conformément à la formule (100). Si le deuxième mode de calcul est retenu, le calcul de G est également réalisé, conformément à la formule (110).

Ensuite, l'indice j de la boucle locale est initialisé à 2 dans une opération 406. Cela est suivi d'une opération 408 de fin de boucle locale laquelle il est testé si j est égal à i. Quand i vaut 2, cela permet de passer directement à l'itération de boucle globale suivante, comme C12 et C21 sont connus.

La boucle locale est alors exécutée, avec les calculs de L1 qui correspond à dans une opération 410, et de qui correspond à c,, dans une opération 412.

On notera que les opérations 410 et 412 peuvent être réalisées en parallèle, ce qui est 15 avantageux. En effet, de par la formule (80), le calcul des termes C,i et Ci; est indépendant.

Ensuite, l'indice de la boucle locale j est incrémenté dans une opération 414, et la boucle locale reprend avec le test de l'opération 408. Lorsque la boucle locale est terminée, 20 c'est-à-dire lorsque tous les termes de la ligne de L et tous les termes de la colonne de U ont été calculés, un test vérifie dans une opération 416 si i est égal au nombre de blocs N de la matrice M.

Si c'est le cas, alors la boucle globale est terminée, et l'opération se termine dans une 25 opération 418 l'appel du calculateur-solveur (12). Sinon, la boucle globale reprend avec l'incrémentation de l'indice i de la boucle globale dans l'opération 402.

Comme mentionné plus haut, la fonction de la figure 4 permet d'obtenir une décomposition de M sous la forme LDU, qui sert de base aux méthodes de résolutions 30 connues en algèbre. Ainsi, le préconditionneur M n'est pas explicitement calculé. Cependant, il serait possible de calculer explicitement le préconditionneur, en appliquant la formule (30). 10 Les Demandeurs ont développé un dispositif mettant en oeuvre un préconditionneur et une méthode de calcul d'une représentation de ce préconditionneur qui sont particulièrement adaptés à la parallélisation des calculs.

Pour mieux comprendre cela, il convient d'expliquer plus en détail un exemple de réalisation de l'opération 320. Cet exemple sera ici basé sur le cas d'une dissection emboîtée à deux niveaux.

Dans ce type d'opération, la matrice A est modifiée pour lui donner la forme d'une matrice B en forme de "flèche". Pour cela, la matrice A est "réordonnée" une première fois pour lui donner la forme de la matrice représentée sur la figure 5, puis les sous-matrices de la matrice de la figure 5 sont elles-mêmes réordonnées de la même manière.

La matrice de la figure 5 comporte trois blocs diagonaux B1, B2 et B3, deux blocs B4 et B5 respectivement le long des bords bas et droit de la matrice B, et est nulle ailleurs. Ce réordonnancement de la matrice A est possible du fait de sa faible densité. Le même réordonnancement est réalisé sur le vecteur t.

Une fois que la matrice de figure 5 a été calculée, il est possible de réappliquer ce même réordonnancement aux blocs B1 et B2, ce qui aboutit à la matrice B de la figure 6. On remarquera sur cette figure que le bloc B3 de la figure 5 est le bloc B77 de la figure 6, et que les blocs B4 et B5 de la figure 5 correspondent respectivement aux blocs B71 à B76 et aux blocs B17 à B67 de la figure 6.

Si on analyse la formule (80) de l'Annexe A, il apparaît par la première ligne de cette formule que les blocs CI l à C17 et C21 à C71 sont connus directement à partir de B. Il apparaît également que, dans l'application de la deuxième ligne de cette formule, de nombreux blocs sont nuls ou connus, ce qui rend possible de calculer certains blocs C,i de manière indépendante les uns des autres, et donc que leur calcul peut être réalisé en parallèle.

Ainsi, B11, B22, B44 et B55. Si l'on calcule ces blocs en parallèle, la même situation se répète, et de nouveaux blocs B33 et B66 peuvent à leur tour être calculés en parallèle. Et ainsi de suite.

D'une manière générale, les Demandeurs ont donc découvert que l'opération 320 peut être utilisée pour produire une matrice B équivalente à la matrice A, et qui comporte des domaines séparés. À partir de ces domaines, il est possible de créer un graphe de dépendance de tâches, dans lequel pour un niveau donné, tous les noeuds représentent des blocs pouvant être calculés en parallèle pour l'application de la formule (80) de l'Annexe A.

10 Une fois tous les blocs liés aux noeuds d'un niveau donné calculés, il devient possible de calculer les blocs liés aux noeuds du prochain niveau dans le graphe, là encore en parallèle. Pour calculer et ordonnancer le graphe de dépendance de tâches, plusieurs techniques peuvent être utilisées, comme cela est connu dans la théorie des graphes.

15 Le graphe de dépendances de tâches est calculé en parcourant la structure de la matrice B. Les noeuds de ce graphe représentent des tâches, c'est à dire des calculs de blocs Cu de la formule (80) de l'Annexe A.

Les dépendances dans ce graphe représentent l'ordre des calculs imposé par la 20 formule (80) de l'Annexe A. Ce graphe peut être ensuite ordonnancé en utilisant un ordonnancement statique ou un ordonnancement dynamique des tâches sur les processeurs.

Par exemple, un ordonnancement dynamique affecte au cours de l'exécution parallèle 25 les tâches prêtes à être exécutées sur les processeurs disponibles. Un ordonnancement statique établit dans une première phase l'ordre des tâches à exécuter sur les différentes unités de calcul parallèles dans le but de minimiser le temps de calcul parallèle.

Dans une deuxième phase leur exécution a lieu. Une distribution statique des données 30 peut utilisée sur les processeurs de type sous-arbre vers sous-cube ou de type bidimensionnel.

C'est cette observation qui a amené les modifications du diagramme de la figure 4, qui sont présentées avec la figure 7.5 Dans cette fonction la première opération 700 correspond à l'opération 400 de la figure 4. Une différence avec l'opération 400 est que l'opération 700 comprend le calcul du graphe de dépendance de tâches décrit plus haut. Ensuite l'indice i, qui correspond au niveau du graphe de dépendance de tâches qui est actuellement parcouru, est incrémenté dans une opération 710.

Une fois l'indice i incrémenté, le pilote 14 appelle l'adaptateur 10 dans une 10 opération 720 pour récupérer tous les noeuds du niveau i du graphe de dépendance de tâches, au moyen d'une fonction Dep_Gr().

Le résultat de cette fonction est stocké dans une liste List, qui est une variable locale qui contient à chaque itération la liste des couples (k,l) qui identifient les blocs 15 indépendants de même niveau du graphe de dépendance de tâches.

Ensuite, l'adaptateur calcule les blocs Ckl pour tous les couples de la list List dans une opération 730. Ce calcul est exécuté en parallèle, comme tous les blocs sont indépendants par rapport à la formule (80) de l'Annexe A, et sont répartis sur tous les 20 processeurs et les coeurs de processeurs disponibles.

On remarquera ici que le calcul des blocs par l'opération 730 est différent de celui des opérations 404, 410 et 412. En effet, si la formule utilisée pour le calcul est la même, les indices des blocs sont totalement indépendants. Là où la fonction de la figure 4 opère en calculant d'abord le terme diagonal, puis les termes de la ligne et la colonne correspondante, ici c'est le graphe de dépendance de tâches qui détermine quels sont les blocs calculés.

30 On notera que dans l'exemple présent, la nature de la matrice B crée une certaine symétrie des indices de chaque niveau dans le graphe de dépendance de tâches. Cependant, l'application d'une autre méthode que la dissection emboîtée pourra limiter cette symétrie, et les blocs pourront être calculés dans un ordre apparemment arbitraire. 25 Enfin dans une opération 740, l'adaptateur 10 vérifie si l'indice i est inférieur à N, le nombre de niveaux du graphe. Si c'est le cas, alors la fonction reprend en 710 avec l'incrémentation de i pour le niveau suivant du graphe de dépendance de tâches. Sinon, la fonction se termine dans l'opération 750 avec l'appel du calculateur-solveur 12, comme avec l'opération 418.

Dans ce qui précède, la condition de filtrage est exprimée par la formule (20) de l'Annexe A. Cette formule est une expression mathématique du fait que la matrice initiale A et le préconditionneur M vérifient une condition de stabilité qui est basée sur la comparaison de leur produit avec un vecteur.

Cependant, la condition de stabilité ne doit pas être limitée à la seule formule (20). Ainsi, les Demandeurs ont également utilisé avec succès la formule (150) de l'Annexe A.

Comme la formule (150) est presque la transposée de la formule (20), l'utilisation de la formule (150) comme condition de stabilité ne change rien au mode de fonctionnement de l'invention. En conséquence, les formules (80) et (90) n'ont qu'à être légèrement modifiées, comme présenté avec les formules (160) et (170). Les formules (100) à (140) pourront être adaptées de manière similaire.

En outre, tous les exemples qui précèdent ont été réalisés pour une condition de stabilité utilisant un vecteur t.

Cependant, lorsqu'un système physique est modélisé, de nombreuses grandeurs sont utilisées. Les expériences des Demandeurs montrent qu'il est avantageux d'utiliser une condition de stabilité utilisant une matrice dont chaque colonne concerne une grandeur physique.

Ainsi, si deux grandeurs physiques caractérisent une mise en équation donnée, il est avantageux d'utiliser comme élément de filtrage t une matrice ayant deux colonnes. Dans la pratique, cela ne change pas la philosophie de l'invention, et les calculs présentés précédemment s'en trouvent peu ou pas modifiés.

En effet, dans ce type de situation, les équations représentées dans la matrice A vont être associées par "mini-blocs" carrés dont le côté est égal au nombre de colonne de la matrice t. Donc, dans le cas décrit au paragraphe précédent, chaque mini-bloc serait un bloc carré de deux fois deux termes de la matrice A initiale.

La raison pour laquelle ces mini-blocs sont mentionnés est qu'ils ne doivent pas être séparés lors de l'opération optionnelle 320, et lorsque la matrice A est découpée en blocs. Un mini-bloc donné doit toujours être contenu dans un unique bloc de la matrice A ou de la matrice B.

Le seul élément qui change légèrement est le calcul de Fki. En effet, la formule (100) de l'Annexe A est adaptée à une matrice t à une seule colonne, c'est-à-dire un vecteur, et les mini-blocs sont donc de taille un fois un, c'est-à-dire des scalaires.

15 Les Demandeurs ont donc généralisé la formule (100) sous la forme de la formule (180), dans laquelle Diag() désigne une fonction qui crée une matrice diagonale dont les éléments sont désignés en arguments de cette fonction, et dans laquelle l'opération "/." désigne la division terme à terme des matrices. Ainsi Al/.A2 est une matrice A3 dont chaque terme A3(i,j) est égal au quotient de A1(i,j) par A2(i,j). 20 Une autre manière de voir cette modification est de remarquer que la formule (100) peut être vue comme un cas particulier de la formule (180), dans laquelle t a une seule colonne.

25 Dans ce qui précède, l'adaptateur 10, le calculateur 12 et le pilote 14 peuvent être réalisés de plusieurs manières.

Tout d'abord, le pilote 14 peut être réalisé de manière intégrée avec l'adaptateur 10 et le calculateur 12, c'est-à-dire que ceux-ci sont agencés pour savoir interagir, au lieu d'être 30 des éléments séparés commandés qui s'ignorent.

De plus, la présentation des éléments du système 2 est principalement fonctionnelle. Ainsi, ces éléments peuvent être séparés physiquement et reliés par des liens de communication, ou mis en oeuvre de manière éloignée dans le temps, ou encore mis en 10 oeuvre sur un même équipement avec le pilote 14 défini par les liaisons intrinsèques entre ces éléments et une interface utilisateur.

En outre, le discrétiseur 8, l'adaptateur 10, le calculateur 12 et le pilote 14 peuvent être mis en oeuvre sous la forme d'éléments analogiques, comme des circuits intégrés ou des cartes filles, ou sous la forme d'éléments numériques, c'est-à-dire sous la forme de programmes mis en oeuvre par un ordinateur, éventuellement de manière éloignée et/ou répartie.

On notera également que, dans ce qui précède, il est souvent indifféremment fait référence à des matrices ou à leur représentation. Il va de soi qu'un ordinateur ne sait pas ce qu'est une matrice, et que c'est bien la représentation numérique de ces matrices, c'est-à-dire les données qui définissent ces matrices qui sont visées. Par matrice ou par représentation matricielle, on entend donc toute structure de données numériques qui permet de traiter la matrice dans le cadre de l'invention.

On notera enfin la visée particulièrement pratique du dispositif de l'invention, qui permet la simulation et la résolution de nombreux problèmes physiques qui n'étaient pas solubles précédemment, par exemple dans l'industrie pétrolière. 19 (10) (20) D11 U12 U1N

DN-1N-1 UN-1N DNN r LN L21 LN1 Die1 D22 LNN-1 DNN (30) 10 L LNN-1 (40) D= D22 (50) U12 U1N UN-1N 0 (60) U= 15 C=L+D+U (70) iJ ZJk=1,Lik#o,UkjO ik kikj Ai./ sii=1ouj=1 Cid = 1A.. _ min(i,f)-1 L F U Si i > 1 ou " > 1 (80) 20 FkiUkiti = Dkk Uk~ t~ (90) Fkj (r, s) = (100) DkkUkjtj(r)/Ukjtj(r) si r = s et Ukjtj(r) # 0 Dkk Ukjtj(r)/Ukjtj(s) si Ukjtj(r) = 0 avec s tel que Is - ri = mink:Ukjtj(k)OIk - rI 0 sinon Gk1 = 2Fkj - FkJDkkFkj (110) Q = Z E-1(Dkk Z)T, E = ((DkkZ)T DkkZ)-1, Z = Ukj tj Q = Z E' ZT,E = (ZTDkkZ)-1,Z = Ukjtj (130) Fkj=P+Q,P=I-QA (140) tTA = tTM (150) Aijsij=fou]=1 Ci j = At) - Ek=1L=,L j)-ikO,1UkjO LtkFtkUkj ~ si i > 1 ou " > 1 trLikFik = ti LikDkk (170) Fkj (r, s) = (180)

Diag (Dkk Uk j tj (r)/. Uk j tj (r)) si r = s et Uk j tj (r) n'a aucune composante nulle Diag(DkkUkjti(r)/. Ukjtj(s)) si Ukjtj(r) = 0 avec s tel que is - ri = min lk - ri k:Ukjt/(k)*0 (120) (160) 0 sinon20

Claims (11)

1. REVENDICATIONS1. Dispositif informatique de calcul polyvalent pour la modélisation et la simulation de systèmes physiques, du type comprenant : - un calculateur-solveur (12), agencé pour recevoir une représentation matricielle de travail (M) et une représentation matricielle initiale (A) correspondant à un système d'équations représentant un système physique, ainsi que des données de résidus, et pour fournir une solution du système d'équations à partir des données de résidus, - un adaptateur (10), agencé pour recevoir une représentation matricielle initiale (A) correspondant à un système d'équations à traiter, ainsi qu'une représentation matricielle filtrante (t) pour ce système d'équations, et agencé pour calculer une représentation matricielle de travail (M) correspondant à un système d'équations soluble par le calculateur-solveur, la représentation matricielle de travail (M) étant contrainte à vérifier avec la représentation matricielle initiale (A) une condition de stabilité ((20), (150)) comprenant une comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante (t) ou sa transposée, et comportant respectivement la représentation matricielle initiale (A), et la représentation matricielle de travail (M), caractérisé en ce que l'adaptateur est agencé pour calculer itérativement par blocs une matrice intermédiaire (C) à partir de la représentation matricielle initiale (A) et de ladite représentation numérique de représentation matricielle filtrante (t), tandis que le calculateur-solveur est agencé pour travailler sur cette matrice intermédiaire, par blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale (A), sans inversion complète de celle-ci, tandis que ledit calcul itératif de l'adaptateur obéit à une règle de calcul ((80) ; (160)) où un bloc courant (i,j) de la matrice intermédiaire (C) est défini par la différence entre le bloc (i,j) correspondant de la représentation matricielle initiale (A) et une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire (C), et un bloc auxiliaire d'une matrice d'approche (F) qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence ((90), (170)), comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante (t) ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire (C), et ledit bloc auxiliaire de la matrice d'approche (F). 21
2. 2. Dispositif selon la revendication 1, dans lequel la condition de stabilité (20) comprend une comparaison de deux produits matriciels, tous deux entre respectivement la représentation matricielle initiale (A), et la représentation matricielle de travail (M), et ladite représentation matricielle filtrante (t), et dans lequel la somme de blocs est réalisée, pour un indice k non nul et strictement inférieur au minimum entre i et j, chaque bloc de cette somme étant définie par le produit matriciel de la forme PQR, * où P est le bloc (i,k) de la matrice intermédiaire (C), * où Q est le bloc auxiliaire (k,j) de la matrice d'approche (F) d'indice, la condition d'équivalence (90) comprenant la comparaison de deux produits matriciels de respectivement ledit bloc auxiliaire (k,j) de la matrice d'approche (F) et l'inverse du bloc diagonal (k,k) déjà calculé de la matrice intermédiaire (C), avec le bloc (k,j) déjà calculé de la représentation matricielle intermédiaire (C) et avec le bloc j de ladite représentation matricielle filtrante (t), et * où R est le bloc (k,j) de la matrice intermédiaire (C).
3. 3. Dispositif selon la revendication 1, dans lequel la condition de stabilité (150) comprend une comparaison de deux produits matriciels, tous deux entre la transposée de la représentation matricielle filtrante (t) et respectivement la représentation matricielle initiale (A), et la représentation matricielle de travail (M), et dans lequel la somme de blocs est réalisée, pour un indice k non nul et strictement inférieur au minimum entre i et j, chaque bloc de cette somme étant définie par le produit matriciel de la forme PQR, * où P est le bloc (i,k) de la matrice intermédiaire (C), * où Q est le bloc auxiliaire (i,k) de la matrice d'approche (F), la condition d'équivalence (170) comprenant la comparaison de deux produits matriciels du bloc i de ladite représentation matricielle filtrante (t) avec le bloc (i,k) déjà calculé de la représentation matricielle intermédiaire (C), et avec respectivement ledit bloc auxiliaire (i,k) de la matrice d'approche (F) et l'inverse du bloc diagonal (k,k) déjà calculé de la matrice intermédiaire (C), et * où R est le bloc (kj) de la matrice intermédiaire (C).
4. 4. Dispositif selon l'une des revendications précédentes, dans lequel le bloc auxiliaire de la matrice d'approche (F) est calculé à partir d'une division terme à terme faisant intervenir le bloc déjà calculé de la matrice intermédiaire (C) d'indice (k,k), et soit lebloc j de ladite représentation matricielle filtrante (t) et le bloc (k,j) déjà calculé de la matrice intermédiaire (C), soit le bloc i ladite représentation matricielle filtrante (t) et le bloc (i,k) déjà calculé de la matrice intermédiaire (C).
5. 5. Dispositif selon l'une des revendications 1 à 3, dans lequel la matrice d'approche est calculée par une méthode de déflation, dans laquelle un premier terme (Z) fait intervenir soit le bloc j de ladite représentation matricielle filtrante (t) et le bloc (k,j) déjà calculé de la matrice intermédiaire (C), soit le bloc i de ladite représentation matricielle filtrante (t) et le bloc (i,k) déjà calculé de la matrice intermédiaire (C), dans laquelle un second terme (Q) fait intervenir le bloc diagonal (k,k) déjà calculé de la matrice intermédiaire (C) et le premier terme, et un troisième terme (P) fait intervenir les premier et seconds termes ainsi que la matrice de la représentation matricielle initiale (A).
6. 6. Dispositif selon l'une des revendications précédentes, dans lequel la représentation matricielle filtrante (t) est un vecteur colonne.
7. 7. Dispositif selon l'une des revendications précédentes, dans lequel l'adaptateur (10) est agencé pour réordonner la représentation matricielle initiale (A) pour produire une représentation matricielle modifiée (B) selon une règle d'ordonnancement agencée pour associer des blocs de la matrice de la représentation matricielle initiale (A) en fonction d'une condition de dépendance ((80),(160)), et pour calculer la représentation matricielle de travail (M) en remplaçant représentation matricielle initiale (A) par la représentation matricielle modifiée (B).
8. 8. Dispositif selon la revendication 7, dans lequel l'adaptateur (10) est agencé pour calculer un graphe à partir de la représentation matricielle modifiée (B) dans lequel pour un niveau donné les calculs des blocs de la matrice intermédiaire (C) peuvent être réalisés de manière indépendante, et pour calculer les blocs d'un même niveau du graphe en parallèle.
9. 9. Dispositif selon l'une des revendications précédentes, comprenant en outre un ensemble de capteurs 4, un numériseur 6, un discrétiseur 8 et un pilote 14, dans lequel le pilote 14 est agencé pour appeler le discrétiseur 8 avec des données tirées du numériseur 6 qui opère sur des données tirées de l'ensemblede capteur 4, pour produire la représentation matricielle initiale (A) et les données de résidus, et agencé pour commander l'adaptateur (10) et le calculateur-solveur (12) en conséquence.
10. 10. Dispositif selon l'une des revendications précédentes, dans lequel le système d'équations est représentatif d'un système physique complexe du monde réel, tel qu'un champ pétrolifère.
11. 11. Procédé de calcul polyvalent pour la modélisation et la simulation de systèmes physiques, du type comprenant : a) recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter représentant un système physique et une représentation matricielle filtrante, b) calculer une représentation matricielle de travail vérifiant avec la représentation matricielle initiale une condition de stabilité comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée, et comportant respectivement la représentation matricielle initiale, et la représentation matricielle de travail, c) recevoir des données de résidus, et résoudre le système d'équations défini par la représentation matricielle initiale, à partir des données de résidus, de la représentation matricielle de travail, et de la représentation matricielle initiale, caractérisé en ce que l'opération b) comprend calculer itérativement par blocs une matrice intermédiaire à partir de la représentation matricielle initiale et de ladite représentation numérique de représentation matricielle filtrante en répétant itérativement, pour chaque bloc courant (i,j) de la matrice intermédiaire, : b1) calculer une somme de blocs chacun défini par un produit faisant intervenir deux blocs déjà calculés de la matrice intermédiaire, et un bloc auxiliaire d'une matrice d'approche qui est contraint à vérifier avec un bloc diagonal déjà calculé de la matrice intermédiaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante ou sa transposée et un bloc déjà calculé de la matrice intermédiaire, et comportant respectivement l'inverse dudit bloc diagonal déjà calculé de la matrice intermédiaire, et ledit bloc auxiliairede la matrice d'approche, b2) calculer la différence entre le bloc (i,j) de la matrice de la représentation matricielle initiale et la somme issue de l'opération b1). et ce que l'opération c) comprend travailler sur la matrice intermédiaire, par 5 blocs, de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci.