Convergence criteria

If an incremental solution strategy based on iterative methods is to be effective, realistic criteria should be used for the termination of the iteration. At the end of each iteration, the solution obtained should be checked to see whether it has converged within preset tolerances or whether the iteration is diverging. The convergence tolerance determines the speed and accuracy of a calculation. If the criterion is too coarse, the solution may be quite inaccurate. On the other hand, a criterion which is too tight results in unnecessary computations. Similarly, an ineffective divergence check can terminate the iteration when the solution is not actually diverging or force the iteration to search for an unattainable solution.

In general, three types of convergence criteria are used to terminate the iterations:

• Criteria based on the residuum \(\boldsymbol{r}\)
• Criteria based on the solution \(\boldsymbol{d}\)
• Criteria based on the ''unbalanced energy''

In numgeo we offer the user two approaches to evaluate the convergence of the solution:

• Dof-based Controls

  Evaluate convergence for each physical field individually (Default)

  ---

  read more...

• Global Controls

  Evaluate convergence without distinguishing between different fields

  ---

  read more...

Recommended

Per default numgeo uses dof-based controls to evaluate convergence of the iterative solution process - which is also the recommended approach.

Norms and notation

In many cases, the Euclidean (\(l_2\)) norm is used to calculate the comparative quantities to assess convergence:

\[ \lVert \sqcup \rVert_2 = \sqrt{\sqcup_{i}^T\sqcup_{i}} = \left( \sum_{i_{dof}}^{n_{dof}} \sqcup_{i_{dof}}^2 \right)^{\frac{1}{2}}. \]

By using the Euclidean norm, the mean error over all degrees of freedom is controlled. In some cases, the infinite norm \(\lvert \sqcup \vert_\infty\) is used instead of the Euclidean norm, which is defined as follows

\[ \lvert \sqcup \vert_\infty = \text{max} \left(|\sqcup|\right) = \text{max} \left( \text{abs} \left(\sqcup,\alpha \right) \right) \]

General considerations

In general, convergence criteria take the following form:

\[ \sqcup^\alpha \leq \epsilon s^\alpha \]

Here, \(\sqcup^\alpha\) is a quantity related to field \(\alpha\) that determines whether the system is in equilibrium, i.e., whether convergence is achieved. Typically, \(\sqcup^\alpha\) is calculated as the Euclidean or infinite norm of the solution increment, the force (or flux), or the energy of field \(\alpha\). The parameter \(\epsilon\) represents the desired accuracy, and \(s^\alpha\) is a scaling factor that ensures the convergence criterion is independent of the problem's dimensions. Generally, \(s^\alpha\) is of the same type as \(\sqcup^\alpha\), but there is no universal consensus on its exact definition, and different programs use various versions.

Dof-based controls

In some cases, especially when dealing with multi-field problems, the flux component is composed of different components \(\boldsymbol{f} = [\boldsymbol{f}^\alpha]^T\) which are conjugated to the \(\alpha\)th variable in the problem. Analogous, the residuum reads \(\boldsymbol{r} = [\boldsymbol{r}^\alpha]^T\). For the example of a three-phase porous media consisting of the solid matrix (\(s\)), the pore water (\(w\)) and the pore air (\(a\)), the variables to be sought of are the solid displacement \(\boldsymbol{u}\), the pore water pressure (\(p^w\)) and the pore air pressure (\(p^a\)), i.e. \(\alpha = \{s,w,a\}\). The vector of solution thus reads \(\boldsymbol{d} = [\boldsymbol{d}^s,\boldsymbol{d}^w,\boldsymbol{d}^a]^T = [\boldsymbol{u}, \boldsymbol{p}^w, \boldsymbol{p}^a]^T\). The flux vector and the residual now read \(\boldsymbol{f}=[\boldsymbol{f}^s,\boldsymbol{f}^w,\boldsymbol{f}^a]^T\) and \(\boldsymbol{r}=[\boldsymbol{r}^s,\boldsymbol{r}^w,\boldsymbol{r}^a]^T\), respectively.

Using the global convergence controls might not be a suitable choice due to the (very) different magnitudes of the different degrees of freedom (e.g. displacements in the order of cm, and pore water pressure in the order of several hundred kPa). In such cases controlling the accuracy of the solution based on each individual type of degree of freedom (e.g. solid displacements or pore water pressure) might be more appropriate. This is especially useful in multi-field solutions.

The dof-based convergence evaluation in numgeo uses all three types of convergence criteria based on

• the residuum \(\boldsymbol{r}^\alpha\)
• the solution \(\boldsymbol{d}^\alpha\)
• the ''internal energy'' of field \(\alpha\)

However, only three of the above criteria must be met (simultaneously) for a solution to be considered converged. In some cases quadratic convergence of the iterations is not possible. In these cases, we use a temporary relaxation of the convergence criteria as described in subsection Nonquadratic convergence. Notice that if the response for field \(\alpha\) is very small, the contribution of this field to the convergence of the overall solution is not considered, i.e. neither the control of the residuum nor of the unbalanced energy will be considered in this case. Therefore, great caution is required in the choice of the criterion for negligible response \(\epsilon^{\Delta d,\alpha}\). A graphical representation of some of the quantities used to evaluate convergence are given in Figure 1

Figure 1: A simplified graphical illustration depicting key quantities in the solution process for the response of a non-linear spring with a single degree of freedom to applied loads.

Solution increment

To limit the maximum error for every type of degrees-of-freedom (e.g. displacements, pore pressures,...), rather than the mean error over all degrees of freedom, the following expression is used:

\[ \label{eq:convergence-solution} \frac{c^{(i),\alpha,\text{max}}}{^{t}\Delta d^{(i),\alpha,\text{max}}} \leq \epsilon^{d,\alpha}, \]

where

• \(c^{(i),\alpha,\text{max}}\) is the largest change in solution (in absolute value) in iteration \(i\):

  \[ c^{(i),\alpha,\text{max}} = \Big\lvert \boldsymbol{c}^{(i),\alpha} \Big\rvert_\infty \]

• \(^{t}\Delta d^{(i),\text{max}}\) is the largest change in solution (in absolute value) in the present increment, i.e. the solution at the end of iteration \(i\) of the present increment minus the solution at the start of the increment:

  \[ ^{t}\Delta d^{(i),\alpha,\text{max}} = \Big\lvert ^{t}\Delta d^{(i),\alpha} \Big\rvert_\infty \]

  Note that \(^{t}\Delta d^{(i),\alpha}\) accounts for changes in the solution resulting from the prescribed Dirichlet boundary conditions.

• \(\epsilon^{d,\alpha}\) is the tolerance criteria used to check for convergence.

Evaluating the accuracy of the solution based on Eq. (\(\ref{eq:convergence-solution}\)) bears the problem, that if the response for field \(\alpha\) is very small, Eq. (\(\ref{eq:convergence-solution}\)) will not be satisfied although the contribution of field \(\alpha\) to the overall solution is negligible. Therefore, the contribution of this field to the convergence of the overall solution is only considered if the following condition holds:

\[ \text{max}\left(c^{(i),\alpha,\text{min}}, ~ ^{t}\Delta d^{(i),\alpha,\text{max}} \right) \geq \epsilon^{\Delta d,\alpha}. \]

Therein, \(\epsilon^{\Delta d,\alpha}\) is the criterion for negligible response for field \(\alpha\).

Unbalanced force

For the dof-based control by means of ''unbalanced force'', the Euclidean norm of the residuum \(\lVert \boldsymbol{r}^{(i),\alpha} \rVert_2\) at iteration \((i)\) is compared to the maximum Euclidean norm of the external forces \(^{t+\Delta t}\boldsymbol{f}^\alpha_{ext}\), the internal forces \(^{t+\Delta t}\boldsymbol{\tilde{f}}^\alpha_{int}\) and the inertia forces \(^{t+\Delta t}\boldsymbol{f}^\alpha_{dyn}\) of field \(\alpha\) at the end of the current increment:

\[ \Big\lVert ^{t}\boldsymbol{r}^{(i),\alpha} \Big\rVert_2 \leq \epsilon^{r,\alpha} \max \left( \Big\lVert~ ^{t+\Delta t} \boldsymbol{f}^\alpha_{ext} ~\Big\rVert_2, \Big\lVert~ ^{t+\Delta t} \boldsymbol{\tilde{f}}^\alpha_{int} ~\Big\rVert_2, \Big\lVert ~ ^{t+\Delta t} \boldsymbol{f}^\alpha_{dyn} ~\Big\rVert_2 \right). \]

Therein, \(\epsilon^{r,\alpha}\) is the criterion for negligible unbalance force. \(^{t+\Delta t} \boldsymbol{\tilde{f}}^\alpha_{int}\) are the internal forces of all degrees of freedom of field \(\alpha\) which are not constrained by Dirichlet boundary conditions. \(^{t+\Delta t} \boldsymbol{f}^\alpha_{ext}\) includes the reaction forces required to enforce Dirichlet boundary conditions.

Unbalanced energy

The energy convergence criterion measures the energy flow to the system resulting from the residual, which is like an error in energy. As a reference energy, the energy resulting from the internal forces is used.

\[ \frac{\Big\lvert \left({^t}\boldsymbol{r}^{(i),\alpha} \right)^T \boldsymbol{c}^{(i),\alpha} \Big\rvert}{ \Big\lvert \Delta {^t}W^{i=0,\alpha}_{int} \Big\rvert } \leq \epsilon^{e,\alpha}, \]

The internal energy is integrated using the trapezoidal rule:

\[ \Delta {^t}W^{(i=0),\alpha}_{int} = \dfrac{1}{2}\left( ^{t_0}\boldsymbol{f}^\alpha_{int} + ^{t}\boldsymbol{f}^{(i=0),\alpha}_{int} \right)^T \Delta \boldsymbol{d}^{(i=0),\alpha} \]

Similar to the control of the solution increment, the unbalance energy of field \(\alpha\) is only considered if the following condition holds:

\[ \Big\lvert \left({^t}\boldsymbol{r}^{(i),\alpha} \right)^T \boldsymbol{c}^{(i),\alpha} \Big\rvert \geq \epsilon^{\Delta e,\alpha}. \]

Therein, \(\epsilon^{e,\alpha}\) is the criterion for negligible unbalance energy for field \(\alpha\).

Constrained dofs

When calculating energy measures, we exclude degrees of freedom constrained by Dirichlet boundary conditions, as their values are given and remain constant throughout the iteration process, resulting in no corrections. Including these degrees of freedom would increase the initial energy used as a reference for the convergence check, making the criterion less strict.

Nonquadratic convergence

In some cases quadratic convergence of the iterations is not possible. To avoid unnecessary iteration, after a certain number \(n^{iter}_{nc}\) of iterations without convergence, numgeo uses looser tolerance criteria \(^{*}\epsilon^{e,\alpha}\) and \(^{*}\epsilon^{\Delta d,\alpha}\) for the dof-based solution controls:

\[ ^{*}\epsilon^{d,\alpha} = r^\epsilon \epsilon^{d,\alpha} ~~~~\text{and}~~~~ ^{*}\epsilon^{e,\alpha} = r^\epsilon \epsilon^{e,\alpha} \]

therein, \(r^\epsilon > 1\) is a relaxation factor, which can be prescribed by the user (default \(r^\epsilon = 2\)). The number of iterations without convergence \(n^{iter}_{nc}\) needed for this relaxation can be modified by the user, the default value is \(n^{iter}_{nc}=9\).

Default values

The default values for the dof based convergence criteria are:

dof	\(\epsilon^{d}\)	\(\epsilon^e\)	\(\epsilon^{r}\)	\(\epsilon^{\Delta d}\)	\(\epsilon^{\Delta e}\)	\(r^\epsilon\)
u	\(1 \cdot 10^{-2}\)	\(1 \cdot 10^{-2}\)	\(1 \cdot 10^{-2}\)	\(1 \cdot 10^{-5}\)	\(1 \cdot 10^{-6}\)	2
w	\(1 \cdot 10^{-2}\)	\(1 \cdot 10^{-2}\)	\(1 \cdot 10^{-2}\)	\(1 \cdot 10^{-5}\)	\(1 \cdot 10^{-5}\)	2
pw	\(5 \cdot 10^{-2}\)	\(5 \cdot 10^{-2}\)	\(5 \cdot 10^{-2}\)	\(1 \cdot 10^{-3}\)	\(1 \cdot 10^{-3}\)	2
pa	\(5 \cdot 10^{-2}\)	\(5 \cdot 10^{-2}\)	\(5 \cdot 10^{-2}\)	\(1 \cdot 10^{-3}\)	\(1 \cdot 10^{-3}\)	2

Global convergence controls

We evaluate the global (overall) convergence without distinguishing between different fields \(\alpha\) (e.g. solid displacement and pore water pressure). This measure gives an idea of the overall (mean) quality of the solution, however, especially in multi-field solutions, this measure could shadow non-converged fields or lead to non-convergence due to the different fields with significantly different magnitudes.

Unbalanced energy

\[ \frac{\left| \left({^t}\boldsymbol{r}^{(i)} \right)^T \boldsymbol{c}^{(i)} \right| }{ \left| \left({^t}\boldsymbol{r}^{(i=0)} \right)^T \Delta \boldsymbol{d}^{(i)} \right| } \leq \epsilon^{e}, \]

Unbalanced forces

For the global convergence control, the Euclidean norm of the residuum \(\boldsymbol{r}\) at iteration \((i)\) is compared to a load norm, which is defined as the maximum of the external \(\boldsymbol{f}_{ext}\), internal \(\boldsymbol{f}^{(i)}_{int}\) and inertia \(\boldsymbol{M^{(i)}\ddot{d}^{(i)}}\) contributions to the residuum:

\[ \lVert ^{t}\boldsymbol{r}^{(i)} \rVert_2 \leq \epsilon^r \text{max} \big( \lVert ^{t}\boldsymbol{f}_{ext} \rVert_2~,~ \lVert ^{t}\boldsymbol{f}_{int} \rVert_2~,~ \lVert ^{t}\boldsymbol{M\ddot{d}} \rVert_2 \big). \]

Therein, \(\epsilon^{r}\) is the criterion for negligible unbalance force. Note that in most cases, the external load vector is only calculated at the beginning of the current load step and does not change during the iterative procedure.

Solution increment

The third global solution control checks if the \(L_2\) norm of the correction \(\boldsymbol{c}^{(i)}\) is smaller than a predefined fraction \(\epsilon^d\) of the L\(_2\) norm of the solution increment \(\Delta \boldsymbol{d}^{(i)}\):

\[ \boldsymbol{c}^{(i)} \leq \epsilon^d \Delta \boldsymbol{d}^{(i)} \]

Default values

The default values for the dof based convergence criteria are:

\(\epsilon^{d}\)	\(\epsilon^e\)	\(\epsilon^{r}\)
\(1 \cdot 10^{-2}\)	\(5 \cdot 10^{-3}\)	\(5 \cdot 10^{-3}\)