Opérateurs de solveur

Le module ops::solver résout le système linéaire A · x = b issu de l’assemblage. Le back-end est un LU creux parallèle (faer), confiné à ops::solver : il consomme le CSC assemblé par nalgebra-sparse. Chaque crate garde son rôle — nalgebra (primitives), nalgebra-sparse (assemblage/stockage/serde), faer (factorise & résout). Voir Parallélisme.

solve(matrix, rhs, …) → NodeField

Un seul point d’entrée. Résout A x = b où A est la Matrice assemblée et b le NodeField de chargement. L’opération :

1. obtient la factorisation de A (cache de la matrice, cf. ci-dessous) ;
2. lit le NodeField de chargement à chacune des lignes de la matrice (les entrées absentes valent 0.0 par défaut) ;
3. effectue la descente/remontée pour ce second membre ;
4. emballe la solution dans un NodeField indexé sur les colonnes de la matrice (les primales : déplacements, températures, multiplicateurs…).

La solution a une zone par support colonne des blocs de la matrice, chaque zone vivant sur le handle POI1 du bloc lui-même (aucun support n’est reconstruit — ces supports sont créés une fois à la construction des sous-modèles et réutilisés d’assemblage en assemblage). Elle est donc same_support avec tout champ posé sur ces supports, et deux résolutions successives partagent les mêmes supports : leur arithmétique (a - b, …) s’aligne zone à zone. Sur un modèle contraint, une zone porte les multiplicateurs lambda, que les relations aient été gardées ou éliminées.

Pour poser un autre champ sur ces mêmes supports (p. ex. projeter les forces externes avant de calculer un résidu f_ext − K·u), la matrice expose ses supports en maillages : k.row_mesh() (côté dual, où vivent le second membre et mul_field) et k.col_mesh() (côté primal, où vit la solution) — handles partagés, dédupliqués. restrict(&f_ext, &k.row_mesh()?) s’aligne alors zone à zone avec K·u ; pour un résidu strict (toute composante soustraite, absente lue à 0 — et non passée brute par l’union), restrict_like(&f_ext, &f_int) reprojette aussi sur les composantes.

Une matrice singulière (p. ex. conditions aux limites oubliées) produit un pivot nul ⇒ solution non finie ⇒ erreur explicite.

K = pyrucast.matrix.stiffness(model, materials)
solution = pyrucast.solver.solve(K, rhs)  # factorise puis résout
T = solution.value(some_node, "T")

# Later solves on the SAME matrix: the factorization is reused.
sol2 = pyrucast.solver.solve(K, autre_rhs)  # descente/remontée seulement
sol3 = pyrucast.solver.solve(
    K, autre_rhs, cache=False
)  # refactorizes, without touching the cache

La matrice se suffit

solve ne prend pas de modèle : tout ce qu’il lui faut est dans la matrice.

Le sens voyage avec le bloc : une somme de matrices (K + K_σ + M) garde les blocs de ses opérandes, et l’archive les écrit. Rien ne peut donc contredire la matrice, et une matrice de contact ne peut pas être résolue comme si ses relations étaient des égalités.

Options

Les voies de la résolution

solve choisit son calcul, sa voie, d’après deux choses : l’option relations, et la présence d’inégalités dans la matrice (un bloc C dont le sens est >= ou <=). Dans la boucle de l’ensemble actif, l’option active_set choisit ensuite comment chaque statut est factorisé. Pour chaque voie, le schéma dit ce qui est calculé, et ce que la matrice garde ensuite dans sa case de cache (section suivante).

solve(matrice, second_membre, relations=…, method=…, active_set=…)
├── relations="lagrange" (défaut)
│   ├── la matrice ne porte aucune inégalité
│   │     calcul : factorisation de la matrice, LU ou Cholesky selon method
│   │     case   : la factorisation
│   └── la matrice porte des inégalités
│         calcul : boucle de l'ensemble actif sur la matrice, toujours en LU
│         ├── active_set="schur" (défaut)
│         │     calcul : une LU de la base (la matrice, toutes les inégalités
│         │              relâchées), faite une fois ; à chaque statut, une
│         │              petite matrice dense, une ligne par relation active
│         │     case   : l'état de l'ensemble actif : les inégalités, la LU de
│         │              la base, les colonnes déjà calculées, le dernier
│         │              statut convergé
│         │     base singulière : calcul de active_set="refactorize"
│         └── active_set="refactorize"
│               calcul : une LU de toute la matrice à chaque statut
│               case   : l'état de l'ensemble actif : les inégalités, le
│                        dernier statut convergé et sa LU
└── relations="eliminate"
    ├── la matrice ne porte aucune relation
    │     calcul : factorisation de la matrice, LU ou Cholesky selon method
    │     case   : la factorisation
    ├── la matrice ne porte que des égalités
    │     calcul : élimination des égalités, puis factorisation du système
    │              réduit K̂ = TᵀKT, LU ou Cholesky selon method
    │     case   : la condensation : T, Tᵀ, K sans ses contraintes, K̂ factorisé
    └── la matrice porte des inégalités
          calcul : élimination des égalités, puis boucle de l'ensemble actif
                   sur le système réduit [[K̂, Bᵀ], [B, 0]], toujours en LU
          ├── active_set="schur" (défaut)
          │     calcul : une LU de la base (le système réduit, toutes les
          │              inégalités relâchées), faite une fois ; à chaque
          │              statut, une petite matrice dense, une ligne par
          │              relation active
          │     case   : la condensation : T, Tᵀ, K sans ses contraintes, le
          │              système réduit, et l'état de l'ensemble actif : les
          │              inégalités, la LU de la base, les colonnes déjà
          │              calculées, le dernier statut convergé
          │     base singulière : calcul de active_set="refactorize"
          └── active_set="refactorize"
                calcul : une LU de tout le système réduit à chaque statut
                case   : la condensation : T, Tᵀ, K sans ses contraintes, le
                         système réduit, et l'état de l'ensemble actif : les
                         inégalités, le dernier statut convergé et sa LU

method ne s’applique qu’aux factorisations directes : la boucle de l’ensemble actif factorise en LU quoi qu’on demande, parce que le système d’un statut n’est jamais symétrique (une relation inactive y devient la ligne λ = 0, sa colonne restant en place).

Les relations sont mises à l’échelle de la structure. Par la voie des multiplicateurs, un système [[K, Cᵀ], [C, 0]] mêle deux ordres de grandeur : la raideur de la pièce (E·h, de l’ordre de 10⁹ pour de l’acier) et les coefficients des relations (1 pour un blocage). Factorisé tel quel, il perd de la précision justement sur les relations : les nœuds bloqués d’une pièce en acier gardaient quelques millionièmes du plus grand déplacement. Avant la factorisation, chaque relation est donc ramenée à l’échelle des inconnues qu’elle lie : sa ligne (imposed) et la colonne de son multiplicateur (lambda) sont multipliées par

\[ s_r = 1{,}5\, \frac{\max_k |K|_k}{\max_k |a_{rk}|}, \]

où \( k \) parcourt les inconnues de la relation, \( a_{rk} \) ses coefficients, et \( |K|_k \) la plus grande valeur, en valeur absolue, de la colonne de l’inconnue \( k \) hors des relations. Une inconnue sans rien à mesurer retombe sur la plus grande valeur de la matrice. Un facteur par relation, et non un seul pour toutes : une relation posée sur une pièce souple, à côté d’une pièce mille fois plus raide, est ramenée à l’échelle de sa pièce, non de l’autre. La solution est ensuite remise à l’échelle : c’est la même, et les relations tiennent à l’arrondi près. La boucle de l’ensemble actif garde ses tolérances dans les unités d’origine (un multiplicateur en newtons, un jeu en mètres). Cela coûte une copie des valeurs de la matrice, faite seulement quand elle porte des relations.

Une relation relâchée tient à l’arrondi près. Dans l’ensemble actif, la ligne d’une relation inactive devient \( d_r\,\lambda_r = 0 \), avec un pivot \( d_r \) de la taille des coefficients de la relation mis à l’échelle. Un pivot de 1 à côté de coefficients de \( 10^{12} \) laissait \( \lambda_r \) à l’arrondi de la résolution, quelques newtons, qui poussaient la pièce par \( C^{\mathsf T} \). Le \( \lambda_r \) rendu vaut ensuite exactement 0.

Factorisation réutilisable (cache transparent)

solve factorise une fois, résout N fois. Ce qui coûte est rangé dans la Matrix, dans une seule case : selon la voie, une factorisation, une condensation ou l’état de l’ensemble actif (schéma ci-dessus). La condensation et l’état de l’ensemble actif contiennent eux-mêmes leurs factorisations. Les résolutions suivantes sur la même matrice reprennent ce contenu et ne font que la descente/remontée, bien moins chère (cas de charge multiples, itérations de Newton, transitoire à matrice constante) ; la boucle de l’ensemble actif repart en plus du dernier statut convergé.

• Le contenu n’est repris que s’il répond à la demande : la même voie, et pour une factorisation, seule ou dans une condensation, la même method. Sinon il est recalculé et remplace l’ancien : alterner deux voies sur une même matrice recalcule à chaque changement.
• Il est oublié dès que la matrice change (add_sub), et n’est pas sauvegardé. cache=False ne le lit ni ne le remplit.
• On ne stocke jamais l’inverse explicite (dense, coûteux, instable) — seulement la factorisation.

L’élimination : relations="eliminate"

Voie alternative pour un modèle contraint : au lieu de border le système par des multiplicateurs de Lagrange, on élimine les contraintes par condensation maître/esclave. Pour chaque relation Σ aₖ·u(nœudₖ, varₖ) = g, un terme esclave s est exprimé par les autres (maîtres) : u_s = (g − Σ_{k≠s} aₖ·u_k)/a_s. Le système se réduit à K̂ û = f̂ avec K̂ = Tᵀ K T, résolu par le même back-end (sur une matrice plus petite et définie, sans degré multiplicateur), puis prolongé u = T·û + u₀.

• les relations sont lues dans la matrice, comme ci-dessus ; le K physique est son bloc point-selle sans les nœuds multiplicateurs ;
• récupère en post-traitement le multiplicateur de chaque relation, aₛ·λ = f − K·u à la ligne de son esclave, et le rend à son nœud : la solution est celle de la voie Lagrange, mêmes inconnues sur les mêmes supports ;
• met en cache la condensation (T, K̂ factorisé) sur la matrice, comme la factorisation ;
• une matrice sans relation retombe sur la factorisation simple.

Les inégalités ne s’éliminent pas. Qu’une relation unilatérale tienne en égalité dépend de la solution : seules les égalités sont éliminées, et chaque inégalité C_I·u ≥ g_I reste une relation du système réduit, écrite sur les inconnues libres, avec son multiplicateur :

⎡ Tᵀ K T   (C_I T)ᵀ ⎤ ⎡û⎤   ⎡Tᵀ(f − K·u₀)⎤
⎣ C_I T       0     ⎦ ⎣λ⎦ = ⎣ g_I − C_I·u₀ ⎦

Les multiplicateurs des égalités ont disparu (C_E·T = 0) ; ceux des inégalités restent, et la même boucle d’ensemble actif décide lesquelles sont saturées. Le λ de chaque inégalité revient à son nœud multiplicateur, comme par la voie Lagrange ; celui de chaque égalité éliminée est la part propre de son esclave, celle des inégalités (C_Iᵀλ) ôtée.

Périmètre : non chaîné, esclaves disjoints — chaque égalité élimine un esclave distinct, jamais réutilisé comme maître ni esclave ailleurs (couvre la périodicité ; erreur explicite sinon). Une inégalité peut porter sur un esclave : elle passe par T comme toute autre ligne.

K = pyrucast.matrix.stiffness(model, materials)
lagrange = pyrucast.solver.solve(K, rhs)  # système augmenté
condense = pyrucast.solver.solve(K, rhs, relations="eliminate")  # reduced system

Voir l’exemple examples/mpc_condensation.py et la page Contraintes.

Les relations unilatérales : l’ensemble actif

Dès que la matrice porte une relation unilatérale (contrainte construite avec sense=">=" / "<=", ou un contact), solve la résout par la méthode actif/inactif (méthode du statut) : chaque relation est soit active (imposée en égalité, son multiplicateur λ du signe d’une poussée), soit inactive (λ = 0, le jeu reste du côté admissible).

Conditions KKT et boucle de statut

Une relation unilatérale Cᵣ·u ≥ gᵣ (ou ≤) n’obéit pas à une équation mais aux conditions de complémentarité de Karush–Kuhn–Tucker : soit elle est active (Cᵣ·u = gᵣ, le multiplicateur λᵣ non nul), soit inactive (le jeu Cᵣ·u − gᵣ est du côté admissible et λᵣ = 0). Les deux ne peuvent être violées à la fois. On ne sait pas a priori quelles relations sont actives ; la méthode du statut itère dessus :

1. partir d’un statut d’essai (toutes actives, ou le statut convergé précédent quand le cache est chaud — un warm start) ;
2. résoudre le système point-selle avec, pour chaque relation inactive, sa ligne de contrainte remplacée par λᵣ = 0 (la matrice garde sa taille) ;
3. vérifier les signes : une relation active dont le λ tire (signe inadmissible pour son sens) est relâchée ; une relation inactive dont le jeu pénètre est activée ;
4. aucun changement de statut ⇒ convergence (la boucle finie classique du statut) ; sinon on recommence.

La convention de signe vient du système point-selle assemblé K·u + Cᵀ·λ = f, C·u = g : contre le multiplicateur KKT μ ≥ 0 d’une contrainte ≥ on a λ = −μ. Donc, dans le champ solution : ≥ active a λ ≤ 0 (relâchée si λ > tol) ; ≤ active a λ ≥ 0 (relâchée si λ < −tol).

Les relations d’égalité sont imposées inconditionnellement, comme sans inégalité.

Deux stratégies de factorisation (active_set)

Les deux stratégies parcourent exactement la même trajectoire de statuts (mêmes tests KKT, même résultat convergé) — elles ne diffèrent que par la façon de factoriser le système d’un statut donné.

"refactorize" — refactorisation par itération. La méthode d’origine : à chaque changement de statut, on refactorise le point-selle creux complet (une LU faer par itération). Robuste (aucune hypothèse sur la structure), mais paie une factorisation creuse à chaque pas.

"schur" (défaut) — complément de Schur / opérateur de Delassus. On factorise une seule fois le socle sans inégalités A (physique K + contraintes d’égalité, toutes les relations unilatérales relâchées), on le met en cache sur la matrice, et on obtient chaque statut par une mise à jour dense.

Le point clé : passer une relation r de l’état relâché à l’état actif ne change qu’une seule ligne de A. Dans le socle, la ligne relâchée porte l’identité λᵣ = 0 (un 1 à la colonne du multiplicateur, notée λcolᵣ) ; l’activer y restaure la vraie ligne de contrainte Cᵣ. Restaurer les k relations actives est donc une mise à jour de rang k :

M = A + Σ_{r actif} e_row(r) · (Cᵣ − e_λcol(r))ᵀ
  = A + U Vᵀ,   U = [e_row(r)],   Vᵣ = Cᵣ − e_λcol(r)

La formule de Sherman–Morrison–Woodbury donne alors la solution du statut sans refactoriser A :

x = A⁻¹·b − X · (I + Vᵀ X)⁻¹ · (Vᵀ A⁻¹ b),   X = [A⁻¹·e_row(r)]

• les colonnes xᵣ = A⁻¹·e_row(r) (une descente/remontée creuse par relation) sont mises en cache paresseusement : calculées la première fois qu’une relation devient active, réutilisées ensuite ;
• le petit système k × k G = I + Vᵀ X est l’opérateur de Delassus restreint aux relations actives — dense, factorisé par une LU dense (nalgebra) à chaque itération (coût k³/3, négligeable jusqu’à quelques milliers de contacts) ; ses entrées se lisent des colonnes cachées : Gᵢⱼ = δᵢⱼ + Cᵢ·xⱼ − xⱼ[λcolᵢ] ;
• une itération de statut ne coûte donc plus aucune factorisation creuse — seulement des descentes/remontées sur A (cachée) et une LU dense k × k. Un re-solve à chargement identique ou proche est quasi gratuit.

Repli automatique sur socle singulier. Le socle A doit être inversible, c.-à-d. la structure doit tenir sans aucun contact (bloquée par ailleurs). Un corps simplement posé sur un appui n’a pas ce luxe : A est singulière (mode rigide) alors que la méthode du statut converge très bien. Comme la LU creuse peut factoriser une matrice singulière en valeurs finies fausses sans erreur, la non-singularité du socle est confirmée par un aller-retour A⁻¹·(A·1) ≈ 1 ; s’il échoue, la voie "schur" retombe automatiquement sur "refactorize" (marqué une fois pour toutes sur la matrice). Aucune régression possible : le résultat est le même, seul le coût change.

K = pyrucast.matrix.stiffness(model, materials)  # model with sense=">="
solution = pyrucast.solver.solve(K, rhs)  # the matrix carries the ">="
reaction = solution.value(mult_node, "lambda")  # 0 if the stop is released

# Forcing the old method (refactorization at every step):
sol2 = pyrucast.solver.solve(K, rhs, active_set="refactorize")

Voir la section « Relations unilatérales » de la page Contraintes pour les conditions de complémentarité et la convention de signe.

Calculs plus gros que la RAM

Sur un gros modèle, c’est la factorisation qui sature la mémoire, pas l’assemblage. Mesuré sur un cube de conduction thermique HEX8, face inférieure imposée :

Premier levier, donc : quand le système éliminé est symétrique défini positif (thermique, élasticité), solve(…, relations="eliminate", method="cholesky") divise le pic par cinq et le temps par huit.

Un cube est le pire cas du remplissage. Une pièce mince ou élancée s’en tient bien en dessous, à nombre de DDL égal.

Déborder sur disque sans être root

Une roue compilée avec la feature spill (Linux ; c’est le cas de la roue Python) sait placer ses grosses allocations dans des fichiers mappés en mémoire plutôt qu’en mémoire anonyme. Le noyau peut écrire ces pages sur disque et les évincer sous pression, puis les recharger au besoin. C’est un swap, sans swap configuré et sans droit root. Les facteurs de faer en profitent sans le savoir.

Les variables sont lues une seule fois, à la première allocation du processus. Il faut donc les poser avant de le lancer, par exemple PYRUCAST_SPILL_DIR=/scratch/moi python calcul.py, et non depuis le script. Un répertoire illisible arrête le processus avec un message. Les fichiers sont anonymes (O_TMPFILE) : rien à nettoyer, même après un plantage.

Le répertoire doit être sur un disque local. /tmp est souvent un tmpfs, c’est-à-dire de la RAM, et y déborder ne sert à rien. Un montage réseau (NFS) rendrait chaque éviction très lente. Il faut aussi la place des facteurs : l’espace est réservé à l’allocation, donc un disque plein fait échouer l’allocation au lieu de planter plus loin.

Même cube à 363k DDL, Cholesky, RAM suffisante, débordement sur ext4, les cinq configurations dans une même série :

La solution est identique au bit près dans les cinq cas.

Le test ne coûte rien de mesurable : sans répertoire de débordement, ou avec un seuil au-dessus du plus gros bloc, on retrouve le temps du binaire compilé sans la feature — l’écart entre les trois premières lignes est du bruit.

Le pic total ne bouge pas, et c’est voulu : tant que la RAM est libre, le noyau garde en cache les pages des fichiers. Ce qui change, c’est qu’elles sont évinçables. La colonne qui décide si un calcul passe ou se fait tuer est le pic anonyme, qui tombe ici de 6,84 Go à 1,13 Go.

Le prix est payé même sans pression mémoire : une page écrite d’un fichier mappé part sur disque au bout d’une trentaine de secondes, pression ou non, et ce délai n’est réglable que par root. D’où un débordement à la demande, par exécution. Ici, ×2,7 à ×2,9 sur un RAID à ~50 Mo/s utiles ; un NVMe en demanderait beaucoup moins.

Choisir le seuil

Le seuil ne connaît pas les types : tout tampon assez gros déborde, quel qu’il soit. Le monter haut ne laisse partir que les tableaux qui comptent vraiment, et garde en RAM les données petites et souvent relues. Sur le cube ci-dessus, avec un seuil de 1 Gio, deux blocs seulement sont partis sur disque — 5,22 Go, qui sont les valeurs du facteur de Cholesky, et 1,08 Go de tampon de factorisation ; la CSR (146 Mo) est restée en mémoire. Le seuil bas, lui, fait déborder vingt-cinq blocs pour 0,9 Go d’anonyme gagnés de plus.

Un seuil de l’ordre du gigaoctet est donc le réglage d’un gros calcul ; un seuil bas ne se justifie que si la RAM manque à ce point. L’écart de temps entre les deux, environ 8 %, est du même ordre que la variabilité d’une exécution à l’autre.

spill_stats() en Python, spill::stats() en Rust, rendent le seuil, le nombre de blocs débordés, le plus gros, ce qui est mappé à l’instant et le maximum d’un coup — en octets, le seuil valant None quand rien ne déborde. PYRUCAST_SPILL_LOG donne la même chose bloc par bloc, au fil de l’eau, dans l’unité qui se lit le mieux.

stats = pyrucast.spill_stats()
if stats["threshold"] is None:
    pass  # PYRUCAST_SPILL_DIR absente : rien ne déborde, tout est en RAM
else:
    print(f"{stats['count']} bloc(s), le plus gros {stats['largest'] / 2**30:.1f} Gio")
    print(f"au plus {stats['peak'] / 2**30:.1f} Gio mappés d'un coup")

Le journal, lui, ressemble à ceci :

pyrucast spill: +5.2 GB, 5.2 GB mapped
pyrucast spill: +1.0 GB, 6.2 GB mapped
pyrucast spill: -1.0 GB, 5.2 GB mapped

Déterminisme

Contrairement au reste des opérateurs (bit-à-bit identiques quel que soit le nombre de threads), le solveur n’est pas bit-à-bit identique à l’ancien LU dense : pivotage et ordering diffèrent. Les résultats restent dans les tolérances numériques usuelles.

Exemples complets

La résolution de bout en bout (assemblage + contraintes + lecture de la solution et des multiplicateurs) est déroulée sur des cas à solution analytique :

• Dirichlet — Poisson 1-D -u'' = 0, u(0)=0, u(1)=1 ;
• Conduction thermique — ligne chauffée et carré ;
• Mécanique — treillis, élasticité, poutres.