5. Nonlinear solver settings#

5.1. Choice of method#

◇ METHODE =/"IMPLEX ",
            /"MODELE_REDUIT ",
            /"NEWTON" (by default),
            /"NEWTON_KRYLOV ",

There are several categories of non-linear solvers available: the Newton method (and its variants), the Newton-Krylov method, the IMPLEX method, the access to the solver non-linear SNES from the PETSc library and the use of methods to reduce pattern (ROM). The solver is selected by the keyword METHODE.

5.2. Newton’s method#

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

When METHODE = » NEWTON « , we use the Newton-Raphson algorithm to solve the problem (see [R5.03.01]).

It is the keyword factor NEWTON (only one occurrence) that allows you to regulate the various settings. The MATR_RIGI_SYME keyword is used to force the symmetry of the stiffness matrix. It is useful when a feature using a non-symmetric matrix is used simultaneously (Simo-Miehe formalism, loadings of pressure type follower, etc.) and the formulation of contact type DISCRETE or LIAISON_UNIL (see [U4.44.11]). In fact, the latter can only function with a symmetric stiffness matrix. Making the stiffness matrix symmetric affects the speed of convergence, which can go so far as to fail but will never produce false results. If the fact of symmetrizing the stiffness matrix makes it impossible for the convergence of Newton’s algorithm, the user is invited to use a contact formulation of the type CONTINUE, which works with a non-symmetric stiffness matrix.

5.2.1. Prediction phase#

The prediction phase (Cf. [R5.03.01]) is intended to calculate an estimate of the displacement field in order to allow Newton’s method to converge more quickly. When the keyword is absent, it is the speed tangent matrix (option RIGI_MECA_TANG) that is used if one has chosen for Newton’s method a PREDICTION =” TANGENTE “, and it is the elastic matrix (option RIGI_MECA) which is used if PREDICTION =” ELASTIQUE “has been chosen.

It is also possible to make a prediction using preliminary results from of another calculation.

For PREDICTION =” EXTRAPOLE “, we calculate the estimate of the displacement increment from of the total increment obtained as a solution to the previous time step (weighted by the time step report). This estimate is projected onto all fields kinematically admissible (i.e. satisfying Dirichlet boundary conditions) according to the norm given by the elastic matrix, which must therefore be calculated. This feature is interesting in the case of using explicit local integration schemes of Runge-Kutta types that do not provide a tangent matrix: in this case the method of Newton uses an elastic matrix, but the number of iterations needed may be high. Using extrapolation can improve performance.

For PREDICTION =” DEPL_CALCULE “propose as a move for the prediction at each step of time, the displacement given by a mechanical story specified under the keyword EVOL_NOLI (with the keyword EVOL_NOLI). The displacement is projected on all fields that are kinematically admissible, as for the EXTRAPOLE method.

Notes:

  • The “EXTRAPOLE” and “DEPL_CALCULE” methods project the solution on all the kinematically admissible fields. For this purpose, conditions are used. at the Dirichlet limits given in the EXCIT keyword. In this case it is not possible to use « kinematic » Dirichlet loads (operand AFFE_CHAR_CINE) but only Dirichlet loads by dualization (operand AFFE_CHAR_MECA). The calculation works, but the risk in this case is that the field of movement is not not kinematically admissible, the user must be attentive to it.

  • The “EXTRAPOLE” and “DEPL_CALCULE” methods do not take into account contact terms at the prediction phase. Therefore, it is advisable to choose the methods “ELASTIQUE” or “TANGENTE” in prediction for contact calculations.

  • It is necessary that the moves used in “EXTRAPOLE” and “DEPL_CALCULE” be derived from calculations using the same mesh because the boundary conditions must be consistent.

  • Because of this projection on the boundary conditions, these two options are incompatible with the features of PILOTAGE.

These two methods have other uses. For example, if we obtained a first solution on the same mesh with other material parameters or other behavior, the fields of trips can be reused in the calculation.

Likewise, it makes it possible to reduce memory space and to keep these results in order to subsequent pursuit. For a big calculation, you can only store trips to all the time in IMPR_RESU. If we want to recalculate the constraints and variables internally, we do a LIRE_RESU in the appropriate format then we use DEPL_CALCULE with ITER_GLOB_MAXI =0 (we perform a single iteration) and ARRET =” NON “(there is no convergence, we do not check balance). However, it is necessary for syntax reasons to give a loading (avoid Dirichlet loads that require linear resolution) as well as a criterion of convergence, even if this information is not taken into account.

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. Correction phase#

The matrix used for the global iterations of the method is the tangent (consistent) matrix [R5.03.01] if MATRICE =” TANGENTE “or the elastic matrix if MATRICE =” ELASTIQUE”.

The REAC_INCR keyword allows you to control how often the matrix is updated. compared to time steps.

Notes:

  • If the matrix used in prediction corresponds to the elastic calculation, it is evaluated only once at the initial moment, at the start of the algorithm, unless the elastic characteristics depend on time or a control variable such as the temperature that depends on time

  • This « elastic » matrix is calculated using Young’s modulus given under the keyword ELAS of the [DEFI_MATERIAU] operator, not the original slope of the traction curve given under the keyword TRACTION (and who is used in the expression behavioral relationships (VMIS_ISOT_TRAC, VMIS_ECMI_TRAC, and VISC_ISOT_TRAC)

The REAC_ITER keyword allows you to control how often the matrix is updated. compared to Newton iterations. So at Newton’s first iteration, we only reassemble the tangent matrix if it is equal to 1: otherwise the matrix used in the prediction phase is retained. By convention, if it is equal to 0 the matrix is not reevaluated during the entire time period.

5.2.3. Matrix switch in correction#

The options PAS_MINI_ELAS and REAC_ITER_ELAS allow you to go from the tangent matrix to the discharge matrix (i.e. considering that the nonlinearities do not change) when the step The time is less than PAS_MINI_ELAS. This discharge matrix is the elastic matrix for plastic-type behavior models; for damage models it is identified to the secant matrix.

Since the convergence with the elastic matrix is slower than that with the tangent matrix, The keyword ITER_GLOB_ELAS under the keyword factor CONVERGENCE defines a number Maximum iterations specific to the use of the elastic matrix and different from that associated with the use of the tangent matrix.

We can define a frequency of updating the discharge matrix with the keyword REAC_ITER_ELAS (analogue of REAC_ITER). If the discharge matrix does not depend on the state of deformation (this is the case for plastic materials but not for models) of damage), take REAC_ITER_ELAS = 0 (since it will be the same during the iterations).

These options can be useful when automatically redividing the time step is not enough to converge a calculation. For example, in the case of softening laws, the tangent matrix can become singular and it is therefore better to use the elastic matrix to converge.

5.3. Newton-Krylov method#

If we choose METHODE =” NEWTON_KRYLOV “, we are using an inaccurate version of the algorithm of Newton-Raphson; the precision of the resolutions of linear systems by an iterative method is adapted during each loading step (see [R5.03.01]). All the options described for Newton’s method (management of matrices, updating, etc.) can be used by the Newton-Krylov method.

5.4. Method IMPLEX#

When METHODE =” IMPLEX “, we use the IMPLEX algorithm to solve the problem (see <R5.03.81>[ref: `R5.03.81 `]). This method can only be used with a limited number of: behaviors.

5.5. Model reduction method (ROM)#

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

When METHODE =” MODELE_REDUIT “, we use a model reduction method to do the calculation non-linear (see [R5.01.05]). It is necessary to have built a base previously reduced (command [DEFI_BASE_REDUITE]).

The model reduction is not compatible with the following features:

  • contact and friction (all methods)

  • piloting,

  • linear search,

  • dualized boundary conditions (AFFE_CHAR_MECA): AFFE_CHAR_CINE must be used,

  • mixed models: the model reduction is only compatible with finite elements on the move (no structure models either)

We can specify the characteristics of the method of solving the problem non-linear incremental using a model reduction method. The prediction strategy can be set as in Newton’s case with the keyword PREDICTION. You can set the updating of the matrices with the keywords REAC_INCR and REAC_ITER.

The reduced base is given by the keyword BASE_PRIMAL. This base must have been built on the same model and the same mesh as the calculation (but the loads, materials and the behavior may be different).

It is possible to activate hyper-reduction using a method such as d EIM (discrete Empirical Interpolation Method) by DOMAINE_REDUIT =” OUI “. In this case, the calculation is reduced to a restricted mesh area (called RID) and built using the [DEFI_DOMAINE_REDUIT] operator.

You must give the group of nodes on which the interface between RID and the rest of the is defined domain using the GROUP_NO_INTERF keyword.

To improve the quality of hyper-reduction results, it is possible to do a correction of the hyper-reduced calculation by a detailed finite element calculation with CORR_COMPLET =” OUI “. To do this, it is necessary to define a group of nodes used to make the connection between the RID and the rest of the template with the GROUP_NO_ENCASTRE keyword (defined for example in command [DEFI_DOMAINE_REDUIT]). This group of nodes makes it possible to impose Dirichlet boundary conditions necessary for the corrected problem to be well defined. These Since limit conditions are defined by penalty, the penalty coefficient can be changed with the keyword COEF_PENA.

Since this correction is always done on the reduced domain, it will not be as accurate. only a real detailed calculation (complete domain). However, it is very useful when we change significantly the boundary conditions between the definition of the reduced base (phase offline in [DEFI_BASE_REDUITE]) and the calculation using this reduced base (phase on-line in STAT_NON_LINE).

5.7. Control of the THM#

◇ SCHEMA_THM = _F (
        ◇ PARM_THETA = float (default: 1.0),
        ◇ PARM_ALPHA = float (default: 1.0),
    ),

The SCHEMA_THM keyword allows you to manage the parameters of THM diagrams for hydraulics (PARM_THETA, see [R7.01.10]) and for finite volumes of type SUSHI (PARM_ALPHA, see [R7.01.34]).

This PARM_THETA parameter is quite distinct from the one defined in the COMPORTEMENT keyword.