Ce commit est pour améliorer quelques routines liées au contact et pour mettre progressivement mes développements sur la version courante.
J'ai refait un petit nettoyage pour l'interface des méthodes de mise à jour de la matrice de raideur tangente.
L'idée principale de la méthode est de pouvoir conserver la matrice de raideur tangente au cours des itérations Newton-Raphson pour la résolution de l'équilibre. Ainsi, étant donné que l'assemblage de la matrice de raideur peut être relativement couteux dans certains cas de figure, nous pouvons gagner du temps CPU en conservant la même matrice de raideur. Cependant, la vitesse de convergence de la méthode de Newton-Raphson modifiée devient linéaire (voir superlinéaire sous certaines conditions) au lieu de quadratique lorsque nous mettons à jour la matrice de raideur à chaque itération. Le fait de conserver la matrice de raideur peut dès lors donner lieu à une décroissance relativement lente du résidu lors des itérations, voir même pire, donner lieu à une croissance du résidu lors des itérations. Dans ce dernier cas de figure, nous rejetons l'itération et nous repartons de l'itération précédente où nous recalculons la matrice de raideur tangente (ce sont ce que l'on appelle des itérations rejetées !).
Auparavant, nous pouvons choisir cette méthode via l'option setNbOfIterationsForStiffUpdate() et l'option setCpuDependency() pour les itérations mécaniques/thermiques. Cette méthode existe toujours sous la dénomination StiffnessUpdateCriterion. Cependant, maintenant ces options sont directement à définir dans l'objet StiffnessUpdateCriterion et non IterationManager. Par conséquent, j'ai dû adapter l'interface des différents cas tests qui utilisent cette méthode dans la batterie.
Après un petit refactoring, pour ne pas changer, j'ai scindé en deux les méthodes de Newton Raphson Modifiée : CpuModifiedNRStiffnessUpdateCriterion et PeriodicModifiedNRStiffnessUpdateCriterion. Par ailleur, la méthode par défaut pour la stratégie de mise à jour de la matrice de raideur tangente est la méthode de Newton Raphson NRStiffnessUpdateCriterion. Toutes ces méthodes dérivent de StiffnessUpdateCriterionBase et ont pour but de décider si nous devons mettre à jour la matrice de raideur tangente lors du calcul de la correction des positions.
La différence entre les deux méthodes précitées réside dans la détermination du paramètre residualDecreaseFactor et maxNbOfModifiedNRIterations. Le premier paramètre définit la valeur maximale autorisée de décroissance entre deux résidus successifs sans mise à jour de la matrice de raideur tangente. Si le résidu ne décroit pas suffisamment, nous décidons de mettre à jour la matrice de raideur tangente (valeur typique est comprise entre 0.1 et 0.9). Le second paramètre permet de définir le nombre maximale d'itérations successives faites avec la même matrice de raideur tangente. Même si le résidu décroit suffisamment au cours des itérations sans mise à jour de la matrice de raideur tangente, nous déclarons qu'au delà de x itérations successives (sur le step courant), nous devons mettre à jour la matrice de raideur tangente. (valeur typique est comprise entre 2 et 5). Le choix de ces deux paramètres sont loin d'être trivial et dépend clairement de l'application.
Dans le cas de la méthode PeriodicModifiedNRStiffnessUpdateCriterion, nous les choisissons tous les deux de manière indépendante. Dans le cas de la méthode CpuModifiedNRStiffnessUpdateCriterion, ces deux paramètres sont déterminés à partir du rapport du temps CPU Real (important pour les cas-test lancés en parallèle !) entre une itération avec mise à jour de la matrice de raideur tangente et une itération sans mise à jour de la matrice de raideur tangente. Ce rapport est borné et doit être une valeur comprise entre 2 et 9. Cette variante permet une automatisation du choix des paramètres précédents.
Finalement, lors de l'intégration temporelle, il existe d'autres situations qui nécessitent un traitement particulier dans le cas la méthode de Newton Raphson modifiée. Le premier cas de figure est la mise à jour des connectivités (Stage, Contact défo-défo, Rupture) pour toutes les itérations et le second cas de figure est le variation du pas de temps pour la première itération (première correction du résidu). Au lieu d'avoir un énorme if/else if et des valeurs booléennes (problème de synchronisation - remise à zéro), j'ai opté pour une structure avec des objets (“Strategy Pattern”). Ainsi, une étape du type beginStep dans l'intégration temporelle va informer le ModifiedNRStiffnessUpdateCriterion que la méthode implémentée pour mettre à jour la matrice de raideur doit être BeginStepStiffnessUpdateMethod. Ainsi, lors de l'étape de décision de mise à jour du pas de temps (qui n'est pas du tout au même moment, en réalité, c'est juste après l'évaluation du résidu), puisqu'il y a déjà une méthode sélectionnée, nous exécuterons cette méthode et une fois exécutée, on va détruire cette méthode. Si nous n'avons pas de méthode pré-établie, la méthode par défaut sera ResidualStiffnessUpdateMethod. A chaque évènement particulier dans l'intégration temporelle ou la résolution de l'équilibre, nous pouvons y associer une méthode de mise à jour de la matrice de raideur tangente. Celà reste de petit objet et dès lors leur destruction et leur construction n'affecte pas de manière significative le temps CPU total et cela permet une plus grande flexibilité pour traiter les cas particuliers.
Voici les méthodes disponibles et les options associées.
mim = metafor.getMechanicalIterationManager()
stiffUpdate = NRStiffnessUpdateCriterion(mim)
stiffUpdate.setDebug(parameters['StiffUpdateDebug']); Default = False
stiffUpdate.setVerbose(parameters['StiffUpdateVerbose']); Default = False
mim.setStiffnessUpdateCriterion(stiffUpdate)
mim = metafor.getMechanicalIterationManager()
stiffUpdate = StiffnessUpdateCriterion(mim)
stiffUpdate.setNbOfIterationsForStiffUpdate(parameters['NbOfIterationsForStiffUpdate']); Default = 1
stiffUpdate.setCpuDependency(parameters['CpuDependency']); Default = False
stiffUpdate.setDebug(parameters['StiffUpdateDebug']); Default = False
stiffUpdate.setVerbose(parameters['StiffUpdateVerbose']); Default = False
mim.setStiffnessUpdateCriterion(stiffUpdate)
mim = metafor.getMechanicalIterationManager()
stiffUpdate = CpuModifiedNRStiffnessUpdateCriterion(mim)
stiffUpdate.setKeepStepStiffness(parameters['KeepStepStiffness']); Default = True
stiffUpdate.setDebug(parameters['StiffUpdateDebug']); Default = False
stiffUpdate.setVerbose(parameters['StiffUpdateVerbose']); Default = False
mim.setStiffnessUpdateCriterion(stiffUpdate)
mim = metafor.getMechanicalIterationManager()
stiffUpdate = PeriodicModifiedNRStiffnessUpdateCriterion(mim)
stiffUpdate.setMaxNbOfModifiedNRIterations(parameters['MaxNbOfModifiedNRIterations']); Default = 1
stiffUpdate.setResidualDecreaseFactor(parameters['ResidualDecreaseFactor']); Default = 1.0
stiffUpdate.setKeepStepStiffness(parameters['KeepStepStiffness']); Default = True
stiffUpdate.setDebug(parameters['StiffUpdateDebug']); Default = False
stiffUpdate.setVerbose(parameters['StiffUpdateVerbose']); Default = False
mim.setStiffnessUpdateCriterion(stiffUpdate)
Le paramètre keepStepStiffness permet de conserver la matrice de raideur tangente d'un pas de temps à l'autre, uniquement si le pas de temps est quasi-identique. Sinon, on force la mise à jour de la matrice de raideur tangente.
J'ai ajouté une méthode pour calculer le pas de temps, qui se base à la fois sur le nombre d'augmentation et sur les nombre d'itérations en moyenne sur chaque résolution de l'équilibre suite à une augmentation. Elle permet ainsi l'estimation du nouveau pas de temps sur base du nombre d'augmentation si on fait peu d'augmentations ou du nombre d'itérations en moyenne sur chaque résolution de l'équilibre suite à une augmentation, si on fait beaucoup d'augmentations.
tscm = NbOfMechAugLagNRIterationsTimeStepComputationMethod(metafor) tscm.setEquivalentIteMethod(ALM_TSCM_AUGLAG) tsm.setTimeStepComputationMethod(tscm)
La raison principale qui m'a poussé à tenter de remettre à jour les méthodes de mise à jour de la matrice de raideur tangente est la combinaison de ce type de stratégie avec la méthode du Lagrangien augmenté. Lors de l'utilisation du Lagrangien augmenté pour traiter les contraintes de contact, nous avons en général une boucle supplémentaire liée aux augmentations, c'est à dire aux mises à jour des lagrangiens augmentés, jusqu'à satisfaire un critère de convergence sur la satisfaction des contraintes de contact. Cet algorithme permet de diminuer la valeur de la pénalité utilisée tout en gardant un contrôle sur la qualité des résultats du point de vue contact, ce qui permet de rendre la matrice de raideur tangente mieux conditionnée (en pratique), dès lors une meilleure convergence lors de la résolution du problème mécanique par la méthode Newton Raphson. Au fur et à mesure, que l'on réalise des augmentations, la résolution du problème mécanique se fait en général de manière plus aisée (moins d'itérations mécaniques), puisque les champs de contact sont de mieux en mieux connus au fur et à mesure des augmentations (meilleur connaissance de la condition aux limites). Par conséquent, on s'attend à ce que la matrice de raideur tangente ne varie presque plus d'une augmentation à l'autre.
Ainsi, il semble intéressant d'utiliser la méthode de Newton Raphson modifiée pour éviter de mettre à jour la matrice de raideur tangente à chaque itération mécanique. En pratique, la difficulté de résolution du problème mécanique réside dans les non-linéarités (contact, géométrie, matériau) présentes au pas de temps courant. Si le problème mécanique est quasi linéaire, le problème se résout en quasi une seule itération comme en élasticité linéaire. Il est possible d'envisager un basculement entre la méthode de Newton Raphson et la méthode de Newton Raphson modifiée au travers des augmentations. Lors du passage à la méthode de Newton Raphson modifiée, il est intéressant aussi de mettre à jour les paramètres ResidualDecreaseFactor et MaxNbOfModifiedNRIterations en fonction du nombre d'augmentations réalisés. On constate que la valeur du premier résidu d'une résolution d'équilibre mécanique décroit au fur et à mesure des augmentations (Attention, ceci n'est pas toujours cas, surtout si on a des status de contact changeant d'une augmentation à l'autre !), ce qui se traduit par une diminution du nombre d'itérations mécaniques pour atteindre la tolérance prescrite sur la procédure itérative de Newton Raphson. Il serait intéressant de faire une moyenne pondérée entre ce rapport entre les premiers résidus entre deux augmentations successives et celui entre deux itérations mécaniques. Ainsi on corrige progressivement les valeurs des deux paramètres et nous pouvons admettre plus de souplesse dans la résolution de l'équilibre mécanique par la méthode de Newton Raphson modifiée puisque le problème devient plus en plus “linéaire” au fur et à mesure des augmentations.
Voici les méthodes disponibles et les options associées.
alm = StandardAugmentedLagrangianManager(metafor)
mim = metafor.getMechanicalIterationManager()
stiffUpdate = AlmCpuModifiedNRStiffnessUpdateCriterion(mim,alm)
stiffUpdate.setAugmentationThreshold(parameters['AugmentationThreshold']); Default = 0
stiffUpdate.setKeepStepStiffness(parameters['KeepStepStiffness']); Default = True
stiffUpdate.setKeepAugmentationStiffness(parameters['KeepAugmentationStiffness']); Default = True
stiffUpdate.setDebug(parameters['StiffUpdateDebug']); Default = False
stiffUpdate.setVerbose(parameters['StiffUpdateVerbose']); Default = False
mim.setStiffnessUpdateCriterion(stiffUpdate)
alm = StandardAugmentedLagrangianManager(metafor)
mim = metafor.getMechanicalIterationManager()
stiffUpdate = AlmPeriodicModifiedNRStiffnessUpdateCriterion(mim,alm)
stiffUpdate.setMaxNbOfModifiedNRIterations(parameters['MaxNbOfModifiedNRIterations']); Default = 1
stiffUpdate.setResidualDecreaseFactor(parameters['ResidualDecreaseFactor']); Default = 1.0
stiffUpdate.setAugmentationThreshold(parameters['AugmentationThreshold']); Default = 0
stiffUpdate.setKeepStepStiffness(parameters['KeepStepStiffness']); Default = True
stiffUpdate.setKeepAugmentationStiffness(parameters['KeepAugmentationStiffness']); Default = True
stiffUpdate.setDebug(parameters['StiffUpdateDebug']); Default = False
stiffUpdate.setVerbose(parameters['StiffUpdateVerbose']); Default = False
mim.setStiffnessUpdateCriterion(stiffUpdate)
Le paramètre keepAugmentationStiffness permet de conserver la matrice de raideur tangente après une augmentation, uniquement si le premier résidu avant l'augmentation est plus grand que le premier résidu après l'augmentation. Sinon, on force la mise à jour de la matrice de raideur tangente.
Le paramètre augmentationThreshold permet de spécifier à partir de quel numéro d'augmentation, nous passons d'une méthode du type Newton Raphson à une méthode du type Newton Raphson modifiée. Si cette valeur est prise à -1, nous commençons d'office par une méthode du type Newton Raphson modifiée avant la première augmentation.
Finalement, l'algorithme est suffisamment malin pour passer à une méthode du type Newton Raphson modifiée si nous ne sommes pas en contact (avant et après impact, après springback), quelque soit la valeur de augmentationThreshold.
J'ai ajouté la contribution de la variation de l'aire nodale pour le calcul de la matrice de raideur tangente en 3D dans le cas du contact déformable-déformable.
[a]:mtFEM/algos/AlmCpuModifiedNRStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/AlmPeriodicModifiedNRStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/AlmStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/CpuModifiedNRStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/ModifiedNRStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/NRStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/PeriodicModifiedNRStiffnessUpdateCriterion.cpp [a]:mtFEM/algos/StiffnessUpdateCriterionBase.cpp [a]:mtFEM/algos/StiffnessUpdateMethods.cpp [a]:mtFEM/algos/AlmCpuModifiedNRStiffnessUpdateCriterion.h [a]:mtFEM/algos/AlmPeriodicModifiedNRStiffnessUpdateCriterion.h [a]:mtFEM/algos/AlmStiffnessUpdateCriterion.h [a]:mtFEM/algos/CpuModifiedNRStiffnessUpdateCriterion.h [a]:mtFEM/algos/ModifiedNRStiffnessUpdateCriterion.h [a]:mtFEM/algos/NRStiffnessUpdateCriterion.h [a]:mtFEM/algos/PeriodicModifiedNRStiffnessUpdateCriterion.h [a]:mtFEM/algos/StiffnessUpdateCriterionBase.h [a]:mtFEM/algos/StiffnessUpdateMethods.h [r]:
[a]:apps/bImp/cyl3DDynAnaModifiedNR.py [a]:apps/qs/cont2AugLagModifiedNR1.py [a]:apps/qs/cont2AugLagModifiedNR2.py [a]:apps/qs/cont2AugLagModifiedNR3.py [a]:apps/qs/cont2AugLagModifiedNR4.py [a]:apps/qs/cont2AugLagModifiedNR5.py [a]:apps/qs/cont2ModifiedNR.py [a]:mtContact/tests/twoToriContactTestAugLag.py [a]:mtContact/tests/twoToriContactTestAugLagModifiedNR.py [a]:mtContact/tests/twoToriContactTestModifiedNR.py [r]:
— gaëtan 2017/04/26 22:37