5. Réglages des solveurs non-linéaires#

5.1. Choix de la méthode#

◇ METHODE = / "IMPLEX",
            / "MODELE_REDUIT",
            / "NEWTON" (par défaut),
            / "NEWTON_KRYLOV",

Il existe plusieurs catégories de solveurs non-linéaires accessibles: la méthode de Newton (et ses variantes), la méthode de Newton-Krylov, la méthode IMPLEX, l’accès au solveur non-linéaire SNES de la librairie PETSc et l’utilisation de méthodes de réduction de modèle (ROM). Le solveur est sélectionné par le mot-clef METHODE.

5.2. Méthode de Newton#

⟐ NEWTON = _F(
      ◇ REAC_INCR = int (défaut: 1),
      ◇ PREDICTION = / "DEPL_CALCULE",
                     / "ELASTIQUE",
                     / "EXTRAPOLE",
                     / "TANGENTE",
      ◇ MATRICE = / "ELASTIQUE",
                  / "TANGENTE" (par défaut),
      ◇ PAS_MINI_ELAS = float,
      ◇ REAC_ITER = int (défaut: 1),
      ◇ REAC_ITER_ELAS = int (défaut: 0),
      ◇ EVOL_NOLI = evol_noli,
      ◇ MATR_RIGI_SYME = / "NON" (par défaut),
                        / "OUI",
  ),

Quand METHODE= »NEWTON », on utilise l’algorithme de Newton-Raphson pour résoudre le problème (voir [R5.03.01]).

C’est le mot-clef facteur NEWTON (une seule occurrence) qui permet de régler les différents paramètres. Le mot-clef MATR_RIGI_SYME sert à forcer la symétrie de la matrice de rigidité. Il est utile lorsqu’on utilise simultanément une fonctionnalité utilisant une matrice non-symétrique (formalisme de Simo-Miehe, chargements de type pression suiveuse, etc) et la formulation de contact de type DISCRETE ou LIAISON_UNIL (cf. [U4.44.11]). En effet, ces dernières ne peuvent fonctionner qu’avec une matrice de rigidité symétrique. Le fait de rendre la matrice de rigidité symétrique a des conséquences sur la vitesse de convergence, pouvant aller jusqu’à son échec mais ne produira jamais de résultats faux. Si le fait de symétriser la matrice de rigidité rend impossible la convergence de l’algorithme de Newton, l’utilisateur est invité à utiliser une formulation de contact de type CONTINUE, qui fonctionne avec une matrice de rigidité non-symétrique.

5.2.1. Phase de prédiction#

La phase de prédiction (Cf. [R5.03.01]) a pour but de calculer une estimation du champ de déplacements afin de permettre à la méthode de Newton de converger plus rapidement. Losque le mot clé est absent, c’est la matrice tangente en vitesse (option RIGI_MECA_TANG) qui est utilisée si l’on a choisi pour la méthode de Newton une PREDICTION=”TANGENTE”, et c’est la matrice élastique (option RIGI_MECA) qui est utilisée si on a choisi PREDICTION=”ELASTIQUE”.

Il est également possible de faire une prédiction en utilisant des résultats préliminaires issus d’un autre calcul.

Pour PREDICTION=”EXTRAPOLE”, on calcule l’estimation de l’incrément de déplacement à partir de l’incrément total obtenu comme solution au pas de temps précédent (pondéré par le rapport des pas de temps). On projette cette estimation sur l’ensemble des champs cinématiquement admissibles (i.e. satisfaisant les conditions aux limites de Dirichlet) selon la norme donnée par la matrice élastique, qui doit donc être calculée. Cette fonctionnalité est intéressante dans le cas de l’utilisation de schémas d’intégration locale explicite de type Runge-Kutta qui ne fournissent pas de matrice tangente: dans ce cas la méthode de Newton utilise une matrice élastique, mais le nombre d’itérations nécessaires peut être élevé. L’utilisation de l’extrapolation peut améliorer les performances.

Pour PREDICTION=”DEPL_CALCULE” propose comme déplacement pour la prédiction à chaque pas de temps, le déplacement donné par une histoire mécanique précisée sous le mot clé EVOL_NOLI (par le mot-clef EVOL_NOLI). Le déplacement est projeté sur l’ensemble des champs cinématiquement admissible, comme pour la méthode EXTRAPOLE.

Remarques:

  • Les méthodes “EXTRAPOLE” et “DEPL_CALCULE” procèdent à une projection de la solution sur l’ensemble des champs cinématiquement admissibles. On se sert pour cela des conditions aux limites de Dirichlet donné dans le mot-clef EXCIT. Dans ce cas, il n’est pas possible d’utiliser des chargements de Dirichlet de type «cinématique» (opérande AFFE_CHAR_CINE) mais uniquement des chargements de Dirichlet par dualisation (opérande AFFE_CHAR_MECA). Le calcul fonctionne mais le risque dans ce cas étant que le champ de déplacement ne soit pas cinématiquement admissible, l’utilisateur doit y être attentif.

  • Les méthodes “EXTRAPOLE’et “DEPL_CALCULE’ne prennent pas en compte les termes de contact à la phase de prédiction. Par conséquent, il est conseillé de choisir les méthodes “ELASTIQUE” ou “TANGENTE” en prédiction pour des calculs de contact.

  • Il est nécessaire que les déplacements utilisés dans “EXTRAPOLE’et “DEPL_CALCULE” soient issus de calcul utilisant le même maillage car les conditions limites doivent être cohérentes.

  • Du fait de cette projection sur les conditions limites, ces deux options sont incompatibles avec les fonctionnalités du PILOTAGE.

Ces deux méthodes ont d’autres utilités. Par exemple, si on a obtenu une première solution sur le même maillage avec d’autres paramètres matériaux ou un autre comportement, les champs de déplacements peuvent être réutilisés dans le calcul.

De même, cela permet de réduire la place mémoire et de conserver ces résultats en vue d’une poursuite ultérieure. Pour un gros calcul, on peut stocker uniquement les déplacements à tous les instants dans IMPR_RESU. Si on veut recalculer les contraintes et variables internes, on fait un LIRE_RESU au format adéquat puis on utilise DEPL_CALCULE avec ITER_GLOB_MAXI=0 (on effectue une seule itération) et ARRET=”NON” (il n’y a pas convergence, on ne vérifie pas l’équilibre). Il est toutefois nécessaire pour des raisons de syntaxe de donner un chargement (éviter les chargements Dirichlet qui imposent une résolution linéaire) ainsi qu’un critère de convergence, même si ces informations ne sont pas prises en compte.

Avertissement

Du fait de l’impossibilité de «projeter» correctement les conditions limites d’un maillage à l’autre, il est désormais fortement déconseillé d’utiliser EXTRAPOLE et DEPL_CALCULE à partir d’un maillage différent du calcul courant.

5.2.2. Phase de correction#

La matrice utilisée pour les itérations globales de la méthode est la matrice tangente (consistante) [R5.03.01] si MATRICE=”TANGENTE” ou la matrice élastique si MATRICE=”ELASTIQUE”.

Le mot-clef REAC_INCR permet de contrôler la fréquence de réactualisation de la matrice par rapport aux pas de temps.

Remarques:

  • Si la matrice utilisée en prédiction correspond au calcul élastique, elle n’est évaluée qu’une fois à l’instant initial, en début d’algorithme, sauf si les caractéristiques élastiques dépendent du temps ou d’une variable de commande comme la température qui dépend du temps

  • Cette matrice « élastique » est calculée en utilisant le module d’Young donné sous le mot clé ELAS de l’opérateur [DEFI_MATERIAU], et non pas la pente à l’origine de la courbe de traction donnée sous le mot clé TRACTION (et qui sert, elle, dans l’expression des relations de comportement VMIS_ISOT_TRAC, VMIS_ECMI_TRAC et VISC_ISOT_TRAC)

Le mot-clef REAC_ITER permet de contrôler la fréquence de réactualisation de la matrice par rapport aux itérations de Newton. Donc à la première itération de Newton, on ne réassemble la matrice tangente que si it vaut 1: sinon on garde la matrice utilisée dans la phase de prédiction. Par convention si it vaut 0 la matrice n’est pas réévaluée durant tout le pas de temps.

5.2.3. Bascule de la matrice en correction#

Les options PAS_MINI_ELAS et REAC_ITER_ELAS permettent de passer de la matrice tangente à la matrice de décharge (i.e en considérant que les non-linéarités n’évoluent pas) lorsque le pas de temps est inférieur à PAS_MINI_ELAS. Cette matrice de décharge est la matrice élastique pour les modèles de comportement de type plastique; pour les modèles d’endommagement elle s’identifie à la matrice sécante.

Comme la convergence avec la matrice élastique est plus lente que celle avec la matrice tangente, le mot clé ITER_GLOB_ELAS sous le mot clé facteur CONVERGENCE permet de définir un nombre d’itérations maximal spécifique à l’utilisation de la matrice élastique et différent de celui associé à l’utilisation de la matrice tangente.

On peut définir une fréquence de réactualisation de la matrice de décharge avec le mot-clé REAC_ITER_ELAS (analogue de REAC_ITER). Si la matrice de décharge ne dépend pas de l’état de déformation (ce qui est le cas pour les matériaux plastiques mais pas pour les modèles d’endommagement), prendre REAC_ITER_ELAS = 0 (puisqu’elle sera la même au cours des itérations).

Ces options peut être utile lorsque le redécoupage automatique du pas de temps ne suffit pas à faire converger un calcul. Par exemple, dans le cas de lois adoucissantes, la matrice tangente peut devenir singulière et il vaut donc mieux utiliser la matrice élastique pour converger.

5.3. Méthode de Newton-Krylov#

Si on choisit METHODE=”NEWTON_KRYLOV”, on utilise une version inexacte de l’algorithme de Newton-Raphson; la précision des résolutions des systèmes linéaires par une méthode itérative est adaptée au cours de chaque pas de chargement (voir [R5.03.01]). Toutes les options décrites pour la méthode de Newton (gestion des matrices, de la réactualisation, etc.) sont utilisables par la méthode de Newton-Krylov.

5.4. Méthode IMPLEX#

Quand METHODE=”IMPLEX”, on utilise l’algorithme IMPLEX pour résoudre le problème (voir [ref:R5.03.81 <R5.03.81>]). Cette méthode n’est utilisable qu’avec un nombre limité de : comportements.

5.5. Méthode de réduction de modèle (ROM)#

◆ MODELE_REDUIT = _F(
      ◇ REAC_INCR = int (défaut: 1),
      ◇ PREDICTION = / "DEPL_CALCULE",
                    / "ELASTIQUE",
                    / "EXTRAPOLE",
                    / "TANGENTE" (par défaut),
      ◇ MATRICE = / "ELASTIQUE",
                  / "TANGENTE" (par défaut),
      ◇ REAC_ITER = int (défaut: 1),
      ◆ BASE_PRIMAL = mode_empi,
      ◇ DOMAINE_REDUIT = / "NON" (par défaut),
                        / "OUI",
      ◇ EVOL_NOLI = evol_noli,
      # Si: (equal_to("DOMAINE_REDUIT", 'OUI'))
          ◆ GROUP_NO_INTERF = grno,
          ◇ CORR_COMPLET = / "NON" (par défaut),
                          / "OUI",
          # Si: (equal_to("CORR_COMPLET", 'OUI'))
              ◇ COEF_PENA = float (défaut: 1000000.0),
              ◆ GROUP_NO_ENCASTRE = grno,
  ),

Quand METHODE=”MODELE_REDUIT”, on utilise une méthode de réduction de modèle pour faire le calcul non-linéaire (voir [R5.01.05]). Il est nécessaire d’avoir construit une base réduite préalablement (commande [DEFI_BASE_REDUITE]).

La réduction de modèle n’est pas compatible avec les fonctionnalités suivantes:

  • le contact et le frottement (toutes méthodes)

  • le pilotage,

  • la recherche linéaire,

  • les conditions limites dualisées (AFFE_CHAR_MECA): il faut utiliser AFFE_CHAR_CINE,

  • les modèles mixtes: la réduction de modèle est uniquement compatible avec des éléments finis en déplacement (pas de modèles de structure non plus)

On peut préciser les caractéristiques de la méthode de résolution du problème incrémental non-linéaire par une méthode de réduction de modèle. La stratégie de prédiction peut être réglée comme dans le cas de Newton avec le mot-clef PREDICTION. On peut régler la réactualisation des matrices avec les mots-clefs REAC_INCR et REAC_ITER.

La base réduite est donnée par le mot-clef BASE_PRIMAL. Cette base doit avoir été construite sur le même modèle et le même maillage que le calcul (mais les chargements, les matériaux et le comportement peuvent être différents).

Il est possible d’activer l’hyper-réduction qui utilise une méthode de type dEIM (discrete Empirical Interpolation Method) par DOMAINE_REDUIT=”OUI”. Dans ce cas, le calcul est réduit sur une zone du maillage restreinte (appelée RID) et construite à l’aide de l’opérateur [DEFI_DOMAINE_REDUIT].

Il faut donner le groupe de noeuds sur lequel est défini l’interface entre le RID et le reste du domaine à l’aide du mot-clef GROUP_NO_INTERF.

Pour améliorer la qualité des résultats en hyper-réduction, il est possible de faire une correction du calcul hyper-réduit par un calcul élément fini détaillé avec CORR_COMPLET=”OUI”. Pour cela, il est nécessaire de définir un groupe de noeuds servant à faire la liaison entre le RID et le reste du modèle avec le mot-clef GROUP_NO_ENCASTRE (défini par exemple dans la commande [DEFI_DOMAINE_REDUIT]). Ce groupe de noeuds permet d’imposer les conditions limites de Dirichlet nécessaire à ce que le problème corrigé soit bien défini. Ces conditions limites étant définies par pénalisation, le coefficient de pénalisation peut être changé par le mot-clef COEF_PENA.

Cette correction se faisant toujours sur le domaine réduit, elle ne sera pas aussi précise qu’un vrai calcul détaillé (domaine complet). Néanmoins, elle est très utile lorsque l’on change de manière significative les conditions limites entre la définition de la base réduite (phase off-line dans [DEFI_BASE_REDUITE]) et le calcul utilisant cette base réduite (phase on-line dans STAT_NON_LINE).

5.6. Recherche linéaire#

◇ RECH_LINEAIRE = _F(
      ◇ METHODE = / "CORDE" (par défaut),
                  / "MIXTE",
                  / "PILOTAGE",
      ◇ RESI_LINE_RELA = float (défaut: 0.1),
      ◇ ITER_LINE_MAXI = int (défaut: 3),
      ◇ RHO_MIN = float (défaut: 0.01),
      ◇ RHO_MAX = float (défaut: 10.0),
      ◇ RHO_EXCL = float (défaut: 0.009),
  ),

La recherche linéaire peut permettre d’améliorer la convergence de la méthode de Newton (Cf.:ref:R5.03.01<R5.03.01> pour plus de détails).

Remarque:

Il est déconseillé d’utiliser la recherche linéaire avec les déformations GROT_GDEP, pour les modélisations COQUE_3D et en présence de contact.

Ce mot-clef facteur RECH_LINEAIRE permet d’activer la recherche linéaire (voir [R5.03.01]).

La mot-clef METHODE permet de choisir la méthode de recherche linéaire, c’est-à-dire l’algorithme de recherche du zéro de la fonctionnelle. La méthode CORDE (par défaut) est la méthode la plus simple, c’est une méthode sécante unidimensionnelle.

La méthode MIXTE est plus élaborée et utilise une méthode sécante avec des bornes variables. Elle est plus efficace lorsque la fonctionnelle n’est pas strictement concave (problèmes avec endommagement ou THM par exemple).

La méthode PILOTAGE est réservée au pilotage de type DEFORMATION, PRED_ELAS et LONG_ARC (voir le paragraphe sur le pilotage ). C’est la seule méthode utilisable avec ce type de pilotage. Pour le pilotage de type DDL_IMPO, on peut utiliser CORDE ou MIXTE.

On donne le nombre d’itérations maximum à effectuer par ITER_LINE_MAXI et la précision à atteindre par RESI_LINE_RELA pour réaliser la convergence de la recherche linéaire.

Pour la méthode CORDE, Il n’est pas nécessaire de spécifier une précision ni un nombre d’itérations très élevés, la pratique montrant que deux ou trois itérations de recherche linéaire sont suffisantes. On peut donc se contenter de demander trois itérations avec la précision par défaut. L’utilisateur ne peut pas mettre plus de 999 itérations de recherche linéaire pour la méthode CORDE.

Par contre, pour la méthode MIXTE, sur des problèmes avec endommagement, plusieurs dizaines d’itérations sont souvent efficaces.

Les mots-clés RHO_MIN, RHO_MAX et RHO_EXCL fixent l’intervalle \(I\) dans lequel on calcule le coefficient \(\rho\): de la recherche linéaire, sous la forme

\[I=[\rho_{min},\rho_{max}]-[-\rho_{exc},\rho_{exc}]\]

5.7. Contrôle de la THM#

◇ SCHEMA_THM = _F(
        ◇ PARM_THETA = float (défaut: 1.0),
        ◇ PARM_ALPHA = float (défaut: 1.0),
    ),

Le mot-clef SCHEMA_THM permet de gérer les paramètres des schémas THM pour l’hydraulique (PARM_THETA, voir [R7.01.10]) et pour les volumes finis de type SUSHI (PARM_ALPHA, voir [R7.01.34]).

Ce paramètre PARM_THETA est bien distinct de celui défini dans le mot-clé COMPORTEMENT.