Calcul de la structure clairsemée pour les matrices d'éléments finis

Question: Quelles méthodes sont disponibles pour calculer avec précision et efficacité la structure de rareté d'une matrice d'éléments finis?

Info: Je travaille sur un solveur d'équation de pression de Poisson, en utilisant la méthode de Galerkin avec une base de Lagrange quadratique, écrite en C, et en utilisant PETSc pour le stockage matriciel clairsemé et les routines KSP. Pour utiliser PETSc efficacement, j'ai besoin de pré-allouer de la mémoire pour la matrice de rigidité globale.

Actuellement, je fais un faux assemblage pour estimer le nombre de non-zéros par ligne comme suit (pseudocode)

int nnz[global_dim]
for E=1 to NUM_ELTS
  for i=1 to 6
    gi = global index of i 
    if node gi is free
      for j=1 to 6
        gj = global index of j
        if node gj is free 
          nnz[i]++

Cela surestime toutefois nnz car certaines interactions nœud-nœud peuvent se produire dans plusieurs éléments.

J'ai envisagé de garder une trace des interactions i, j que j'ai trouvées, mais je ne sais pas comment faire cela sans utiliser beaucoup de mémoire. Je pourrais également faire une boucle sur les nœuds et trouver le support de la fonction de base centrée sur ce nœud, mais je devrais alors rechercher tous les éléments de chaque nœud, ce qui semble inefficace.

J'ai trouvé cette question récente, qui contenait des informations utiles, en particulier de Stefano M, qui a écrit

mon conseil est de l'implémenter en python ou C, en appliquant certains concepts théoriques des graphes, c'est-à-dire considérer les éléments de la matrice comme des arêtes dans un graphe et calculer la structure de rareté de la matrice d'adjacence. La liste des listes ou le dictionnaire des clés sont des choix courants.

Je cherche plus de détails et de ressources à ce sujet. Certes, je ne connais pas beaucoup la théorie des graphes, et je ne connais pas toutes les astuces CS qui pourraient être utiles (j'aborde cela du côté mathématique).

Merci!

Votre idée de garder une trace des interactions i, j que vous avez trouvées peut fonctionner, je pense que c'est le "truc CS" auquel vous et Stefano M faites référence. Cela revient à construire votre matrice clairsemée au format liste de listes .

Je ne sais pas combien de CS vous avez, donc je m'excuse si cela vous est déjà connu: dans une structure de données de liste liée , chaque entrée stocke un pointeur vers l'entrée suivante et l'entrée précédente. Il n'est pas cher d'ajouter et de supprimer des entrées, mais pas aussi simple que d'y trouver des éléments - vous devrez peut-être les parcourir tous.

Ainsi, pour chaque nœud i, vous stockez une liste chaînée. Ensuite, vous parcourez tous les éléments; si vous trouvez que deux nœuds i et j sont connectés, vous allez regarder dans la liste chaînée de i. Si j n'est pas déjà là, vous l'ajoutez à la liste, et ajoutez également i à la liste de j. C'est plus simple si vous les ajoutez dans l'ordre.

Une fois que vous avez rempli votre liste de listes, vous connaissez maintenant le nombre d'entrées non nulles dans chaque ligne de la matrice: c'est la longueur de la liste de ce nœud. Ces informations sont exactement ce dont vous avez besoin pour préallouer une matrice clairsemée dans la structure de données de matrice de PETSc. Ensuite, vous pouvez libérer votre liste de listes car vous n'en avez plus besoin.

Cependant, cette approche suppose que tout ce que vous avez est la liste des nœuds que contient chaque élément.

Certains packages de génération de maillage - Triangle par exemple - peuvent générer non seulement une liste d'éléments et les nœuds qu'ils contiennent, mais également une liste de chaque arête de votre triangulation. Dans ce cas, vous ne risquez pas de surestimer le nombre d'entrées non nulles: pour les éléments linéaires par morceaux, chaque arête vous donne exactement 2 entrées de matrice de rigidité. Vous utilisez quadratique par morceaux, donc chaque bord compte pour 4 entrées, mais vous avez l'idée. Dans ce cas, vous pouvez trouver le nombre d'entrées non nulles par ligne avec un passage dans la liste des bords à l'aide d'un tableau ordinaire.

Avec cette approche, vous devez lire un gros fichier supplémentaire à partir du disque dur, ce qui pourrait être plus lent que d'utiliser la liste des éléments si votre calcul réel n'est pas si gros. Néanmoins, je pense que c'est plus simple.

Si vous spécifiez votre maillage en tant que DMPlex et votre disposition de données en tant que PetscSection, alors DMCreateMatrix () vous donnera automatiquement la matrice correctement préallouée. Voici des exemples PETSc pour le problème de Poisson et le problème de Stokes .

A voté

Personnellement, je ne connais aucun moyen bon marché de le faire, donc je surestime simplement le nombre, c'est-à-dire que j'utilise une valeur raisonnablement grande pour toutes les lignes.

Par exemple, pour un maillage parfaitement structuré composé d'éléments hexagonaux linéaires à 8 nœuds, les nnz par ligne dans les blocs diagonaux et hors diagonale sont dof * 27. Pour la plupart des maillages hexadécimaux générés automatiquement et non structurés, le nombre dépasse rarement le ddl * 54. Pour les tets linéaires, je n'ai jamais eu besoin d'aller au-delà de dof * 30. Pour certains maillages avec des éléments très mal formés / à faible rapport d'aspect, vous devrez peut-être utiliser des valeurs légèrement plus grandes.

La pénalité est que la consommation de mémoire locale (sur rang) se situe entre 2x-5x, vous devrez donc peut-être utiliser plus de nœuds de calcul sur votre cluster que d'habitude.

En fait, j'ai essayé d'utiliser des listes consultables, mais le temps nécessaire pour déterminer la structure de densité était plus important que l'assemblage / la résolution. Mais mon implémentation était très simple et n'utilisait pas d'informations sur les bords.

L'autre option consiste à utiliser des routines comme DMMeshCreateExodus comme indiqué dans cet exemple.

Vous cherchez à énumérer toutes les connexions uniques (gi, gj), ce qui suggère de les placer toutes dans un conteneur associatif (sans duplication), puis de compter sa cardinalité - en C ++, ce serait un std :: set <std :: pair <int, int>>. Dans votre pseudocode, vous remplaceriez "nnz [i] ++" par "s.insert [pair (gi, gj)]", puis le nombre final de nonzeros est s.size (). Il doit s'exécuter en temps O (n-log-n), où n est le nombre de non-zéros.

Puisque vous connaissez probablement déjà la gamme des gi possibles, vous pouvez "splay" la table par l'index gi pour améliorer les performances. Cela remplace votre ensemble par un std :: vector <std :: set <int>>. Vous remplissez cela avec "v [gi] .insert (gj)", alors le nombre total de nonzeros provient de la somme de v [gi] .size () pour tous les gi. Cela devrait fonctionner en temps O (n-log-k), où k est le nombre d'inconnues par élément (six pour vous - essentiellement une constante pour la plupart des codes pde, sauf si vous parlez de méthodes hp).

(Remarque - je voulais que ce soit un commentaire sur la réponse sélectionnée, mais c'était trop long - désolé!)

Partir d'une matrice clairsemée $E^T$ des éléments de taille$\times$dofs.