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

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 (Lagrange), une zone porte les multiplicateurs — les réactions s’y lisent directement.

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

Factorisation réutilisable (cache transparent)

solve factorise une fois, résout N fois : la factorisation est mise en cache dans la Matrix (état dérivé, non sérialisé, à mutabilité intérieure). La première résolution factorise et met en cache ; les suivantes sur la même matrice ne font que la descente/remontée — bien moins cher (cas de charge multiples, itérations de Newton, transitoire à matrice constante). Le cache est invalidé automatiquement dès que la matrice change (add_sub). On ne stocke jamais l’inverse explicite (dense, coûteux, instable) — seulement la factorisation.

Options de solve :

• method — méthode directe ("lu" par défaut ; l’énum laisse la place à d’autres back-ends sans changer les appels) ;
• cache — réutiliser/peupler le cache (True par défaut ; False factorise à neuf sans toucher le cache).

solve_eliminate(matrix, model, rhs) → NodeField

Voie alternative pour un modèle contraint : au lieu de border le système par des multiplicateurs de Lagrange (ce que fait solve sur la matrice augmentée), 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 LU creux (sur une matrice plus petite et définie, sans degré multiplicateur), puis prolongé u = T·û + u₀.

• lit la structure des contraintes via le seam méthode-neutre Constraint::relations() (partagé avec la voie Lagrange) ; le K physique est extrait du bloc point-selle assemblé (nœuds multiplicateurs filtrés) ;
• récupère en post-traitement la réaction (équivalent du multiplicateur), −(K·u − f) à la ligne duale de chaque esclave (= aₛ·λ) ;
• met en cache la condensation (T, K̂ factorisé) sur la matrice, comme la factorisation LU ; mêmes options method / cache ;
• un modèle sans contrainte retombe sur un solve simple.

Périmètre v1 : non chaîné, esclaves disjoints — chaque relation élimine un esclave distinct, jamais réutilisé comme maître ni esclave ailleurs (couvre la périodicité ; erreur explicite sinon).

K = pyrucast.matrix.stiffness(model, materials)
lagrange = pyrucast.solver.solve(K, rhs)  # système augmenté
condense = pyrucast.solver.solve_eliminate(K, model, rhs)  # reduced system — same field

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

solve_unilateral(matrix, model, rhs) → NodeField

Solveur actif/inactif (méthode du statut) pour un modèle portant des relations unilatérales (contraintes construites avec sense=">=" / "<=") : chaque relation est soit active (imposée en égalité, λ = la réaction), 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 λᵣ porte la réaction), 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é du modèle sont imposées inconditionnellement, comme par solve ; un modèle sans inégalité retombe sur un solve simple ;
• la structure des contraintes est lue via le seam méthode-neutre Constraint::relations() (partagé avec les voies Lagrange et élimination) ;
• options : method / cache (comme solve), active_set (stratégie de factorisation, ci-dessous), max_iter (borne de la boucle de statut, 100), tol (tolérance de signe sur λ et sur le jeu, 1e-10).

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_unilateral(K, model, rhs)  # "schur" by default
reaction = solution.value(mult_node, "lambda_T")  # 0 if the stop is released

# Forcing the old method (refactorization at every step):
sol2 = pyrucast.solver.solve_unilateral(K, model, 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_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 de réaction) 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.