Commit 2016-04-25

Ce commit porte principalement sur l'amélioration du Lagrangien Augmenté.

Je vous met à profit les derniers développements sur l'algorithme du Lagrangien Augmenté.

J'ai constaté à ma grande surprise qu'il faut forcer le calcul d'une itération mécanique après une augmentation. Après chaque augmentation (modification des lagrangiens augmentés à configuration cinématique fixée), on résout l'équilibre mécanique par la méthode de Newton Raphson. Cependant, il se peut que la modification du vecteur des forces externes ne soit pas suffisamment conséquente pour changer la norme moyenne du résidu (alors que la norme max du résidu est nettement plus sensible à ces changements là …) et dès lors forcer une nouvelle itération mécanique, c'est à dire une correction des positions et obtenir une nouvelle configuration cinématique. Donc, l'algorithme continue de modifier les lagrangiens augmentés à configuration cinématique donnée jusqu'à espérer avoir une variation suffisante des forces externes pour avoir une itération mécanique. Si la configuration donnée ne change pas, les gaps normal et tangentiel ne peuvent pas évoluer et nous ne pouvons pas converger vers les tolérances imposées à l'algorithme du lagrangien augmenté. En résumé, l'algorithme stagne et l'intégration temporelle risque de s'arrêter tout simplement !

Il est vrai que les articles de référence en la matière travaillent en général en terme de norme énergétique (ce qui nécessite d'office le calcul d'une itération mécanique !).

Dans le jeux de données Python, on a

alm = AutomaticAugmentedLagrangianManager(metafor)
alm.setForceFirstIterationAfterAugmentation(True)
alm = StandardAugmentedLagrangianManager(metafor)
alm.setForceFirstIterationAfterAugmentation(True)

Par défaut, ce paramètre est mis à False !

J'ai développé une autre approche pour mettre à jour les lagrangiens augmentés sur base de l'histoire de l'évolution de la pression de contact en fonction du gap normal et le cisaillement de contact en fonction du gap tangentiel. Ces lagrangiens augmentés “extrapolés” sont mis à jour si besoin via le critère de frottement de Coulomb ou de Tresca. Au lieu de mettre à jour les coefficients de pénalité pour chaque nœud de contact à chaque augmentation ou au cours de chaque itération mécanique, on utilise cette autre approche développée par Wriggers and Zavarise en contact sans frottement qui possède un ordre de convergence super-linéaire, mais qui n'est plus spécialement monotone en convergence sans garde fou !

Avec ces derniers gardes fous, l'algorithme semble être encore plus efficace que l'approche standard car on tient de l'historique au cours des augmentations !

augLagAugmentation = AugLagExtrapolationAugmentation(alm)
augLagAugmentation.setUseStabilityCriterion(True)
augLagAugmentation.setUseExtrapolationCriterion(True)
augLagAugmentation.setExtrapolationFactor(_factor)
augLagAugmentation.setResetPreviousAugmentedLagrangians(True)

La dernière option permet d'activer la procédure dès la première augmentation en indiquant False, sinon il faut au moins une augmentation selon la procédure d'Uzawa avant de lancer celle par extrapolation.

Par défaut, tous les critères sont activés et le facteur d'extrapolation est 0.1 (c'est à dire on autorise une variation de maximum 10 % de la valeur courante de la pression de contact et du cisaillement de contact.).

Cette procédure est activée si on a le nœud de contact est en contact au cours de deux résolutions successives de l'équilibre mécanique (Le lagrangien normal est mis à jour en fonction de la valeur à l'augmentation k du couple (contact pressure, normal gap) et l'augmentation k-1 du même couple), sinon les lagrangiens augmentés sont mis à jour selon la méthode d'Uzawa.

Finalement, il est possible de mettre à zéro les lagrangiens augmentés d'un pas de temps à l'autre (très déconseillés !), quelque soit la méthode d'augmentation :

augLagAugmentation.setResetAugmentedLagrangians(True)

J'ai ajouté une nouvelle méthode de calcul de longueur de contact utilisée lors des critères géométriques normalisés. Dans les estimateurs de pas de temps critique dans les schémas explicites, la notion de dimension équivalente d'un élément apparaît naturellement dans la définition du pas de temps critique. Cependant, ce choix reste vague (distance minimale entre deux nœuds de la maille ou autre) et dès lors, il propose une analyse des valeurs propres à partir de laquelle il dérive une longueur équivalente (plus conservative !). Leur raisonnement est indépendant à priori du type de maille et de la dimension de l'espace. C'est pour ça, que j'ai repris leur développement et je l'ai codé dans Metafor pour le cas 2D, 3D et axisymmétrique.

En résumé, on a les méthodes suivantes :

• _method = ALM_AIC_LENGTH (Par défaut) : on calcule une longueur équivalente basée sur la surface de contact (en 2D, c'est la moyenne des demi-longueurs des courbes adjacentes au nœud et en 3D, c'est la moyenne de la racine carrée des quart aires des faces adjacentes au noeud.). On normalise par la même valeur pour les gaps normal et tangentiel.
• _method = ALM_CRITICAL_LENGTH : on calcule une longueur équivalente basée sur l'estimateur du pas de temps critique en explicit des éléments volumiques adjacents au nœud de contact. On normalise par la même valeur pour les gaps normal et tangentiel.
• _method = ALM_GEO_LENGTH : Pour le gap normal, on calcule une distance normale des mailles volumiques par rapport au plan tangent de contact ce qui donne une épaisseur des éléments volumiques sous-jacents. On divise cette distance par deux pour obtenir la longueur caractéristique. Pour le gap tangentiel, on procède de la même manière que dans la première méthode mais on travaille dans le plan tangent de contact ! Dès lors, ce critère est par définition anisotrope.

Par exemple, pour un élément quadrangulaire rectangulaire $a \times b$ (avec $b > a > 0$), on a

• _method = ALM_AIC_LENGTH : $Lc_N = Lc_T = \frac{b}{2}$
• _method = ALM_GEO_LENGTH : $Lc_N = \frac{a}{2}$ et $Lc_T = \frac{b}{2}$
• _method = ALM_CRITICAL_LENGTH : $Lc_N = Lc_T = \frac{a}{\sqrt{1+\frac{b^2}{a^2}}}$ (= Moitié de la diagonale si on a un carré $b = a$.)

Dans le jeux de données Python, on choisit les options comme ceci :

augLagGeoCriterion = AugLagNormalisedGeoCriterion(alm)
augLagGeoCriterion.setCharacteristicLengthMeth(_method)

Sur un pas de temps donné, la maille peut subir de grande déformation d'où le fait que l'on recalcule à chaque fois la longueur caractéristique au cours des augmentations.

L'algorithme du Lagrangien augmenté permet à priori de résoudre un problème de contact sur un pas de temps plus élevé que l'équivalent en Méthode de Pénalité. Pour mettre à jour le pas de temps, j'ai proposé plusieurs méthodes :

• _equivalentIteMethod = ALM_TSCM_FIRST : Nombre d'itérations mécaniques avant la première augmentation (DEFAULT)
• _equivalentIteMethod = ALM_TSCM_MAX : Nombre d'itérations mécaniques maximales au cours des augmentations
• _equivalentIteMethod = ALM_TSCM_MEAN : Nombre d'itérations mécaniques moyennes au cours des augmentations
• _equivalentIteMethod = ALM_TSCM_WEIGHTEDMEAN : Nombre d'itérations mécaniques moyennes pondérées au cours des augmentations (Les poids sont le rapport entre le nombre d'itération mécanique divisée par le nombre d'itération mécanique maximale.)

Cette valeur d'itération “équivalente” est insérée dans l'algorithme de mise à jour du pas de temps proposé par JPP dans sa thèse.

Dans le jeux de données Python, on a

tscm = NbOfAugLagMechNRIterationsTimeStepComputationMethod(metafor)
tscm.setTimeStepDivisionFactor(_timeStepDivisionFactor)
tscm.setEquivalentIteMethod(_equivalentIteMethod)
tscm.setMaxNbOfIteInAdaptOfTimeStep(_maxNbOfIteInAdaptOfTimeStep)
tscm.setNbOptiIte(_nbOptiIte)
tsm.setTimeStepComputationMethod(tscm)

Par exemple, si on a le profil suivant de convergence sur un pas de temps de donné (12 augmentations pour 30 itérations mécaniques) :

on obtient l'itération mécanique équivalente suivante :

Il est possible de demander d'afficher des informations sur le status des nœuds de contact d'une interaction.

ci = RdContactInteraction(1)
ci.setShowInformation(True)

Voici un exemple de sortie.

mechanical iteration 0  : Rmean = 1.000000e+000  Rmax = 2.244123e+002 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 4 / 11
Number of Nodes in Sticking/Sliding Contact = 4 / 0
mechanical iteration 1  : Rmean = 2.192696e-001  Rmax = 4.097902e+001 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 4 / 11
Number of Nodes in Sticking/Sliding Contact = 4 / 0
mechanical iteration 2  : Rmean = 2.016662e-001  Rmax = 2.821996e+001 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 3  : Rmean = 1.084675e-001  Rmax = 9.532517e+000 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 4  : Rmean = 2.831632e-002  Rmax = 2.978541e+000 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 5  : Rmean = 1.060746e-003  Rmax = 1.218047e-001 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0
mechanical iteration 6  : Rmean = 1.744038e-005  Rmax = 1.792596e-003 
RdContactInteraction #2
Actual/Total Number of Nodes in Contact = 3 / 11
Number of Nodes in Sticking/Sliding Contact = 3 / 0

Cela peut permettre de comprendre la raison pour laquelle un pas de temps ne converge pas (Status collant glissant qui change ou status en contact ou hors contact qui change).

Il est désormais possible d'estimer le nombre de conditionnement réciproque d'une matrice carrée (Matrix de Metafor) quelconque via les routines LAPACK. Les matrices utilisées pour la résolution du système linéaire sous-jacente au calcul de l'équilibre implicite se trouvent sous un forme creux liée au solveur utilisé. Cependant, il n'existe pas (sauf erreur de ma part) d'estimateur du nombre de conditionnement réciproque dans les routines LAPACK, donc nous devons transformer une matrice creuse sous un format plein pour pouvoir effectuer le calcul proprement dit (ce n'est pas top du tout mais bon faute de mieux …). J'ai réuni un peu de documentation sur ces estimateurs, si qqn veut avoir les infos.

On peut demander de calculer à chaque résolution du système linéaire, le nombre de conditionnement réciproque. Pour se faire, il suffit de faire :

solvermanager = metafor.getSolverManager()
solver = DSSolver()
solver.setShowRCond(True)
solvermanager.setSolver(solver)        

Également, il est possible de faire un extracteur de cette grandeur (Nombre de conditionnement réciproque de la matrice FreeJacobian du schéma d'intégration temporelle) via les commandes suivantes :

valuesmanager.add(1, MiscValueExtractor(metafor, EXT_RCOND), 'RCond')

Dans Matlab, on compare ce nombre de conditionnement réciproque avec l'epsilon machine multipliée par la taille du système à résoudre. Pour rappel, ce nombre est proche de 1 lorsque le système est bien conditionné et ce nombre est proche de 0 lorsque le système est mal conditionné.

Finalement, cet extracteur peut aider à quantifier l'influence du coefficient de pénalité sur la résolution du système linéaire.

Par la même occasion, j'ai ajouté plein de nouveau extracteurs de variables de contact. J'ai ajouté directement les infos dans la partie sur les extracteurs de contact. Finalement, la documentation est à jour à ce niveau là.

[a]:mtFEMBase/AugLagTimeStepComputationMethod.cpp
[a]:mtFEMBase/AugLagTimeStepComputationMethod.h
[a]:mtContact/src/AugLagExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagExtrapolationAugmentation.h
[a]:mtContact/src/AugLagNormalExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagNormalExtrapolationAugmentation.h
[a]:mtContact/src/AugLagTangentialExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagTangentialExtrapolationAugmentation.h
[a]:mtContact/src/AugLagFrictionalExtrapolationAugmentation.cpp
[a]:mtContact/src/AugLagFrictionalExtrapolationAugmentation.h
[a]:mtContact/src/NormalAugmentedLagrangianValueExtractor.cpp
[a]:mtContact/src/NormalAugmentedLagrangianValueExtractor.h
[r]:mtContact/src/GapValueExtractor.cpp => [a]:mtContact/src/NormalGapValueExtractor.cpp
[r]:mtContact/src/GapValueExtractor.h => [a]:mtContact/src/NormalGapValueExtractor.h
[a]:mtContact/src/Shear1ContactValueExtractor.cpp
[a]:mtContact/src/Shear1ContactValueExtractor.h
[a]:mtContact/src/Shear2ContactValueExtractor.cpp
[a]:mtContact/src/Shear2ContactValueExtractor.h
[a]:mtContact/src/Tangent1ForceValueExtractor.cpp
[a]:mtContact/src/Tangent1ForceValueExtractor.h
[a]:mtContact/src/Tangent2ForceValueExtractor.cpp
[a]:mtContact/src/Tangent2ForceValueExtractor.h
[a]:mtContact/src/TangentAugmentedLagrangian1ValueExtractor.cpp
[a]:mtContact/src/TangentAugmentedLagrangian1ValueExtractor.h
[a]:mtContact/src/TangentAugmentedLagrangian2ValueExtractor.cpp
[a]:mtContact/src/TangentAugmentedLagrangian2ValueExtractor.h
[a]:mtContact/src/TangentAugmentedLagrangianValueExtractor.cpp
[a]:mtContact/src/TangentAugmentedLagrangianValueExtractor.h
[a]:mtContact/src/TangentGap1ValueExtractor.cpp
[a]:mtContact/src/TangentGap1ValueExtractor.h
[a]:mtContact/src/TangentGap2ValueExtractor.cpp
[a]:mtContact/src/TangentGap2ValueExtractor.h
[a]:mtContact/src/TangentGapValueExtractor.cpp
[a]:mtContact/src/TangentGapValueExtractor.h

[r]:

[r]:
[a]:apps/qs/cont2AugLag.py

— Gaëtan WAUTELET 2016/04/25