Brevet FR2965080A1 — publié le 20120323

INRIAl22.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. Des méthodes dites directes comme l'élimination de Gauss permettent de résoudre ces 20 systèmes linéaires sans inversion. Mais ces méthodes directes sont très gourmandes en temps de calcul et en utilisation mémoire, ce qui les rend inutilisables en pratique pour des matrices de très grande taille. Pour résoudre ce problème, des méthodes itératives sont largement utilisées 25 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. 30 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 10 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 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 limite fortement l'usage et l'utilité des préconditionneurs satisfaisant une propriété de filtrage. L'invention vient améliorer la situation. À 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 la représentation matricielle initiale, ainsi qu'une représentation matricielle filtrante pour le système d'équations correspondant, 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 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é d'une part pour renuméroter la représentation matricielle initiale et la représentation matricielle filtrante pour produire une représentation matricielle modifiée et une représentation matricielle filtrante modifiée selon une règle d'ordonnancement agencée pour associer des blocs de la matrice de la représentation matricielle initiale en fonction d'une condition de dépendance, et d'autre part pour calculer récursivement la représentation matricielle de travail à partir de la représentation matricielle modifiée et de ladite représentation matricielle filtrante modifiée, tandis que le calculateur-solveur est agencé pour travailler récursivement sur la représentation matricielle de travail de façon à fournir une solution du système d'équations de la représentation matricielle initiale, sans inversion complète de celle-ci. Ledit calcul récursif de l'adaptateur comprend le calcul de la représentation matricielle de travail sous la forme d'un produit matriciel PQR : * où P est la somme entre d'une part une matrice diagonale calculée à partir d'une matrice auxiliaire et d'une matrice d'approche, et d'autre part une matrice triangulaire inférieure par blocs dont seuls les termes de la dernière ligne de blocs et non diagonaux sont non nuls et sont égaux aux termes de même indice dans la représentation matricielle modifiée, * où Q est l'inverse de la matrice diagonale, * où R est la somme entre la matrice diagonale et une matrice triangulaire supérieure par blocs dont seuls les termes de la dernière colonne de blocs et non diagonaux sont non nuls, et sont égaux aux termes de même indice dans la représentation matricielle modifiée, * la matrice auxiliaire étant une matrice diagonale par blocs, dont chaque bloc d'indice i est : - défini égal au bloc diagonal d'indice i de la représentation matricielle modifiée lorsque ce bloc ne vérifie pas la condition de récursivité, et - sinon calculé par un appel récursif à l'adaptateur avec le bloc diagonal d'indice i de la représentation matricielle modifiée comme représentation matricielle modifiée et avec un sous-ensemble de la représentation matricielle filtrante modifiée tiré de l'indice i comme représentation matricielle filtrante modifiée, * la matrice d'approche est une matrice diagonale par blocs dont le dernier bloc est nul, et dont chaque bloc non nul d'indice i est contraint à vérifier avec le bloc diagonal d'indice i de la matrice auxiliaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux respectivement ladite représentation matricielle filtrante modifiée ou sa transposée et respectivement un bloc de la matrice triangulaire supérieure par blocs ou un bloc de la matrice triangulaire inférieure par blocs, et comportant respectivement l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire, et ledit bloc d'indice i de la matrice d'approche, * les blocs diagonaux de la matrice diagonale étant égaux aux blocs de même indice de la matrice auxiliaire, sauf pour le dernier, qui est défini comme la différence entre le dernier bloc de la matrice auxiliaire et la somme pour un indice k non nul et inférieur au nombre de blocs diagonaux de la représentation matricielle modifiée de produits matriciels sous la forme WXY où W est le bloc non nul de la k-ème colonne de la matrice triangulaire inférieure par blocs, X est le k-ème bloc de la matrice d'approche, et Y est le bloc non nul de la k-ème ligne de la matrice triangulaire supérieure par blocs. 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 : b 1) renuméroter la représentation matricielle initiale et la représentation matricielle filtrante pour produire une représentation matricielle modifiée et une représentation matricielle filtrante modifiée, et calculer récursivement la représentation matricielle de travail à partir de la représentation matricielle modifiée et de la représentation filtrante, b2) pour chaque bloc diagonal de la représentation matricielle modifiée : b2a) déterminer si le bloc diagonal courant de la représentation matricielle modifiée vérifie une condition de récursivité, b2b) si la condition de récursivité est vérifiée, calculer un bloc courant d'une matrice auxiliaire en réitérant l'opération b) avec le bloc diagonal courant comme représentation matricielle modifiée, et avec un sous-ensemble de la représentation matricielle filtrante modifiée tiré de l'indice i (t;) comme représentation matricielle filtrante modifiée, b2c) si la condition de récursivité n'est pas vérifiée, définir le bloc courant de la matrice auxiliaire comme égal au bloc diagonal courant, b3) calculer des blocs d'une matrice d'approche diagonale dont chaque bloc non nul d'indice i est contraint à vérifier avec le bloc diagonal d'indice i de la matrice auxiliaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux respectivement ladite représentation matricielle filtrante modifiée ou sa transposée et respectivement un bloc d'une matrice triangulaire supérieure par blocs ou un bloc d'une matrice triangulaire inférieure par blocs, et comportant respectivement l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire, et ledit bloc d'indice i de la matrice d'approche, lesdites matrice triangulaire supérieure par blocs et matrice triangulaire inférieure par blocs étant tels que seuls les blocs de la dernière colonne et non diagonaux et respectivement la dernière ligne et non diagonaux sont non nuls et définis égaux aux blocs correspondant de la représentation matricielle modifiée, b4) calculer une matrice diagonale par blocs dont les blocs sont égaux aux blocs de même indice de la matrice auxiliaire, sauf pour le dernier, qui est défini comme la différence entre le dernier bloc de la matrice auxiliaire et la somme pour un indice k non nul et inférieur au nombre de blocs diagonaux de la représentation matricielle modifiée de produits matriciels sous la forme XYZ où X est le bloc non nul de la k-ème colonne de la matrice triangulaire inférieure par blocs, Y est le k-ème bloc de la matrice d'approche, et Z est le bloc non nul de la k-ème ligne de la matrice triangulaire supérieure par blocs, et b5) calculer la représentation matricielle de travail à partir d'un produit matriciel dont un premier terme est égal à la somme de la matrice triangulaire inférieure par blocs et la matrice de l'opération b4), un deuxième terme est l'inverse de la matrice de l'opération b4), et un troisième terme est égal à la somme de la matrice triangulaire supérieure par blocs et la matrice de l'opération b4). L'opération c) comprend travailler récursivement sur la représentation matricielle de travail 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 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 de matrice obtenue après une première opération de renumérotation selon une opération de la figure 3, - la figure 5 représente un exemple de matrice obtenue après une deuxième opération de renumérotation selon la même opération, et - la figure 6 représente un exemple d'une fonction de calcul d'un préconditionneur selon une opération de la figure 3 à partir du résultat obtenu après l'opération dont le résultat est représenté sur la figure 5. Les dessins et la description ci-après contiennent, pour l'essentiel, des éléments de 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 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 10 plupart n'admettent que des solutions par discrétisation et méthode numérique de type différences finies ou volumes finis. 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 15 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 de A pour calculer x. Mais l'inversion des matrices ou l'utilisation de méthodes directes du type élimination de Gauss sont des problèmes complexes, qui monopolisent des 20 puissances de calcul qui croissent de manière exponentielle avec la taille de la matrice. Pour cela, des méthodes itératives basées sur des sous-espaces de Krylov, comme 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 25 matrices qui permettent d'approcher rapidement la matrice inverse de A. En utilisant le préconditionneur M, la méthode itérative résout le système linéaire M-1Ax=M-1 y. Dans ce mode de résolution, on calcule des opérations de type M-1 v et A v, où v est un 30 vecteur, sans calculer de manière explicite l'inverse de M. 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, les méthodes les plus connues pour produire un préconditionneur satisfaisant une propriété de filtrage utilisent des matrices A qui sont tridiagonales par blocs. 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 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. Les Demandeurs ont également conçu une méthode utilisant tout type de matrice A pour produire un préconditionneur satisfaisant une propriété de filtrage. Cette méthode repose sur la factorisation de la matrice A sous la forme d'un produit LDU. Cependant, la parallélisation des calculs au sein de cette méthode peut être améliorée, et cette méthode n'est pas imbriquée. La figure 1 représente un dispositif informatique de calcul polyvalent 2 selon 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 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 ê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. 10 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 15 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. 20 25 30 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, l'adaptateur 10 renumérote les éléments de la matrice A pour permettre un traitement ultérieur en parallèle. Par renumérotation, on entend la réécriture des équations qui définissent la matrice A, de manière à lui donner une forme plus facile à manipuler. Dans la pratique, cela se traduit par la détermination d'une matrice de permutation, qui permet de passer de la forme originale de la matrice A à sa forme renumérotée. L'opération 320 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 renumérotée qui comprend des blocs nuls. Ensuite, dans une opération 340, l'adaptateur 10 traite la matrice A renumérotée, 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. 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. Les Demandeurs ont développé un dispositif mettant en oeuvre un préconditionneur et une méthode de calcul imbriqué d'une représentation de ce préconditionneur qui sont particulièrement adaptés à la parallélisation des calculs. 15 Pour mieux comprendre cela, il convient d'abord 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 20 matrice B en forme de "flèche". Pour cela, la matrice A est renumérotée une première fois pour lui donner la forme de la matrice représentée sur la figure 4, puis les sous-matrices de la matrice de la figure 4 sont elles-mêmes renumérotées de la même manière. 25 La matrice de la figure 4 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. Cette renumérotation de la matrice A est possible du fait de sa faible densité. La même renumérotation est réalisée sur le vecteur x, y et t. 30 Une fois que la matrice de figure 4 a été calculée, il est possible de réappliquer cette même renumérotation aux blocs B1 et B2 et B3. La figure 5 représente une matrice B dans laquelle la renumérotation a été réalisée sur les blocs Bl et B2. On remarquera sur 10 cette figure que le bloc B3 de la figure 4 est le bloc B77 de la figure 5, et que les blocs B4 et B5 de la figure 4 correspondent respectivement aux blocs B71 à B76 et aux blocs B17 à B67 de la figure 5. Une fois que l'opération 320 est terminée, la matrice A a donc été renumérotée d'une manière telle qu'elle présente la forme représentée sur la figure 4. Cependant, chaque bloc diagonal de la figure 4 présente également la même forme de flèche que la matrice A, comme cela est représenté sur la figure 5. L'opération de renumérotation de l'opération 320 pourrait à nouveau être répétée sur les blocs diagonaux de la figure 5, puis sur les blocs diagonaux résultants, jusqu'à ce qu'elle ne produise plus de domaines indépendants. Cette propriété est importante, car comme on le verra plus bas, un bloc diagonal peut être remplacé par un préconditionneur qui l'approche selon la méthode de l'invention, ce qui permet d'imbriquer les calculs, et de réduire ainsi la taille des éléments de travail, tout en augmentant la parallélisation des calculs. Le calcul du préconditionneur M part d'une matrice A sous la forme de la formule (30) 20 de l'annexe A. Les formules (40), (50) et (60) de l'annexe A donnent la composition des éléments respectifs composant la matrice A. Cette forme de matrice correspond à celle de la matrice B de la figure 4, forme qui se retrouve pour tous les blocs diagonaux pour lesquels la renumérotation a été réalisée 25 comme expliqué plus haut. Les Demandeurs ont découvert qu'à partir d'une matrice respectant la formule (30) et de sa factorisation exacte selon la formule (70), il est possible de construire par similarité un préconditionneur M selon la formule (80) de l'Annexe A. Pour cela, la matrice S est définie similairement à la matrice T par la formule (90) de l'Annexe A, où la matrice F représente une approximation de l'inverse de la matrice T. 30 Les compositions des matrices F et C de la formule (90) sont respectivement données dans les formules (100) et (110) de l'Annexe A. Pour que le préconditionneur M satisfasse l'équation (20), il suffit que les matrices S, F, et C vérifient les deux conditions exprimées dans la formule (120) de l'Annexe A. Le développement de ces deux conditions à partir des formules (100) et (110) conduit aux conditions exprimées dans la formule (130) de l'Annexe A. La première condition est fondamentale, car elle exprime le fait que la matrice C doit satisfaire la condition de filtrage par rapport à la matrice D. Ainsi, il est possible de choisir C = D. Mais lorsqu'un bloc diagonal Di, a la même forme en flèche que la matrice de la formule (30), il est également possible d'appliquer à nouveau la méthode, et utiliser un préconditionneur C1 qui approche D' et qui satisfait la condition de filtrage. Ainsi, cette première condition permet de rendre imbriqué le calcul du préconditionneur. La deuxième condition permet de simplifier la condition de la formule (70). En effet, dans celle-ci, la matrice T est définie par rapport à son inverse, ce qui la rend complexe à déterminer. L'utilisation de la matrice F permet de s'affranchir de ce problème. Le fait que F vérifie la deuxième condition de la formule (130) peut être vu comme une condition d'équivalence de chaque bloc F1 avec l'inverse du bloc diagonal Cil. La formule (140) est un premier mode de calcul de la matrice F, et exprime le fait qu'elle regroupe des termes Fl, qui approchent le bloc diagonal C;,-1 pour le N-ème composant de filtrage t. Par j-ème composant, on entend l'ensemble des éléments de t dont les indices correspondent aux indices de la matrice Par exemple, si les matrices D11 à DG_1)G_1) ont p colonnes, et que la matrice Di1 a q colonnes, alors le j-ème composant comprend tous les termes d'indices compris entre p+l et p+q. Lorsque le vecteur (UiNtN) n'a aucune composante nulle, le calcul de Fi; est aisé. La formule (140) permet de tenir compte des cas où ce vecteur présente des composantes nulles. Les Demandeurs ont également découvert un autre calcul pour la matrice Fii, dans laquelle, la matrice Fii est remplacée par une matrice Gii, qui est définie avec la formule (150) de l'annexe A. Pour calculer la matrice Gii, on calcule d'abord la matrice Fii correspondant grâce à la formule (140), puis on applique la formule (150). Dans l'état actuel des recherches des Demandeurs, la formule (140) est préférée à la formule (150), pour des raisons de stabilité. Il est également possible de définir la matrice F satisfaisant la condition de filtrage (130) par la formule (160) combinée à la formule (170) de l'Annexe A, qui est 15 issue des méthodes de déflation de l'algèbre linéaire. Si on analyse les formules (80) et (90) en combinaison avec les conditions de la formule (130) de l'Annexe A, il apparaît donc que le calcul de M dépend du calcul des matrices C et F, que le calcul de C peut impliquer une répétition de la méthode ou une 20 copie d'une partie de la matrice D, et que les composantes de la matrice F peuvent être calculées de manière indépendante, c'est-à-dire en parallèle. Plusieurs moyens de parallélisation apparaissent donc : - l'exécution des diverses boucles imbriquées peuvent être réalisées en parallèle (par 25 exemple celle du bloc B1 et celle du bloc B2 de la figure 4, et ainsi de suite avec les blocs B11, B22, B33, B44, B55, B66 et B77 de la figure 5), - au sein d'une même boucle, les calculs des composantes de la matrice F peuvent également être réalisés en parallèle, et - pour le calcul du terme SNN, les calculs peuvent également être réalisés en parallèle. 30 La figure 6 représente un exemple de fonction NFFQ permettant de calculer le préconditionneur M. Pour cela, chaque terme de la matrice est calculé, et affecté aux matrices C,F et S comme il convient. Dans une opération 600, la fonction NFFQ reçoit une matrice B et un vecteur de filtrage t comme entrées, et utilise comme paramètres un nombre N qui est le nombre de blocs diagonaux dans la matrice B, ainsi qu'un indice i initialisé à 0. Comme on l'a vu plus haut, dans le cadre de la figure 3, la matrice B est la matrice A renumérotée, qui est utilisée comme entrée pour la fonction NFFQ. Dans ce qui suit, les mots "élément" ou "terme", doivent être compris comme désignant un bloc. En effet, la matrice B est prédécoupée en blocs rectangulaires. Seuls les blocs diagonaux de B doivent être carrés. Les valeurs de côtés des blocs sont définies par les dimensions respectives des blocs diagonaux de la matrice B, et sont données par la procédure de renumérotation de la matrice A. Ensuite, une première boucle est exécutée afin de lancer toutes les instances imbriquées de la fonction NFFQ. Pour cela, l'indice i est incrémenté, dans une opération 602, puis une fonction Nest() détermine dans une opération 604 si le bloc diagonal Dii correspondant a une forme semblable à celle de la formule (30) de l'Annexe A. Cela définit une condition de récursivité. Lorsque ce n'est pas le cas, alors le terme Ci; est défini comme identique au terme Dii dans une opération 606. Lorsque le terme Dii a une forme semblable à celle de la formule (30) de l'Annexe A, cela signifie que la fonction NFFQ peut être utilisée pour calculer Ch, et la fonction NFFQ est à nouveau instanciée dans une opération 608, avec comme entrées le bloc Dii et le sous-vecteur ti.30 Ensuite, un test vérifie dans une opération 610 si l'indice i est inférieur à N+1. Si c'est le cas, alors la boucle est relancée en 602. Sinon, la boucle se termine avec le calcul des éléments de la matrice F et de la matrice S. On notera ici que chaque appel à la fonction NFFQ dans l'opération 608 peut être lancé en faisant appel à un nouveau processeur ou à un nouveau coeur de processeur, c'est-à-dire en parallèle avec les calculs liés aux autres appels à cette fonction pour d'autres blocs. De plus, bien que la fonction décrite ici soit présentée de manière itérative, les tests de l'opération 604 sur les N éléments Dii peuvent être réalisés en parallèle, ainsi que les opérations 606 ou 608 qui suivent ces tests. Une fois la boucle terminée, et que toutes les instances de la fonction NFFQ ont été initialisées, chacune de ces instances calcule les (N-1) termes Fii où i varie entre 1 et (N-1), au moyen de la formule (140) de l'Annexe A. Si le deuxième mode de calcul est retenu, le calcul de Gii est également réalisé, conformément à la formule (150). Ici encore, on notera que ces calculs peuvent être réalisés en parallèle, car ils sont indépendant les uns des autres. Une fois ces termes calculés, il ne reste plus qu'à calculer le terme SNN conformément à la formule (180) de l'Annexe A, dans laquelle les termes de la somme peuvent ici encore être calculés en parallèle. Pour les termes S11 à S(N_1)(_1), le développement de la formule (90) montre qu'ils sont chacun égaux au terme Cii de même indice. Enfin, la matrice M qui approche la matrice B en entrée est retournée comme résultat dans une opération 616, sous sa forme décomposée comme dans l'équation (80).. Cela est nécessaire pour permettre de réaliser les calculs dans les fonctions NFFQ d'ordre supérieur. Ainsi, les fonctions NFFQ les plus profondes réalisent leurs calculs, puis les résultats remontent de couche en couche, jusqu'à la première instance, et le calculateur-solveur (12) est appelé. Il est avantageux de conserver la matrice S calculée dans chaque instance de la fonction NFFQ. En effet, la résolution d'un système Mu = v peut être faite en utilisant la décomposition de M selon la formule (80) de l'Annexe A, et en résolvant deux systèmes linéaires successifs par substitution. Chaque bloc de M ayant été calculé de manière imbriquée peut être résolu de la même manière, ce qui permet d'accélérer le calcul en imbriquant la résolution des systèmes linéaires. 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 (190) de l'Annexe A. Comme la formule (190) est presque la transposée de la formule (20), l'utilisation de la formule (190) comme condition de stabilité ne change rien au mode de fonctionnement de l'invention. En conséquence, les formules (120) à (140) n'ont qu'à être légèrement modifiées, comme présenté avec les formules (200) à (220). Les formules (160) et (170) 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. 30 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.25 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 320, et lorsque la matrice A est découpée en blocs. Un minibloc 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 F. En effet, la formule (140) 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. Les Demandeurs ont donc généralisé la formule (140) sous la forme de la formule (230), 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). Une autre manière de voir cette modification est de remarquer que la formule (140) peut être vue comme un cas particulier de la formule (230), dans laquelle t a une seule colonne. La formule (220) pourra être modifiée à l'identique. 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.30 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 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 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 25 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. ANNEXE A Ax=y (10) At = Mt (20) D11 U1N A= (30) L = LN 1 LN2 D.yx (40) L 0 D22 (50) DNN.. 0 = (60) 0 A =(L+ T) T-I(U + T), auecT tel que T = D - LT:"-IU (70) M = (L (80) 5 = C - LFU (90) 21 F(N_ (N-1) (100) F11 F= C.] (110) [(S-1 - F)Ut= 0 (C - D) t = 0 (120) Cuti Dut. avec t iN t)J = Fu tN avec . N - 1 Fit (ri (Ili.,NtN)(r) si r = s et (UivtN)(.r) 0 tNXr)./ ,v) s) si UiNtN)(r) = 0 avec s tel que (U.,, o 0 sinon Gii 2eii ii ii zi = uiNtN I- CiiHi + Hi (170) SiVN CNN - LNiFiiuÈN (180) 20 tTA = tTM (190) 15 (130) (140) (150) (160) EtrL(s-I .F) = 0 tr(C -.D) = 0 tr C.. = tT Du avec t = 1 N' ut m Cr; = t'itrL Ni FI avec i = 1. . N -1 (200) (210) 25 (220) (tLN,)(s) = o avec r tel que (tZLNi r * 0 0 sinon (g LNA71)(r)/ZNi)( si r = s et (t,;,., Njr) * O FÉi(r,$) = (230) IDie(Ç'UiNtN) H/. (UiNtN)(r)) si r s et (UaVtNX'r) n'a aucune composante nulle DagerUiNtNXr)/.(UiNtN)(s» si (UiNtN) (r) = 0 avec s tel que (UNt,N)(s) 0 0 sinon 22 REVENDICATIONS1.- Dispositif informatique de calcul polyvalent pour la modélisation et la simulation de phénomènes 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 phénomène 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 (12), la représentation matricielle de travail (M) étant contrainte à vérifier avec la représentation matricielle initiale (A) une condition de stabilité ((20), (190)) 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 (10) est agencé d'une part pour renuméroter la représentation matricielle initiale (A) et la représentation matricielle filtrante (t) pour produire une représentation matricielle modifiée (B) et une représentation matricielle filtrante modifiée (t) 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 (30), et d'autre part pour calculer récursivernent la représentation matricielle de travail (M) à partir de la représentation matricielle modifiée (B) et de ladite représentation matricielle filtrante modifiée (t), tandis que le calculateur-solveur (12) est agencé pour travailler récursivement sur la représentation matricielle de travail (M) 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 récursif de l'adaptateur (10) comprend le calcul de la représentation matricielle de travail (M) sous la forme d'un produit matriciel PQR, * où P est la somme entre d'une part une matrice diagonale (S) calculée à partir d'une matrice auxiliaire (C) et d'une matrice d'approche (F), et d'autre part une matrice triangulaire inférieure par blocs (L) dont seuls les termes de la dernière ligne de blocs et non diagonaux sont non nuls et sont égaux aux termes de même indice dans la représentation matricielle modifiée (B), * où Q est l'inverse de la matrice diagonale (S),* où R est la somme entre la matrice diagonale (S) et une matrice triangulaire supérieure par blocs (U) dont seuls les termes de la dernière colonne de blocs et non diagonaux sont non nuls, et sont égaux aux termes de même indice dans la représentation matricielle modifiée (B), * la matrice auxiliaire (C) étant une matrice diagonale par blocs, dont chaque bloc d'indice i est : - défini égal au bloc diagonal d'indice i (D;; de la représentation matricielle modifiée (B) lorsque ce bloc ne vérifie pas la condition de récursivité, et - sinon calculé par un appel récursif à l'adaptateur (10) avec le bloc diagonal d'indice i (D;;) de la représentation matricielle modifiée (B) comme représentation matricielle modifiée et avec un sous-ensemble de la représentation matricielle filtrante modifiée tiré de l'indice i (t;) comme représentation matricielle filtrante modifiée, * la matrice d'approche (F) est une matrice diagonale par blocs dont le dernier bloc est nul, et dont chaque bloc non nul d'indice i est contraint à vérifier avec le bloc diagonal d'indice i de la matrice auxiliaire (C) une condition d'équivalence ((130), (210)), comprenant une expression de comparaison de deux produits matriciels comportant tous deux respectivement ladite représentation matricielle filtrante modifiée (t) ou sa transposée et respectivement un bloc de la matrice triangulaire supérieure par blocs (U) ou un bloc de la matrice triangulaire inférieure par blocs (L), et comportant respectivement l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire (C), et ledit bloc d'indice i de la matrice d'approche (F), * les blocs diagonaux de la matrice diagonale (S) étant égaux aux blocs de même indice de la matrice auxiliaire (C), sauf pour le dernier, qui est défini comme la différence entre le dernier bloc de la matrice auxiliaire (C) et la somme pour un indice k non nul et inférieur au nombre de blocs diagonaux de la représentation matricielle modifiée (B) de produits matriciels sous la forme WXY où W est le bloc non nul de la k-ème colonne de la matrice triangulaire inférieure par blocs (L), X est le k-ème bloc de la matrice d'approche (F), et Y est le bloc non nul de la k-ème ligne de la matrice triangulaire supérieure par blocs (U). 2.- Dispositif selon la revendication 1, dans lequel la condition de stabilité ((20)) comprend une comparaison de deux produits matriciels comportant tous deux ladite 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 condition d'équivalence ((130)) comprend une expression de comparaison de deux produits matriciels comportant tous deux ladite représentation matricielle filtrante modifiee (t) et un bloc de la matricetriangulaire supérieure par blocs (U), et respectivement l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire (C) et ledit bloc d'indice i de la matrice d'approche (F). 3.- Dispositif selon la revendication 1, dans lequel la condition de stabilité ((190)) comprend une comparaison de deux produits matriciels comportant tous deux la transposée de ladite 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 condition d'équivalence ((210)) comprend une expression de comparaison de deux produits matriciels comportant tous deux la transposée de ladite représentation matricielle filtrante modifiée (t) et un bloc de la matrice triangulaire inférieure par blocs (L), et respectivement l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire (C) et ledit bloc d'indice i de la matrice d'approche (F). 4.- Dispositif selon l'une des revendications précédentes, dans lequel le bloc d'indice i de la matrice d'approche (F) est calculé à partir d'une division terme à terme faisant intervenir l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire (C), et soit le bloc d'indice i de la représentation matricielle filtrante modifiée (t) et le bloc d'indice i de la matrice triangulaire supérieure par blocs (U), soit le bloc d'indice i de la transposée de ladite matricielle filtrante modifiée (t) et le bloc d'indice i de la matrice triangulaire inférieure par blocs (L). 5.- Dispositif selon lune des revendications 1 à 3, dans lequel la matrice d'approche (F) est calculée par blocs par une méthode de déflation, dans laquelle un premier terme (Z;) fait intervenir le bloc d'indice N de ladite représentation matricielle filtrante modifiée (t) et le bloc non nul de la i-ème ligne de la matrice triangulaire supérieure par blocs (U), dans laquelle un second terme (Hi) fait intervenir le bloc d'indice i de la matrice auxiliaire (C) et le premier terme, le bloc d'indice i de la matrice d'approche étant défini comme la différence entre le deuxième terme (Hi) et le produit matriciel du bloc d'indice i de la matrice auxiliaire (C) avec le deuxième terme (Hi), cette différence étant ajoutée de la matrice identité (I). 6.- Dispositif selon l'une des revendications précédentes, dans lequel la représentation matricielle filtrante (t) est un vecteur colonne. 7.- 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'ensemble de 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. 8.- 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. 9.- Procédé de calcul polyvalent pour la modélisation et la simulation de phénomènes physiques, du type comprenant : a) recevoir une représentation matricielle initiale correspondant à un système d'équations à traiter représentant un phénomène 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 : b1) renuméroter la représentation matricielle initiale et la représentation matricielle filtrante pour produire une représentation matricielle modifiée et une représentation matricielle filtrante modifiée, et calculer récursivement la représentation matricielle de travail à partir de la représentation matricielle modifiée et de la représentation filtrante modifiée, b2) pour chaque bloc diagonal de la représentation matricielle modifiée : b2a) déterminer si le bloc diagonal courant de la représentation matricielle modifiée vérifie une condition de récursivité, b2b) si la condition de récursivité est vérifiée, calculer un bloc courant d'une matrice auxiliaire en réitérant l'opération b) avec le bloc diagonal courant comme représentation matricielle modifiée, et avec un sous-ensemble de la représentation matricielle filtrante modifiée tiré de l'indice i (t;) comme représentation matricielle filtrante modifiée, b2c) si la condition de récursivité n'est pas vérifiée, définir le bloc courant de la matrice auxiliaire comme égal au bloc diagonal courant, b3) calculer des blocs d'une matrice d'approche diagonale dont chaque bloc non nul d'indice i est contraint à vérifier avec le bloc diagonal d'indice i de la matrice auxiliaire une condition d'équivalence, comprenant une expression de comparaison de deux produits matriciels comportant tous deux respectivement ladite représentation matricielle filtrante modifiée ou sa transposée et respectivement un bloc d'une matricetriangulaire supérieure par blocs ou un bloc d'une matrice triangulaire inférieure par blocs, et comportant respectivement l'inverse dudit bloc diagonal d'indice i de la matrice auxiliaire, et ledit bloc d'indice i de la matrice d'approche, lesdites matrice triangulaire supérieure par blocs et matrice triangulaire inférieure par blocs étant tels que seuls les blocs de la dernière colonne et non diagonaux et respectivement de la dernière ligne et non diagonaux sont nuls et définis égaux aux blocs correspondant de la représentation matricielle modifiée, b4) calculer une matrice diagonale par blocs dont les blocs sont égaux aux blocs de même indice de la matrice auxiliaire, sauf pour le dernier, qui est défini comme la différence entre le dernier bloc de la matrice auxiliaire et la somme pour un indice k non nul et inférieur au nombre de blocs diagonaux de la représentation matricielle modifiée5 de produits matriciels sous la forme XYZ où X est le bloc non nul de la k-ème colonne de la matrice triangulaire inférieure par blocs, Y est le k-ème bloc de la matrice d'approche, et Z est le bloc non nul de la k-ème ligne de la matrice triangulaire supérieure par blocs, et b5) calculer la représentation matricielle de travail à partir d'un produit matriciel dont un premier terme est égal à la somme de la matrice triangulaire inférieure par blocs et la matrice de l'opération b4), un deuxième terme est l'inverse de la matrice de l'opération b4), et un troisième terme est égal à la somme de la matrice triangulaire supérieure par blocs et la matrice de l'opération b4), et en ce que l'opération c) comprend travailler récursivement sur la représentation matricielle de travail de façon à fournir une solution du système d'équations de la représentation matricielle initiale sans inversion complète de celle-ci.25