d4.06.40 Un exemple de matrice distribuée pour PETSc#
Résumé :
Pour utiliser les solveurs de la librairie PETSc en parallèle, on doit fournir une matrice assemblée distribuée sur les processeurs disponibles selon un modèle de distribution imposé par PETSc (répartition par groupe de lignes contiguës). Le modèle de distribution des données Code_Aster repose sur une répartition des mailles entre les processeurs. Il conduit à une matrice assemblée qui ne suit pas le modèle de distribution PETSc. On explique ici en suivant un exemple physique simple les différentes numérotations des degrés de liberté: numérotation locale et globale pour Code_Aster puis numérotation PETSc.
Problème modèle#
On considère un problème 2D d’élasticité statique, distribué sur deux processeurs.
Géométrie#
Le domaine élastique est un carré de points extrêmes (-50,50) et (50,50). On maille ce domaine en 4 quadrangles.
|
|
|
|
Propriétés du matériau#
\(E=1,0{10}^{11}N/{m}^{2}\)
\(\nu =0,3\)
Modélisation 2D#
Le modèle est affecté sur les 4 mailles quadrangulaires (groupe de mailles 2D all) et sur les 2 segments qui constituent le bord supérieur du domaine (groupe de mailles 1D up).
MODEL=AFFE_MODELE(MAILLAGE=MA,
AFFE=_F(GROUP_MA=('all','up',),
PHENOMENE='MECANIQUE',
MODELISATION='C_PLAN',),
PARALLELISME=_F(DISTRIBUTION='SOUS_DOMAINE',),
);
On a défini également un type de répartition des mailles (par sous-domaine).
Conditions limites et chargement#
On exerce une pression répartie sur le bord supérieur du carré, d’extrémités (-50,50) (50,50):
PRESSION=AFFE_CHAR_MECA(MODELE=MODEL,
PRES_REP=_F(GROUP_MA='up',
PRES=1000000000,),);
La base du carré est encastrée : on applique la condition DX =0, DY =0 sur le segment d’extrémités \((-50,-50)\) \((50,-50)\) .
Cet encastrement est appliqué de deux façons:
modélisation A: avec AFFE_CHAR_CINE
ENCASTR=AFFE_CHAR_CINE(MODELE=MODEL, MECA_IMPO=_F(GROUP_NO='bottom', DX=0, DY=0,),);
modélisation B: avec AFFE_CHAR_MECA
ENCASTR=AFFE_CHAR_MECA(MODELE=MODEL, DDL_IMPO=_F(GROUP_NO='bottom', DX=0, DY=0,),);
Distribution du problème#
Lors de la création du modèle, on a choisi une répartition des calculs élémentaires par sous-domaine. On précise au solveur que la matrice est distribuée:
MECA_STATIQUE(MODELE=MODEL,
CHAM_MATER=AFMAT,
EXCIT=(_F(CHARGE=PRESSION,),
_F(CHARGE=ENCASTR,),),
SOLVEUR=_F(METHODE=”PETSC”,
MATR_DISTRIBUEE=”OUI”,
ALGORITHME=”GMRES”,),);
Le calcul est exécuté sur 2 processeurs.
Généralités sur le système linéaire#
On veut résoudre un problème d’élasticité linéaire. On rappelle sa formulation variationnelle au paragraphe 3.1 .
Le problème discret associé au problème continu est différent selon la méthode utilisée pour appliquer la condition limite de Dirichlet sur la base du domaine. On le décrira ultérieurement pour chaque modélisation. On regroupe dans cette section la partie commune aux modélisations A et B.
Problème continu#
Les contraintes vérifient l’équation \(-\nabla \cdot \sigma =0\) dans le domaine élastique. Le déplacement \(u\) est fixé sur la frontière inférieure du domaine \({\Gamma}_{D}\) par la condition de Dirichlet \(u={u}_{0}\) . Ici, \({u}_{0}=0\) . Sur la frontière supérieure \({\Gamma}_{N}\) , on applique le chargement de Neumann \(\sigma \cdot n=p\) , pour une pression \(p\) uniforme.
On écrit la formulation variationnelle du problème. On définit:
l’espace affine \(V=\lbrace v,v\in {({H}^{1}(\Omega ))}^{2,}v={u}_{0}\mathit{sur}{\Gamma}_{D}\rbrace\) ,
l’espace vectoriel associé (déplacements cinématiquement admissibles) \({V}^{0}=\lbrace v,v\in {({H}^{1}(\Omega ))}^{2,}v=0\mathit{sur}{\Gamma}_{D}\rbrace\)
On cherche \(u\in V\) , tel que \(\forall v\in {V}^{0}\) , \(a(u,v)=(f,v)\) , où \(a\) est la forme bilinéaire de l’élasticité linéaire et \(f\) correspond au chargement «pression» sur \({\Gamma}_{N}\) .
Numérotation des équations de la matrice Code_Aster globale#
La description de la numérotation de la matrice (globale et locale) est faite par les objets du NUME_DDL. Le NUME_EQUA contient les informations globales. On rappelle rapidement le rôle des champs dont on va afficher la valeur et on renvoie à [D4.06.07] pour la description de référence.
.NUME.REFN est un tableau de noms, contenant le nom du maillage, de la grandeur discrétisée et du solveur. Ici REFN contient:
MA |
DEPL_R |
GMRES |
.NUME.NEQU est le nombre total d’équations \(N\)
.NUME.DELG est un tableau de taille \(N\) , qui contient un marqueur des degrés de liberté de Lagrange avec la convention:
\(\mathit{DELG}(i)=\lbrace \begin{array}{ccc}0& :& \mathit{ddl}\mathit{physique}\\ -1& :& \mathit{Lagrange}1\\ -2& :& \mathit{Lagrange}2\end{array}\)
.NUME.DEEQ est un tableau entier de taille \(\mathrm{2N}\) L’équation ieq est décrite par un couple d’entiers ( ino , icmp ) en positions ( ieq -1) × 2+1 et ( ieq -1) × 2+2 du tableau:
ino >0 est le numéro du nœud du maillage support du degré de liberté ieq .
Si icmp >0, le degré de liberté correspond à la composante icmp de la grandeur sous-jacente.
Si icmp <0, alors l’équation est une des deux équations de dualisation du blocage de la composante icmp de cette grandeur sur le nœud ino >0.
ino =0 indique que l’équation ieq est l’équation de dualisation d’une relation linéaire. Dans ce cas, icmp est forcément nul.
.NUME.PRNOest une collection de taille 1 + le nombre de charges dualisées dans le modèle. On accède à ces objets par le pointeur de noms .LILI. Le premier objet de la collection contient les degrés de liberté associés au maillage (par convention .LILI(1)=&MAILLA) et les objets suivants les degrés de liberté tardifs créés par dualisation. Pour chaque nœud, on trouve
le numéro d’équation de la première composante associée à ce nœud,
le nombre de composantes portées par ce nœud
un vecteur d’entiers codés précisant quelles sont les composantes de la grandeur présentes sur ce nœud (voir D4.06.05 pour le système de codage)
Remarque :
Attention, il y a un niveau supplémentaire d’indirection (vecteur .NUEQ ) qui permet de trouver la position des valeurs des termes de la matrice dans le tableau .VALE *.* **On omet volontairement ici ce niveaucar il n’intervient pas dans l’exemple considéré. Il est seulement nécessaire* pour la sous-structuration statique. En effet, on a alors besoin d’assurer que les degrés de liberté d’interface se trouvent dans la matrice après les degrés de liberté internes à un sous-domaine. Cela est nécessaire pour calculer le complément de Schur sans avoir explicitement une structure de matrice par blocs. La sous-structuration statique correspond à la commande MACR_ELEM_STAT. Un macro-élément est équivalent à un sous-domaine. On calcule entre autres sa rigidité et il peut être intégré comme un élément à un modèle.
Numérotation des équations de la matrice Aster locale#
On se place maintenant dans le cas où la matrice Code_Aster est assemblée par deux processeurs.
Dans notre exemple, le processeur 0 possède les mailles M7, M8 et M10 et le processeur 1 les mailles M9, M11 et M12. Chaque processeur assemble les contributions élémentaires correspondant aux mailles qu’il possède. Les caractéristiques de ces matrices locales sont définies dans le .NUML du NUME_DDL. Les tableaux .NUML.NULG et .NUML.NUGL sont inverses l’un de l’autre. Ils contiennent la correspondance entre la numérotation locale à chaque processeur des degrés de liberté et la numérotation globale. Ainsi, si \({i}_{l}\) est le numéro local d’un degré de liberté et \({i}_{g}\) son numéro global, on a:
\(\mathrm{NULG}({i}_{l})={i}_{g},\mathrm{NUGL}({i}_{g})={i}_{l}\)
Si le degré de liberté \({i}_{g}\) n’est pas présent sur un processeur, \(\mathrm{NUGL}({i}_{g})=0\) .
Le système linéaire pour la modélisation A#
Les conditions limites sont appliquées par élimination. On explique le principe de l’élimination dans le paragraphe 4.1 , et on donne les valeurs des différents descripteurs de la matrice dans le paragraphe 4.2 .
Le système discret#
Le système discret s’écrit \(KU=F\) où la matrice de rigidité \(K\) est de taille \(N\times N\) .
On partitionne l’ensemble des indices des degrés de liberté \(S=⟦1:N⟧\) , en deux sous-ensembles:
\(D\) pour les indices des degrés de liberté fixés par la condition limite de Dirichlet (encastrement). Ici, \(D=\lbrace 1,2,5,6,11,12\rbrace\) , ce qui correspond à fixer \(\mathit{DX},\mathit{DY}\) sur les nœuds \(\mathit{N1},\mathit{N3},\mathit{N6}\) et \({U}_{0}=0\) .
\(I\) pour les indices des degrés de liberté non contraints, inconnus: \(I=S\setminus D\)
On applique cette partition au système qui s’écrit:
\(\left(\begin{array}{cc}{K}_{II}& {K}_{ID}\\ {K}_{DI}& {K}_{DD}\end{array}\right)\left(\begin{array}{c}{U}_{I}\\ {U}_{D}\end{array}\right)=\left(\begin{array}{c}f\\ 0\end{array}\right)\) , avec \({U}_{D}={U}_{0}\) .
On remplace ce système par:
\(\left(\begin{array}{cc}{K}_{II}& {0}_{\mathit{ID}}\\ {0}_{\mathit{DI}}& {\mathit{Id}}_{DD}\end{array}\right)\left(\begin{array}{c}{U}_{I}\\ {U}_{D}\end{array}\right)=\left(\begin{array}{c}f-{K}_{\mathit{ID}}{U}_{0}\\ {U}_{0}\end{array}\right)\) .
Ce qui correspond aux mises à jours suivantes à faire sur la matrice et le second membre:
\(\begin{array}{ccc}{F}_{I}& \leftarrow & f-{K}_{\mathit{ID}}{U}_{0}\\ {F}_{D}& \leftarrow & {U}_{0}\\ {K}_{\mathit{DD}}& \leftarrow & {\mathit{Id}}_{\mathit{DD}}\\ {K}_{\mathit{ID}}& \leftarrow & {0}_{\mathit{ID}}\\ {K}_{\mathit{DI}}& \leftarrow & {0}_{\mathit{DI}}\end{array}\) .
Les descripteurs de la matrice#
Matrice globale
Il n’y a pas de degré de liberté de Lagrange: le tableau .NUME.DELG est rempli de 0. Le tableau .NUME.DEEQ est de taille \(\mathrm{2N}=36\) . Voilà son contenu (1 ere ligne) et son interprétation (2 eme ligne):
La collection .PRNO n’a qu’un objet: les seuls degrés de liberté sont les degrés de liberté physiques associés aux nœuds du maillage. Voilà son contenu nœud par nœud:
1 2 6 0 0 0 0 # nœud N1
3 2 6 0 0 0 0 # nœud N2
5 2 6 0 0 0 0 # nœud N3
7 2 6 0 0 0 0 # nœud N4
9 2 6 0 0 0 0 # nœud N5
11 2 6 0 0 0 0 # nœud N6
13 2 6 0 0 0 0 # nœud N7
15 2 6 0 0 0 0 # nœud N8
17 2 6 0 0 0 0 # nœud N9
Ainsi la première équation posée sur le nœud N3 du maillage porte le numéro d’équation 5. Il y a 2 équations posées sur ce nœud. Pour coder les composantes de la grandeur DEPL_R on utilise un vecteur de taille 5. Ce vecteur a comme valeurs
Seul le premier terme est non nul: il vaut \(6={2}^{1}+{2}^{2}\) .
Cela qui signifie que les composantes sélectionnées de la grandeur sont les deux premières dans l’ordre du catalogue, c’est à dire les composantes DX et DY.
Matrice distribuée#
On reporte ensuite les valeurs du .NUML pour les processeurs 0 et 1.
La correspondance inverse est donnée par le tableau .NUML.NUGL:
Grâce au .DEEQ on peut retrouver quels nœuds sont supports des degrés de liberté sur chaque processeur:
On voit que les composantes sur un même nœud sont numérotées de façon contiguës.
Remarque :
Lorsque la matrice est distribuée, les équations de blocage cinématiques sont multipliées par le nombre de processeurs utilisés pour le calcul. La contrainte \(U={U}_{0}\) est remplacée sur p processeurs par \(\mathit{pU}={\mathit{pU}}_{0}\) . Dans la routine asmchc , on met des 1 sur la diagonale de la matrice et cette opération est effectuée sur tous les processeurs. Ces termes s’accumulent dans la matrice globale.
Le système linéaire pour la modélisation B#
Les conditions limites sont dualisées. La méthode de (double) dualisation est décrite dans la documentation [R3.03.01]. On précise le système discret au paragraphe 5.1 et la valeur des descripteurs de la matrice au paragraphe 5.2 .
Le système discret#
La condition limite sur \({\Gamma}_{D}\) est imposée comme une contrainte sur les déplacements, par des multiplicateurs de Lagrange. On note:
\(U\) les degrés de liberté «physiques»;
\(\Lambda\) les multiplicateurs de Lagrange.
On résout un système linéaire augmenté, de la forme:
\(\left(\begin{array}{cc}{K}_{\mathit{UU}}& {C}_{U\Lambda }^{T}\\ {C}_{\Lambda U}& {D}_{\Lambda \Lambda }\end{array}\right)\left(\begin{array}{c}U\\ \Lambda \end{array}\right)=\left(\begin{array}{c}f\\ {u}_{0}\end{array}\right)\)
où \({C}_{\Lambda U}\) est la matrice des contraintes.
Les descripteurs de la matrice assemblée#
Matrice globale#
Le nombre d’équations du problème est
.
Calcul du nombre total de degrés de liberté :
Il y a 9 nœuds et chaque nœud porte les composantes DX et DY de la grandeur DEPL_R, soit au total 18 équations portant sur 18 degrés de liberté « physiques ».
A ces degrés de liberté « physiques », on ajoute des degrés de liberté de type « Lagrange » qui permettent de prendre en compte par dualisation la condition limite d’encastrement. Il y a 3 nœuds encastrés dans le problème. Sur chaque nœud encastré, on bloque deux composantes DX et DY. Pour chaque composante bloquée, on utilise 2 multiplicateurs de Lagrange. Au total, on définit donc pour ce problème \(3\times 2\times 2=12\) degrés de liberté de type Lagrange.
Valeurs du NUME_EQUA:
Dans le .PRNO, le premier objet décrit les nœuds du maillage, qui portent des degrés de liberté physique (comme dans la modélisation A):
3 2 6 0 0 0 0 # nœud 1
7 2 6 0 0 0 0 # nœud 2
11 2 6 0 0 0 0 # nœud 3
15 2 6 0 0 0 0 # nœud 4
17 2 6 0 0 0 0 # nœud 5
21 2 6 0 0 0 0 # nœud 6
25 2 6 0 0 0 0 # nœud 7
27 2 6 0 0 0 0 # nœud 8
29 2 6 0 0 0 0 # nœud 9
Le second décrit l’encastrement:
1 1 134217728 0 0 0 0
5 1 134217728 0 0 0 0
2 1 134217728 0 0 0 0
6 1 134217728 0 0 0 0
9 1 134217728 0 0 0 0
13 1 134217728 0 0 0 0
10 1 134217728 0 0 0 0
14 1 134217728 0 0 0 0
19 1 134217728 0 0 0 0
23 1 134217728 0 0 0 0
20 1 134217728 0 0 0 0
24 1 134217728 0 0 0 0
Chaque nœud porte une seule composante qui est déterminée par l’entier 134217728=227, c’est donc la composante LAGR dont il s’agit (multiplicateur de Lagrange).
Le .DEEQ donne pour chaque équation le nœud support du degré de liberté et la composante de la grandeur. On affiche ses valeurs dans le tableau suivant. Lorsque le degré de liberté est un degré de liberté physique, (DX,DY) indique quelle composante de la grandeur il discrétise. Pour les degrés de liberté de type Lagrange, (Lag DX,Lag DY) indique quelle composante est bloquée par le multiplicateur de Lagrange.
Remarque:
Les deux degrés de liberté de Lagrange qui permettent de bloquer une composante sur un degré de liberté physique doivent être numérotés de façon à encadrer le degré de liberté physique : l’un doit se trouver avant et l’autre après. Dans Code_Aster, l’algorithme de numérotation assure cette contrainte mais les degrés de liberté de Lagrange ne sont pas forcément numérotés au plus près de la composante bloquée. C’est le cas sur cet exemple mais ce n’est pas une obligation.
Matrice distribuée#
La partition est la même que dans la modélisation A. Les degrés de liberté de Lagrange sont tous affectés au processeur 0. On décrit la numérotation de la matrice dans le tableau suivant:
Dans ce tableau, la valeur -1 indique qu’il s’agit d’un Lagrange 1 et la valeur -2 qu’il s’agit d’un Lagrange 2.
Remarque :
Le .DEEQ permet aussi de savoir sur quels degrés de liberté sont des multiplicateurs de Lagrange, mais pas s’il s’agit d’un Lagrange 1 ou 2. La valeur -1 indique que l’on bloque la composante DX et la valeur -2 que l’on bloque la composante DY . Pour reconstruire toute l’information, il faut utiliser les deux objets (.DELG et .DEEQ ).
Construction de la matrice PETSc#
Modèle de données PETSc#
PETSc propose un type de données Mat. Ce type de données suppose, en parallèle, que la matrice est répartie sur les processeurs disponibles par groupe de lignes contiguës.
Il faut donc redistribuer la matrice Code_Aster . On construit une nouvelle numérotation globale pour PETSc. Dans cette numérotation, un degré de liberté qui n’appartient qu’à un processeur est affecté à ce processeur et un degré de liberté partagé par plusieurs processeurs appartient au processeur de rang inférieur.
Remarque :
PETSc stocke la matrice en entier (pas de stockage symétrique).
PETSc utilise la convention C d’indiçage des tableaux (de 0 à N-1 pour un tableau de taille N).
Descripteur de matrice pour la modélisation B#
La correspondance entre la numérotation locale de la matrice et sa numérotation globale pour PETSc est donnée par le tableau NUML.NLGP. On indique dans le tableau suivant les trois numérotations pour la modélisation B: \({i}_{l}\) est l’indice local, \({i}_{g}\) l’indice global pur la numérotation Code_Aster et \({i}_{\mathit{gp}}\) l’indice global pour PETSc.
