This overview of the TFEL/Math library mostly focuses on two important features:

• The TFEL/Math library provides a linear algebra engine with mathematical objects and operations on those objets required to express the constitutive equations in an efficient and natural manner, i.e. as close as possible to the mathematical expressions. Section 1 is devoted to describe this linear algebra engine.
• The TFEL/Math library provides efficient and robust non linear solvers for small sized problems described in Section 2.

# 1 Mathematical objects

The mathematical objects provided by the TFEL/Math library are based on scalar values which can be standard C++ numeric types (float, double or long double), or so-called quantities, i.e. an object describing a value with an unit. The later allows to distinguish physical quantities such as a strain from a stress and prevent illegal operations between them. Quantities are described in depth in Section 1.1.

The TFEL/Math library then introduces some mathematical objects of rank 1, i.e. objects which can be stored in a contiguous memory location and access by a single index, such as:

• fixed size vectors (implemented by the tvector class)
• symmetric tensors (implemented by the stensor class)
• unsymmetric tensors (implemented by the tensor class)

Higher order objects are then defined as the derivatives of objects of lower ranks. For instance, the derivative of a symmetric tensor with respect to an unsymmetric tensor is implemented by the t2tost2 class. The latter can also be viewed a linear application between unsymmetric tensor and symmetric tensor or as a fourth order tensor. In practice, there are four kind of fourth order tensor defined in the library respectively implemented by the st2tost2, st2tot2, t2tost2, and t2tot2 classes.

Currently only tensors up to the fourth order are specialized which allows to provide some specific features.

Higher order tensors may be represented by a generic class which represent derivatives of two mathematical objects. Currently, the number of operations allowed by this generic class is rather limited and will be extended as needed in future versions.

The main concepts behing the implementation of the mathematical objects are described in Section 1.2.

## 1.1 Quantities

Quantities were introduced in TFEL/Math following ideas of the book D. Abrahams and A. Gurstovoy  to allow dimensional analysis at compile time. Basically, a quantity is the combinaison of a value and an unit. The type of the value is called the base type of the quantity.

The unit is encrypted in a type called Unit so that the C++ type-system can be used to detect at compile-time (without runtime-checks) if operations are legal. According to the International System of Units (SI), an unit is decomposed along $$7$$ basic units which are the mass (kg), the length (l), the time (s), the electric current (A), the temperature (K), the luminous intensity (cd) and the amount of substance (mole). To be able to describe fractional unit, such as the fracture toughness which is has the unit of a stress times the square root of a length, the Unit class is parametrized by $$14$$ integers. The TFEL/Math library provides some convenient aliases for the most common units.

The NoUnit unit is a special case where all those integers are set to 0 for the numerator parts and 1 (by convention) for the denomiator parts.

### 1.1.1 The qt class

A quantity is presented by the qt class which is parametrized by a class describing the unit and a numeric base type (generally float, double or long double) which stores the value of the quantity. By default, the double type is used. This class can be used as follows:

constexpr qt<Mass> m1(100.);
constexpr qt<Mass> m2(100.);
constexpr qt<Mass> m3 = m1 + 0.5 * m2;
constexpr qt<Acceleration> a(2);
constexpr qt<Force> f = m1 * a;

The previous code sample illustrates how to declare a new quantity and how to perform operations on quantities. It also highlights that those operations can be used in a constexpr context, i.e. that at compile-time.

### 1.1.2 The qt_ref and const_qt_ref classes

The library also provides two wrapper classes called qt_ref and const_qt_ref which allows to wrap a value into an object which acts respectively as an mutable an immutable quantity, as illustrated in the following code sample:

constexpr auto eps = 1e-14;
constexpr auto value = [] ()constexpr {
constexpr qt<Mass> m(100.);
constexpr qt<Acceleration> a(2);
auto vf = double{};
auto f = qt_ref<Force>(vf);
f = m * a;
return vf;
}();
TFEL_TESTS_STATIC_ASSERT(std::abs(value - 200.) < eps);

While a bit contrieved, this example shows that the qt_ref class can also be used in a constexpr context.

The qt_ref and const_qt_ref classes are parametrized by two template parameters which describe respectively the unit and a numeric base type.

The qt_ref and const_qt_ref classes are essentially used to build views of mathematical objects from a raw memory area, as described in Section 1.7.

### 1.1.3 The Quantity class

For the sake of simplicity, the qt, qt_ref and const_qt_ref were described as classes. This is not actually the case in the current implementation which defines them through type aliases to a more general Quantity class which is parametrized by three template arguments: the unit, the basic numeric type and a third argument describing how the value associated with the quantity is handled.

### 1.1.4 Operations on quantities

Common operations on quantities, such as additions, substraction, multiplications, division and negations are supported.

Scaling and multiplication by a standard numeric value works as expected, as demonstrated by the following code:

constexpr qt<Mass> m(100.);
// multiplication by a raw numeric type
auto m2 = 2. * m;
// scaling
m2 *= 4;

A quantity with unit NoUnit is implicitly convertible to its base type:

constexpr qt<NoUnit> q(1.2);
constexpr double q_value = q;
const auto cos_q = std::cos(q);

The latter line shows that standard mathematical functions can be called directly.

Fractional power of a quantity are supported through the power function as follows:

// declaration of a stress symmetric tensor
stensor<3u,qt<Stress>> s;
// contracted product of s which has the unit
// of the square of a stress
const auto s2 = s|s;
// norm of s
const auto norm = power<1,2>(s2);

About the power function

The power function is a very convenient function which takes two template parameters N and D defining the exponent as the fraction $$N/D$$. The default value of D is $$1$$, so the following syntax is valid:

auto v = power<3>(12);

which computes $$12^{3}$$.

The power function is optimised for the following values of the denominator:

• $$1$$: In this case, the Nth power of the argument is computed using by multiplications if $$N$$ is positive. If $$N$$ is negative, the Nth power of the inverse of the argument is computed.
• $$2$$: In this case, the square root of the Nth power of the argument is computed using the std::sqrt function.
• $$3$$: In this case, the square root of the Nth power of the argument is computed using the std::cbrt function.

## 1.2 A generic framework to build tensorial-like objects

This section first introduces the notion of array policy and indexing policy which describes respectively how data are handled (stored) and how those data are accesssed.

Common objects of rank $$1$$ are then defined, including tiny vectors, symmetric tensors and unsymmetric tensors. Higher order objects are defined recursively as derivative of objects of lower ranks.

### 1.2.1 Description of the ArrayPolicy concept

Array policies describes how data are handled and accessed. For the latter, the description of the access partern is delegated to an indexing policy which will be described later in Section 1.3. Thus, this section mostly focuses on how array policies describe how the data are handled.

A class matching the ArrayPolicy concept must provide:

• An alias to a type called value_type which corresponds to the logical type hold by the array (see below).
• An alias to a type called const_value_type which corresponds to a constant value to a an object of the logical type hold by the array.
• An alias to a type called storage_type which is the type actually stored in memory.
• An alias to a type called reference which is the type returned by the non-constant access operators of the array. A object of the reference type must be convertible and assignable to a value_type object and provide the following operators:
• operator=(const other_type&)
• operator+=(const other_type&)
• operator-=(const other_type&)
• operator*=(const other_type&)
• operator/=(const other_type&)
where other_type denotes for each operator a type for which the considered operator in the value_type class is valid.
• An alias to a type calle const_reference which is the type returned by the constant access operators of the array. A object of the const_reference type must be convertible and assignable to a value_type object.
• A static function called make_const_reference which turns a a constant reference to a storage_type object into a const_reference object.
• A static function called make_reference which turns a reference to a storage_type object into a reference object.
• An alias type a type called IndexingPolicy which refers an indexing policy. The IndexingPolicy concept is described in depth in Section 1.3.

It is worth illustrating those concepts in two commonly used cases:

• Standard array policies.
• Array policies used by views.

### 1.2.2 Standard array policies

The StandardArrayPolicy class is based in the following statements:

• storage_type is the same as value_type.
• reference is the same as as value_type&.
• const_reference is the same as as const value_type&.
• make_reference is trivial an only returns a const reference to the stored value.
• make_const_reference is trivial an only returns an immutable reference to the stored value.

Standard array policies are which are used to implement concrete mathematical objects such as tvector, tmatrix, stensor, tensor, etc.

Note

Contrary to view array policies described in the next paragraph, the storage_type is not restricted to be a basic numeric type (float, double or long double).

### 1.2.3 View array policies

Views generally maps a memory location made of basic numeric types (float, double or long double) to a mathematical objects.

To support quantities, the type stored (storage_type) must be distinguished from the logical type of the mapped object (value_type). The reference type is then a simple wrapper around a numeric value which acts as the desired quantity.

For example, let us consider the case of a view mapping a stress symmetric tensor from an array of double precision numbers. In this case, storage_type is the base type of value_type, so:

• storage_type is an alias to double.
• value_type is an alias to qt<Stress, double>.
• reference is an alias to qt_ref<stress, double>.
• const_reference is an alias to const_qt_ref<stress, double>.
• make_reference simply builds a reference for the address of the stored value. So does make_const_reference.

This mostly describes the implementation of the ViewArrayPolicy class.

However, those rules only applies when storage_type is a quantity. Otherwise, i.e. when value_type is a basic numeric type, the implementation of the ViewArrayPolicy class is equivalent to the StandardArrayPolicy class.

## 1.3 Description of the IndexingPolicy concept

A class matching the IndexingPolicy concept must provide:

• An alias to a type called size_type to the type used to index data.
• An alias to a type called RowMajorIndexingPolicy which can be used to import data from raw C pointers.
• A static member data called arity of type size_type which describes the number of indices required to describe the object.
• A static boolean data member called hasFixedSizes. If this data member is true, the class must be stateless, i.e. empty.
• A static boolean data member called areDataContiguous which states if the data are continous
• A member function called size which takes no argument returning the number of data values accessible by the indexing policy. This member function must be constexpr if hasFixedSizes is true.
• A member function called size which takes a size_type argument returning range in the given dimension. This member function must be constexpr if hasFixedSizes is true.
• A member function called getUnderlyingArrayMinimalSize which returns the minimal size of an array required by the indexing policy. This member function must be constexpr if hasFixedSizes is true.
• A member function called getIndex which returns the memory offset of an element of the object with respect to the location of the first element of the object. The number of argument of this method must be equal to the arity member. This member function must be constexpr if hasFixedSizes is true.

### 1.3.1 An example of indexing policy Figure 1: Storage of the elements of a matrix using the row-major format

Let us consider an $$N\,\times\,M$$ matrix stored in a contiguous memory location using the row major format, i.e. all elements of the same row are stored continuously as depicted on Figure 1. The position of $$i,j$$ element from the start of the memory area is given by $$i \, M + j$$.

Let us now consider a view on a $$L\,\times\,K$$ submatrix in a which starts at row $$i0$$ and column $$j0$$. The position of the $$i,j$$ element of the submatrix in the matrix is also given by $$i \, M + j$$ if we start from the first element of the submatrix.

This formula is implemented by the getIndex method of the FixedSizeRowMajorMatrixIndexingPolicy which is used by many classes provide by the library, such as tmatrix (matrices whose size is known at compile time), st2tost2, t2tot2, etc.. and views to those objects in a matrix.

### 1.3.2 Higher order indexing policies for fixed-size mathematical objects

Let us consider an mathematical object $$o_{3}$$ resulting from the derivation of an mathematical object $$o_{1}$$ of arity $$a_{1}$$ with respect to an mathematical object $$o_{2}$$ of arity $$a_{2}$$. $$o_{3}$$ has an arity of $$a_{1}+a_{2}$$.

The storage of those objects are described respectively by the indexing policies $$p_{1}$$ and $$p_{2}$$. We assume that the object $$o_{1}$$ can be stored in a memory area of size $$s_{1}$$ and that the object $$o_{2}$$ can be stored in a memory area of size $$s_{2}$$. Then the object $$o_{3}$$ can be stored in a memory location $$s_{1}\,s_{2}$$.

Then, an indexing policy $$p_{3}$$ suitable to describe $$o_{3}$$ may compute the position of the derivative of the component $$o_{1}{\left(i_{0},\ldots,i_{a_{1}-1}\right)}$$ with respect to the component $$o_{2}{\left(j_{0},\ldots,j_{a_{2}-1}\right)}$$ is given by:

$p_{3}{\left(i_{0},\ldots,i_{a_{1}-1},j_{0},\ldots,j_{a_{2}-1}\right)}= p_{1}{\left(i_{0},\ldots,i_{a_{1}-1}\right)}\,s_{2}+p_{2}{\left(j_{0},\ldots,j_{a_{2}-1}\right)}$

This choice is implemented in the FixedSizeIndexingPoliciesCartesianProduct class which is used to by the FixedSizeArrayDerivative class to describe derivatives of two arbitrary mathematical objects.

## 1.4 Low rank mathematical objects

The library is based on a few low rank mathematical objects:

• tvector: a tiny vector of fixed size. This template class is parametrized by the size of the vector and the value hold.
• stensor: a symmetric tensor. This template class is parametrized by the space dimension and the scalar type.
• tensor: a unsymmetric tensor. This template class is parametrized by the space dimension and the scalar type.

All those objects are represented by objects of rank one using a vector-like notations. For example, a $$3D$$ symmetric tensor is represented as follows:

$\underline{s}= \begin{pmatrix} s_{\,11}\quad s_{\,22}\quad s_{\,33}\quad \sqrt{2}\,s_{\,12}\quad \sqrt{2}\,s_{\,13}\quad \sqrt{2}\,s_{\,23} \end{pmatrix}^{T}$

This notations has the property that the contracted product of two symmetric tensors is the scalar product of this vector form (hence the $$\sqrt{2}$$).

In a similar manner, an unsymmetric tensor is represented as follows: $\underline{s}= \begin{pmatrix} s_{\,11}\quad s_{\,22}\quad s_{\,33}\quad s_{\,12}\quad s_{\,21}\quad s_{\,13}\quad s_{\,31}\quad s_{\,23}\quad s_{\,32} \end{pmatrix}^{T}$

## 1.5 Higher order objects defined as derivatives

The library provides a generic class called FixedSizeArrayDerivative which allows to create higher order objects as being the derivative of two objects of lowest ranks.

This class is currently very limited but will be extended in future versions of the library.

### 1.5.1 The case of fourth order tensors

Fourth order tensors can be defined as derivatives of two tensors or as linear mappings from the second order tensors to second order tensors.

As there are two kinds of second order tensors (i.e. symmetric and non symmetric tensors), there are four kinds of fourth order tensors defined in the TFEL/Math library, which satisfy the following concepts:

• ST2toST2Concept: linear mapping from symmetric tensors to symmetric tensors.
• ST2toT2Concept: linear mapping from symmetric tensors to non symmetric tensors.
• T2toST2Concept: linear mapping from non symmetric tensors to symmetric tensors.
• T2toT2Concept: linear mapping from non symmetric tensors to non symmetric tensors.

An end user will mostly use the following implementations of those concepts: st2tost2, st2tot2, t2tost2 and t2tot2 respectively. Those classes have the following template arguments:

• The space dimension (1, 2 or 3).
• The type used to perform the computation.

### 1.5.2 The tmatrix case

A fixed size matrix, implemented by the tmatrix class can be seen as the derivative of a tiny tensor with respect to a tiny vecor.

### 1.5.3 The derivative_type metafunction

The TFEL/Math library provides a very convenient type alias called derivative_type which automatically selects the correct type as the derivative of two objects of fixed sizes. This alias also works with scalars. This type alias takes quantities into account if required.

Here are a few examples:

• derivative_type<double, double> is an alias to double.
• derivative_type<qt<Stress>, double> is an alias to qt<Stress>.
• derivative_type<qt<Stress>, qt<Time>> is an alias to qt<StrainRate>.
• derivative_type<tvector<2u, double>, tvector<3u, double>> is an alias to tmatrix<2u, 3u, double>.
• derivative_type<stensor<2u, double>, stensor<2u, double>> is an alias to st2tost2<2u, double>>.
• derivative_type<stensor<2u, double>, tvector<3u, double>> is an alias to FixedSizeArray Derivative<stensor<2u, double>, tvector<3u, double>>.

## 1.6 Expressions templates

One may expect that the addition of two vector results in a new vector. This naive approach may lead to poor performances, due to temporaries objects and data copying .

Expression templates is a C++ template metaprogramming technique which introduces additional classes which represent the actions to be performed on some objects and lazily delay the execution of those actions until the result is explicitly requested.

The objects of those classes are placeholders, also called handlers within the library, that are meant to be assigned to an object whose type is the expected result of the operation treated.

To illustrate this technique, let us consider the addition of three vectors $$\vec{a}$$, $$\vec{b}$$ and $$\vec{c}$$ and its assignment to a vector $$\vec{d}$$:

d = a + b + c;

The addition of the vectors $$\vec{a}$$ and $$\vec{b}$$ produces an intermediate object $$\vec{e1}$$ of type Expr1 which keeps a reference to those two vectors. Similarly, the addition of those three vectors defines another object $$\vec{e2}$$ of type Expr2 which stands for the addition of $$\vec{e1}$$ and the vector $$\vec{c}$$.

Figure 2 shows how the access operator of e2 is implemented. The assignment of an object of type Expr2 to a vector d is implemented as a standard for loop:

for(size_type i=0;i!=d.size();++i){
d(i) = e2(i);
}

The temporary objects e1 and e2 are meant to be eliminated by the compiler optimisation process. Thanks to function inlining, the compiler is able to produce a code that is equivalent to what would have been obtained with the following instructions:

for(size_type i=0;i!=d.size();++i){
d[i] = a[i] + b[i] + c[i];
}

When dealing with objects whose size is known at compile-time, the TFEL/Math library also performs an additional optimisation technique known as loop unrolling. For example, if $$\underline{a}$$, $$\underline{b}$$ and $$\underline{c}$$ are three $$1D$$ symmetric tensors, the code:

c = a + 2 * b;

is equivalent to:

c = a + 2 * b;
c = a + 2 * b;
c = a + 2 * b;

### 1.6.1 The eval function

One possible caveat of the expression template technique can be illustrated by the following example. Let us consider two vectors a and b:

a = 1;
b = 2;
auto c = a + b;
// Here, c evaluates to 3, as expected.
a = 2;
// However, at this stage, c evaluates to 4 !

Another caveat is that is sometimes more efficient to evaluate the result of an operation once and use the result of this evaluation rather than performing the evaluation of the operation several times.

To avoid those two caveats, the eval function allows the evaluation of an expression, as follows:

const auto c = eval(a + b);

## 1.7 Views

Views allows to map memory area to mathematical objets.

Typical usage of views in MFront

A typical usage of views is given by the example of the integration of a behaviour using an implicit scheme. In such a scheme, a non linear solver tries to determine the increments $$\Delta\,Y$$ of a set of internal state variables packed in a vector $$Y$$ which are the zero of residual denoted $$F$$. The derivative of the residual $$F$$ with respect to $$\Delta\,Y$$ is called the jacobian and is denoted $$J$$.

If one considers a simple plastic law with isotropic hardening, the vector of internal state variables typically contains the elastic strain, a symmetric tensor, and the equivalent plastic strain, a scalar.

In the general case, the vector of internal state variables, the residual and the jacobian can be decomposed as follows:

$Y= \begin{pmatrix} y_{1} \\ \vdots \\ y_{i} \\ \vdots \\ y_{n} \\ \end{pmatrix} \quad\quad F= \begin{pmatrix} f_{y_{1}} \\ \vdots \\ f_{y_{i}} \\ \vdots \\ f_{y_{n}} \\ \end{pmatrix} \quad\quad J = \frac\partial F}\partial \Delta\,Y}} = \begin{pmatrix} \frac\partial f_{y_{1}}}\partial \Delta\,y_{1}}} & \ldots & \ldots & \ldots & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \frac\partial f_{y_{i}}}\partial \Delta\,y_{j}}} & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & \frac\partial f_{y_{N}}}\partial \Delta\,y_{N}}} \\ \end{pmatrix}$

MFront automatically define views to the objects $$y_{i}$$, $$f_{y_{i}}$$, $$\frac\partial f_{y_{i}}}\partial \Delta\,y_{j}}$$ which allows to compute the residual blocks and the jacobian blocks in a very intuitive ways using tensorial operations. Hence, the user never manipulate directly the vectors $$Y$$, $$\Delta\,Y$$ and $$F$$ nor the jacobian matrix $$J$$ but views which acts as tensorial objects.

### 1.7.1 The map function

The map function is a small utility function which simplifies the creation of views from either raw pointers or from tiny vectors (i.e. objects of type tvector).

For example, Listing 1 shows how a vector containing the elastic strain and the equivalent plastic strain can be decomposed by blocks. This decomposition is illustrated on Figure 3.

Listing 1: Usage of the map function

auto Y = tvector<7, double>{0};
auto eel = map<stensor<3u, double>>(Y);
auto& p = Y; Figure 3: “Exemple of decomposition by blocks of a memory area using views”

The map function allows to define offset at compile-time: this allows to checks at compile-time that the memory area is large enough to store the mapped object when mapping a fixed size object (i.e. an object whose size is known at compile-time) from a memory area hold by a tiny vector.

### 1.7.2 The map_array function

The map_array returns an object which acts like a fixed size of mathematical objects. It takes one template argument which describes an arry of of mathematical objects. This template argument is used to determine the number of object mapped and the kind of object mapped.

It can be used as follows:

auto a = map_array<tvector<2u, stensor<2u, double>>>(ptr);

where ptr is a pointer to a suitable memory location.

### 1.7.3 The map_derivative function

The map_derivative function allows to create a view of the derivative of two math objects in a matrix as illustrated by the following example which create a view of the object resulting of the derivation of a symmetric stress tensor with respect to a scalar whose first element is located in element $$1,1$$ of the matrix:

auto r = tmatrix<5, 3>{0};
map_derivative<1, 1, stensor<2u, qt<Stress>>, double>(r) = stensor<2u, Stress>::Id();

The result of this operation is the matrix:

$\begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}$

# 2 Solvers for fixed size non linear systems

The TFEL/Math library provides several non linear solvers which are mainly used by MFront implicit schemes. Those solvers mostly targets systems of small size and focuse on robustness and flexibility. Each solver implements a classical algorithm by also provides several customisation points which allows to create many variants of this basic algorithms.

The solvers available are based on the following algorithms:

• The standard Newton-Raphson algorithm is implemented in the TinyNewtonRaphsonSolver class. This algorithm is coupled with a Powell’ dog leg step in the TinyPowellDogLegNewtonRaphsonSolver class.
• The first and second Broyden algorithms, also known respectively as the good and bad Broyden algorithms are respectively implemented in the TinyBroydenSolver and TinyBroyden2Solver classes. The first Broyden algorithm is coupled with a Powell’ dog leg step in the TinyPowellDogLegBroydenSolver class.
• The Levenberg-Marquardt algorithm is implemented in the TinyLevenbergMarquardtSolver class.

Those classes implements the curiously recurring template pattern (CRTP) to avoid the use of virtual calls . The derived class must provide a method called computeResidual which must compute the residual for the current estimate of the unknowns and, if required by the solver, the jacobian.

A typical usage of those classes is given by the following example:

struct NewtonRaphsonSolver
: public tfel::math::
TinyNewtonRaphsonSolver<2u, double, NewtonRaphsonSolver> {
NewtonRaphsonSolver() {
this->zeros = {0., 0.};
this->epsilon = 1.e-14;
this->iterMax = 20;
}

bool solve() { return this->solveNonLinearSystem(); }

auto getCurrentEstimate() const noexcept { return this->zeros; }

bool computeResidual() noexcept {
constexpr double a = 1.;
constexpr double b = 10.;
auto& x = this->zeros;
auto& f = this->fzeros;
auto& J = this->jacobian;
f = {1 - x(0), 100 * (x(1) - x(0) * x(0))};
J = {-1., 0.,  //
-200 * x(0), 100.};
return true;
}  // end of computeResidual

};  // end of struct NewtonRaphsonSolver

which solves the non linear system:

$\vec{f}{\left(x,y\right)} = \begin{pmatrix} 0 \\ 0 \\ \end{pmatrix} \quad\text{with}\quad \vec{f}{\left(x,y\right)} = \begin{pmatrix} 1-x \\ 10\,{\left(y-x^{2}\right)} \\ \end{pmatrix}$

whose obvious root is $${\left(1,1\right)}$$.

This previous example shows that:

1. The current estimate of the unknowns is stored in a data member called zeros.
2. The residual is stored in a data member called fzeros.
3. The jacobian is stored in a data member called jacobian.

Those names have been chosen in the early versions of MFront and are kept for back-ward compatibility.

This section describes the generic framework used to implement those algorithms.

## 2.1 The TinyNonLinearSolverBase class

All the available solvers are based on the TinyNonLinearSolverBase which provides two main methods solveNonLinearSystem and solveNonLinearSystem2 for the end users and many methods that can be overloaded to customize the behaviour of the algorithm.

In pratice, the NewtonRaphsonSolver class presented in the previous example inherits from TinyNewtonRaphson Solver<2u, double, NewtonRaphsonSolver> which itself inherits from TinyNonLinearSolverBase<2u, double, NewtonRaphsonSolver>.

Note

For the sake of clarity, one template parameter of the TinyNewtonRaphsonSolver and TinyNonLinearSolverBase class have been omitted.

This template parameter describes a data structure containing the so-called workspace of the solver, i.e. all the data members required by the solver. The default value of this template parameter allocates those data members on the stack.

While a bit contrieved, this design ensures that all algorithms share the same customization points, a constraint which is not easy to enforce when relying on CRTP (compared to standard approach based on virtual calls).

As a consequence, the TinyNewtonRaphsonSolver class only implements some methods specific to the Newton-Raphson algorithm while the global algorithmic structure in handled by the TinyNonLinearSolverBase class.

More precisely, the TinyNonLinearSolverBase provides a method called solveNonLinearSystem. This method internally calls a method called solveNonLinearSystem2 which indeed implements the core of the resolution algorithm. The solveNonLinearSystem2 method is described in depth in Section 2.1.2.

The point of the solveNonLinearSystem is to handle failures of the resolution algorithms. The current strategy to handle those failures is described in depth in Section 2.1.3.

### 2.1.1 Some data members of the TinyNonLinearSolverBase class

The TinyNonLinearSolverBase class provides the following data members:

• fzeros: residual vector.
• zeros: current estimate of the unknowns.
• delta_zeros: last computed correction.
• iter: current iteration number.
• is_delta_zeros_defined: boolean stating if an increment of the unknowns has already been computed.
• epsilon: criterion value
• iterMax: maximum number of iterations.

It is worth metionning that a few variables must be initalized, by the base class before calling this method, such as:

• The initial estimate of the unknowns. The method processNewEstimate must be called must be called one this initial estimate is known.
• The values epsilon and iterMax, as well as other numerical parameters required by the core algorithm and declared in derived class.

The iter and has_delta_zeros members are automatically initialized at the beginning of the solveNonLinearSystem method. This method also calls the processNewEstimate method.

### 2.1.2 The solveNonLinearSystem2 method

This method is called internally by the solveNonLinearSystem method. It could be called directly if the required initialization are performed beforehand.

The algorithm implemented by solveNonLinearSystem2 method is depicted in Figure 4. Figure 4: “Flowchart for the resolution of non linear systems proposed by the solveNonLinearSystem2 of the TinyNonLinearSolverBase class”

The computeNewCorrection is the only method that must be implemented in the derived class. It totally defines the resolution algorithm. This method thus has no default implementation.

This flowchart also shows that the solveNonLinearSystem method has many customization points which either defines the underlying algorithm and or that can be use to enhance directly affects the performance and the robustness of the algorithm:

• executeInitialisationTaskBeforeBeginningOfCoreAlgorithm which is started at the start of the method. The default implementation does nothing. This method is mostly meant to be used by the derived class to initialize members specific to the resolution algorithm.
• computeResidual: is meant to compute the residual for the current estimate of the unknowns. This method may also compute the jacobian if required.
• computeResidualNorm which must computes the norm of the residual. By default, this method computes the euclidian norm of the residual. This method.
• checkConvergence which checks if convergence is reached. By default, this method checks if the residual norm is lower than the criterion value epsilon. This method can also be used to implement active-sets algorithms for non-smooth functions.
• processNewCorrection: this method meant to set bounds on some components of the current correction or to implement a line search along the search direction returned by the computeNewCorrection method.
• processNewEstimate: this method is meant to update variables dependant on the current estimate of the unknowns.
• By convention, the computeNewCorrection method may use the following helper methods:
• updateOrCheckJacobian: This method can be used to compute the jacobian or part of the jacobian numerically. If the jacobian was computed in computeResidual, this method can be used to compare it to a numerical approximation.
• solveLinearSystem: This method can be used to solve a linear system if required. By default, an LU method with partial decomposition is used (implemented by the TinyMatrixSolve class).

The algorithm also provides methods which are meant to display informations about the state of the resolution:

• reportInvalidResidualEvaluation which is called method when the evaluation of the residual failed or when the norm of the residual is not finite.
• reportStandardIteration: which is called once the residual is known.

Note that other methods for reporting the current status of the algorithm, such as reportFailure and reportSuccess are also available but are called in the solveNonLinearMethod.

Usage in MFront

It is worth hilighting how the implicit DSL’s provided by MFront may override those methods:

• The computeResidual method first computes the thermodynamic forces using the code provided by the @ComputeThermodynamicForces code bock (or @ComputeStress for mechanical behaviours) and then evaluates the residual using the code given by the @Integrator code block and the jacobian if it is computed analytically.
• The updateOrCheckJacobian method is overloaded if:
• A numerical approximation of the jacobian or blocks of the jacobian are required.
• A comparison of the analytical jacobian to a numerical approximation is required.
• The checkConvergence method is overloaded to take into account additional convergence checks defined in the @AdditionalConvergenceChecks code block. This method can be used to implement active-sets algorithms or interior points methods for multi-surface plasticity or the fixed-point method used by the StandardElastoViscoplasticity brick for porous viscoplasticity.
• The processNewCorrection is used either to set bounds on some components of the current correction based on the physical bounds of the variables or to limit the magnitude of the corrections of some integration variables. See the setMaximumIncrementValuePerIteration method for integration variables or the @MaximumIncrementValuePerIteration keyword.
• The processNewEstimate method is overloaded if some quantities depend on the current estimation of the integration variables. A typical example is given by the case of porous viscoplastic behaviours for which the elastic properties depend on the porosity.

### 2.1.3 The solveNonLinearSystem method

As highlighted by Figure 4, the solveNonLinearSystem2 methods may fail for many reasons. A very common one is that the current estimate of the unknowns are unphysical, leading to a failure in the evaluation of the residual.

The solveNonLinearSystem method implements a simple algorithm which can be seen as a hand crafted line search method which greatly improve the robustness of the non linear solvers. When used correctly, this method may also be used to increase the performances of the non linear solvers (i.e. reduce the total number of iterations). “Flowchart for the resolution of non linear systems proposed by the solveNonLinearSystem of the TinyNonLinearSolverBase class”

The idea of this hand-crafted linesearch is simply to take the last correction to the unknowns and divide it by two and restart the core algorithm. In other words, the search direction is leaved unchanged, but the norm of the correction is reduced. This operation can be repeated several times to find a suitable estimate of the unknowns.

Of course, this only works a correction is known, i.e. if the has_delta_zeros flag is true. Otherwise, this means that the initial guess of the unknowns is incorrect. In this case, we just divide this initial guess by zero. This choice may seem arbitrary but does makes sense in most cases of MFront implicit schemes where the unknowns are almost always the increment of the state variables: if we divide the increments of the state variables, their estimates at the middle of the time step tends to theirs values at the beginning of the time step, which is generally physically acceptable1.

Pratical importance

While very simple, the strategy described in this section is in practice extremly powerful and can be used to easily build very robust and efficient algorithms based on physical considerations. Let us consider a few example:

• If one considers a plastic behaviour, one may want to reject a prediction of the stress which too largely exceeds the current estimate of the plastic limit.
• If one considers a viscoplastic behaviour, one may want to reject a prediction of the stress which would lead to an unphysical viscoplastic strain rate. This is particularly useful with viscoplastic behaviour having an exponential sensitivity to the stress.
• If one considers a plasic behaviour yield surface with sharp corners, such as the approximations of the Mohr-Coulomb criterion given by Abbo and Sloan [4, 5], the solver may fail because the flow direction may oscillate between too iterations. Too solve this issue, we can simply reject iterations that lead to a large change in the flow direction.

The reader may find an example of such algorithms in the case of a perfect plastic behaviour based on the Hosford stress criterion on this page of the MFront gallery: http://tfel.sourceforge.net/hosford.html

An important caveat to this strategy is link to the use of an active-set method to describe multi-surface plasticity. The active-set method makes a priori assumptions on which plastic mechanisms are active and solve the non-linear equations with those assumptions. After convergence, those assumptions are checked and some mechanisms may be activated or desactivated and the non linear solver is restarted. The activation of a plastic mechanism can lead the computeResidual method to fail (as described earlier, rejecting steps leading to a prediction well beyond the current plastic limit is generally a good strategy). The trouble here is that the last correction computed by the solver is very small since the algorithm had converged. Thus, the strategy implemented by the solveNonLinearSystem method would divide a correction that is already almost null. To avoid this caveat, MFront automatically resets the is_delta_zeros_defined member to false. More precisely, here is the implementation of the checkConvergence method when the user has defined additional convergence checks:

bool checkConvergence(const NumericType error) {
auto converged = error < this->epsilon;
auto mfront_internals_converged = converged;
if((mfront_internals_converged) && (!converged)){
this->is_delta_zeros_defined = false;
}
}
1.
Abrahams, David and Gurtovoy, Aleksey. C++ template metaprogramming: Concepts, tools, and techniques from boost and beyond. Boston : Addison-Welsley, 2004. ISBN 0321227255 9780321227256.
2.
Veldhuizen, Todd. Techniques for Scientific C++. 1999.
3.
Coplien, James O. Curiously recurring template patterns. C++ Report. 1995. Vol. 7, no. 2, p. 24–27.
4.
Abbo, A. J. and Sloan, S. W. A smooth hyperbolic approximation to the Mohr-Coulomb yield criterion. Computers & Structures. February 1995. Vol. 54, no. 3, p. 427–441. DOI 10.1016/0045-7949(94)00339-5. Available from: http://linkinghub.elsevier.com/retrieve/pii/0045794994003395
5.
Nagel, Thomas, Minkley, Wolfgang, Böttcher, Norbert, Naumov, Dmitri, Görke, Uwe-Jens and Kolditz, Olaf. Implicit numerical integration and consistent linearization of inelastic constitutive models of rock salt. Computers & Structures. April 2017. Vol. 182, p. 87–103. DOI 10.1016/j.compstruc.2016.11.010. Available from: http://www.scopus.com/inward/record.url?eid=2-s2.0-85006482432&partnerID=MN8TOARS http://linkinghub.elsevier.com/retrieve/pii/S0045794916306319

1. This is the reason why, by default, the initial guess of the unknowns in MFront is simply a null vector. The user may specify an initial guess for the unknowns using the @Predictor code block, although this is seldom used.↩︎