5. Construction and resolution of systems

An important part of the non-linear algorithm consists in solving linear systems, which are built in four steps:

• Calculation and assembly of loads;

• Calculation and assembly of the second members;

• Calculation and assembly of the matrix;

• Linear system resolution;

In this part, we will describe the different phases and the main routines to consider.

5.1. Systems to be solved

Currently, there are several different systems that are resolved in the op0070 operator. They will always have the following form (except for post-processing, which uses eigenvalue search operators):

(5.1)\[\begin{split} \ left [\ begin {array} {cc} K& {B} ^ {T}\\ B& 0\ end {array}\ right]\ mathrm {.} \ left\ {\ begin {array} {c}\ delta r\\\ delta s\ end {array}\ right\}\ mathrm {=}\ left\ {\ begin {array} {c} {c} {L} {L} _ {r} _ {r}\\\ {L} _ {s}\ end {array}\ right\}\end{split}\]

\(\left[\begin{array}{c}K\end{array}\right]\) is the stiffness matrix or a linear combination of matrices and \(\left[\begin{array}{c}B\end{array}\right]\) is the Dirichlet Lagrange matrix.

\(\left\{\begin{array}{c}\delta r\end{array}\right\}\) is the solution increment of the unknown nodal values of the system and \(\left\{\begin{array}{c}\delta s\end{array}\right\}\) is the increment of the Lagrange parameters associated with the dualization of the boundary conditions. Concretely, the « unknown nodal values » of the system can be displacements, speeds, accelerations, pressures, temperatures, etc. The second two members are also nodal values. Details of the systems are given in the reference documentation [R5.03.01] and [R5.05.05].

5.2. The Merimo routine

For the calculation of tangent matrices (options FULL_MECA, FULL_MECA_ELAS,, RIGI_MECA_TANG,,,, RIGI_MECA_ELAS,, RIGI_MECA, RIGI_MECA_IMPLEX) and internal forces (option RAPH_MECA), we systematically use the merimo routine. This routine prepares the input fields (numerous in this case), it is called both for the calculation of elementary vectors for internal forces (see § 5.5.2) but also for the calculation of tangent elementary stiffness matrices (see § 5.3).

5.3. Calculating matrices

Here we are considering elementary matrices (MEELEM) and assembled matrices (MEASSE). Just like for loading (§ 5.5.1), this calculation uses a two-stage system:

Creation of a**local list* of matrices to be calculated or assembled;

• Calculation or assembly of matrices;

For the phase of creating the list of matrices to be calculated or assembled, there are several routines (unlike loads that only use the nmchar routine). These routines are as follows:

Operations	Routine
Construction of the local list for matrices used in explicit dynamics	ndxprm
Construction of the local list for the matrices used in correction	nmcoma
Construction of the local list for the matrices used in prediction	nmprma
Construction of the local list for matrices used in post-processing (vibratory modes or buckling modes)	nmflma
Construction of the local list for constant matrices throughout the calculation	nminmc

For the phase of calculating or assembling the matrices, we find routines similar to the case of loading.

Operations	Routine
Addition of a matrix to be assembled and/or calculated in the local list Initializing the local list	nmcmat Nmcmat
Switching, calculation and/or assembly of matrices	nmxmat
Calculating matrices	nmcalm
Assembling the matrices	nmassm
Calculation and assembly of stiffness matrices	nmrigi

However, there is one notable exception: the calculation of non-linear tangent matrices uses a direct call by nmrigi because this requires additional parameters (call to merimo, see 5.2) that the other matrices do not need. So matrix type MERIGI uses nmrigi in place of nmcalm and nmassm.

If you want to modify a matrix, you must therefore:

• Add its calculation or assembly in the right place (ndxprm, nmcoma, nmprma, nmflma or nminmc) as appropriate, or even add a new routine;

• Add elementary calculation or assembly in nmcalm and nmassm;

The nmxmat and nmcmat routines for handling the local list should not be step changed.

5.4. Calculating the resulting matrix MATASS

Depending on the type of resulting matrix that one wishes to obtain (in prediction, in correction or for initial acceleration), three global routines are used that will take care of this construction. These are the nmprac, nmcoma, and nmprma routines. These three routines also manage the display (option for calculating the stiffness matrix) and the statistics (time spent in operations). They have the same basic diagram:

• Calculation or assembly of elementary matrices (MEELEM) into assembled matrices (MEASSE), (see § 5.3);

• Management of matrix update parameters (stiffness, mass, damping), nmchrm routine;

• Management of the type of the stiffness matrix (name of the option), nmchoi routine;

• Construction of the resulting matrix by linear combination of the other matrices. In dynamics, the various contributions in the resulting matrix (mass, damping and stiffness matrices) are combined using a multiplying coefficient depending on the time scheme used and the time step, these coefficients are recovered via the access routine NDYNRE (see § 2.4.2.17). It’s the nmmatr routine;

• Factorization (or not) of the resulting matrix by calling the preres routine;

Operations	Routine
Calculation of the resulting assembled matrix for calculating the initial acceleration	nmprac
Calculation of the resulting assembled matrix for the prediction phase	nmprma
Calculation of the resulting assembled matrix for the correction phase	nmcoma
Management of matrix update parameters (stiffness, mass, damping)	nmchrm
Stiffness matrix type management (option name)	nmchoi
Choice to re-create the numbering	nmrenu
Calculating the resulting matrix MATASS	nmmatr
Taking into account follower loadings with their multiplier function	ascoma
Taking into account the modification of the resulting matrix by touch/friction (method DISCRETE)	nmasfr

5.5. Calculation of the second limb

The calculation of the second member is the most different part from one linear system to another, it will contain:

• The given loads (see § 5.5.1)

• Internal forces arising from the integration of behavior (§ 5.5.2);

• Quantities related to Dirichlet limit conditions such as support reactions (§ 5.5.3);

• Quantities linked to order variables (§ 5.5.4);

• The different quantities related to touch/friction (no details in this document);

• Inertia and damping contributions for dynamics (§ 5.5.6);

In the case of dynamics, for multi-step schemes (Newmark complete with MODI_EQUI =” OUI “) we will also add the contributions of the previous time steps.

5.5.1. Calculation of loads

There are three ways to assess loads and two phases. The modes are as follows:

• <ACCI> load mode for evaluating the initial acceleration in dynamics;

• <FIXE> load mode for the evaluation of fixed loads that do not depend on the solution;

• <VARI> load mode for evaluating variable loads depending on the solution (follower loadings);

There are also two phases. Either we are in correction or we are in prediction. These phases are only useful for very specific cases: modal damping, method IMPLEX and modal impedance.

The main routine that calculates loads is nmchar. In this routine, depending on the functionalities activated (list_func_acti) and depending on the calculation mode/phase, we will cause the calculation of the elementary vectors (hat variable VEELEM) and their assembly into assembled vectors (hat variable VEASSE). To do this, a number of utility routines are used that manage lists local to NMCHAR.

Operations	Routine
Addition of a load to be assembled or calculated, with an optional option Initialization of the load list (initializations local lists of nmchar)	nmcvec Nmcvec
Switching, calculation or assembly of loads	nmxvec
Load calculation	nmcalv
Assembling the loads	nmassv

Some loads do not have an elementary calculation phase but only the creation of a vector assembled directly. If you want to add a loading, you must therefore:

• Define which phase and which mode will activate it;

• Add its calculation or its assembly in nmchar according to the activated functionality;

• Add elementary calculation and/or assembly in nmcalv and nmassv;

The nmxvec, nmcvec routines should not be modified.

5.5.2. Calculation of quantities related to internal forces

The calculation of the internal forces can be carried out in two ways:

• By integrating the law of behavior, this is calculation option RAPH_MECA;

• By re-using the constraints calculated elsewhere, this is calculation option FORC_NODA;

This distinction is important because it does not involve the same costs (the second case is much faster). Integrating behavior is a costly operation that involves modifying internal variables and constraints, and often solving a non-linear system locally (at each Gauss point). As the calculation of constraints is also necessary when evaluating tangent matrices, the behavior is also integrated when calculating tangent matrices (option FULL_MECA). In static terms, the prediction phase calculates the internal forces using option FORC_NODA. In dynamics, behavior is systematically integrated.

Operations	Routine
Calculation of elementary vectors for internal forces (integration of behavior)	nmfint
Assembly of elementary vectors for internal forces (integration of behavior)	nmaint

The calculation of option FORC_NODA is done by the load evaluation mechanism (§ 5.5.1).

5.5.3. Calculation of quantities linked to dualized limit conditions

The dualization of the Dirichlet boundary conditions (see [R3.01.01]) results in the calculation of two assembled vector-type quantities:

• The support reactions, by the product of the matrix of boundary conditions dualized with the Lagranges corresponding to degrees of freedom \(\left[B\right]\mathrm{.}\left\{\lambda \right\}\), this is the vector identified by CNDIRI;

• The value of the degrees of freedom imposed by dualization \(\left[B\right]\mathrm{.}\left\{u\right\}\). At convergence, this quantity will be equal to the given limit conditions \(\left\{{u}^{d}\right\}\), the vector is identified by CNBUDI;

Operations	Routine
Calculation of elementary vectors for support reactions	nmdiri
Assembly of elementary vectors for support reactions	nmadir
Calculation and assembly of quantities imposed by dualization	nmbudi

5.5.4. Calculation of quantities linked to order variables

The calculation of these quantities is done by the load evaluation mechanism (§ 5.5.1), type CNVCF0, CNVCF1 and CNVCPR.

5.5.5. Calculating constant quantities during calculation

A certain number of vectors are constant during the calculation, they then use a mechanism similar to loads (§ 5.5.1), using the nmxvec, nmcvec, nmcalv and nmassv routines, except that we use the nminvc routine instead of nmchar.

Operations	Routine
Calculation and assembly of constant vectors during the calculation	nminvc

5.5.6. Calculation of quantities related to dynamics — Inertia and damping forces

These quantities are calculated using the load evaluation mechanism (§ 5.5.1), type CNDYNA. These quantities do not exist in the elementary state since they are the result of a matrix/vector product.

Operations	Routine
Calculating dynamic quantities for the second member	ndfdyn
Calculation of inertial forces	Miner
Calculation of damping forces	nmhyst

According to the time diagrams, it is necessary to use specific multiplier coefficients in front of each of these quantities, these coefficients are recovered via the access routine NDYNRE (see § 2.4.2.17).

5.5.7. Resulting second members

All that remains is to build the various second members by linear combination of all the quantities detailed above. The main routines for building these second members are seven in number:

Operations	Routine
Calculating the second limb for correction	nmassc
Calculating the second member for prediction — Option DEPL_CALCULE or EXTRAPOLE	nmassd
Calculation of the second limb for initial acceleration	nmassi
Calculating the second limb for prediction	nmassp
Second limb calculation for prediction — Static case	nsassp
Second limb calculation for prediction — Implicit dynamic case	ndassp
Calculating the second limb for prediction — Explicit dynamic case	nmassx

All of these routines are based on the same principles:

• Retrieving the required quantities, already assembled or calculated locally in these routines, via the hat variable VEASSE (see § 5.5.1 to § 5.5.6);

• Retrieval of the various multiplier coefficients. In dynamics, according to the time patterns, it is necessary to use specific multiplier coefficients in front of each of these quantities, these coefficients are recovered via the access routine NDYNRE (see § 2.4.2.17);

• Construction of one or two second members by linear combination (use of standard vtaxpy routines). There are two second members in the case of piloting;

To improve readability, standardized utility routines are used to group together the most frequent cases.

Operations	Routine
Calculation of the components of the second member vector for Dirichlet loads, fixed over the time step, for the normal case and the pilot case	nmasdi
Calculation of the components of the second member vector for Neumann-type loads, fixed over the time step, for the normal case and the pilot case	nmasfi
Calculation of the components of the second member vector for Neumann-type loads, variables (followers) during the time step, for the normal case (no controlled case)	nmasva
Calculation of the components of the second member vector for the specific loads of the dynamics, variables (followers) during the time step, for the normal case (no controlled case)	ndasva
Calculation of the components of the second member vector for specific dynamic loads, calculation of the initial acceleration, variables (followers) during the time step, for the normal case (no controlled case)	nmacva

5.6. System resolution

The resolution is done via the nmresd routine, which takes as input the assembled resulting matrix (see § 5.4) and the appropriate second members (two second members if you have control).

Operations	Routine
System resolution	nmresd
System resolution on physical degrees of freedom (call for resolution)	nmreso
System resolution on generalized degrees of freedom	nmresg

The resolution routine will therefore return one or two solution fields DEPSO1 and DEPSO2, stored in the hat variable SOLALG.