TFEL/Math
libraryThis overview of the TFEL/Math
library mostly focuses on two important features:
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.TFEL/Math
library provides efficient and robust non linear solvers for small sized problems described in Section 2.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:
tvector
class)stensor
class)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.
Quantities were introduced in TFEL/Math
following ideas of the book D. Abrahams and A. Gurstovoy [1] 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.
qt
classA 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.
qt_ref
and const_qt_ref
classesThe 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;
}();std::abs(value - 200.) < eps); TFEL_TESTS_STATIC_ASSERT(
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.
Quantity
classFor 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.
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
4; m2 *=
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
3u,qt<Stress>> s;
stensor<// 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
functionThe
power
function is a very convenient function which takes two template parametersN
andD
defining the exponent as the fraction \(N/D\). The default value ofD
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, theNth
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 thestd::sqrt
function.- \(3\): In this case, the square root of the
Nth
power of the argument is computed using thestd::cbrt
function.
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.
ArrayPolicy
conceptArray 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:
value_type
which corresponds to the logical type hold by the array (see below).const_value_type
which corresponds to a constant value to a an object of the logical type hold by the array.storage_type
which is the type actually stored in memory.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&)
other_type
denotes for each operator a type for which the considered operator in the value_type
class is valid.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.make_const_reference
which turns a a constant reference to a storage_type
object into a const_reference
object.make_reference
which turns a reference to a storage_type
object into a reference
object.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:
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
orlong double
).
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.
IndexingPolicy
conceptA class matching the IndexingPolicy
concept must provide:
size_type
to the type used to index data.RowMajorIndexingPolicy
which can be used to import data from raw C
pointers.arity
of type size_type
which describes the number of indices required to describe the object.hasFixedSizes
. If this data member is true
, the class must be stateless, i.e. empty.areDataContiguous
which states if the data are continoussize
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.size
which takes a size_type
argument returning range in the given dimension. This member function must be constexpr
if hasFixedSizes
is true.getUnderlyingArrayMinimalSize
which returns the minimal size of an array required by the indexing policy. This member function must be constexpr
if hasFixedSizes
is true.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.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.
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.
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} \]
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.
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:
1
, 2
or 3
).tmatrix
caseA fixed size matrix, implemented by the tmatrix
class can be seen as the derivative of a tiny tensor with respect to a tiny vecor.
derivative_type
metafunctionThe 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>>
.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 [2].
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]; }
About loop unrolling
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:2 * b; c = a +
is equivalent to:
0] = a[0] + 2 * b[0]; c[1] = a[1] + 2 * b[1]; c[2] = a[2] + 2 * b[2]; c[
eval
functionOne possible caveat of the expression template technique can be illustrated by the following example. Let us consider two vectors a
and b
:
0] = 1;
a[0] = 2;
b[auto c = a + b;
// Here, c[0] evaluates to 3, as expected.
0] = 2;
a[// However, at this stage, c[0] 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);
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 = {\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,Y}} = \begin{pmatrix} {\displaystyle \frac{\displaystyle \partial f_{y_{1}}}{\displaystyle \partial \Delta\,y_{1}}} & \ldots & \ldots & \ldots & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & {\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\,y_{j}}} & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \ldots & \ldots & \ldots & \ldots & {\displaystyle \frac{\displaystyle \partial f_{y_{N}}}{\displaystyle \partial \Delta\,y_{N}}} \\ \end{pmatrix} \]
MFront
automatically define views to the objects \(y_{i}\), \(f_{y_{i}}\), \({\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \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.
map
functionThe 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[6];
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.
map_array
functionThe 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.
map_derivative
functionThe 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};
1, 1, stensor<2u, qt<Stress>>, double>(r) = stensor<2u, Stress>::Id(); map_derivative<
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} \]
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:
TinyNewtonRaphsonSolver
class. This algorithm is coupled with a Powell’ dog leg step in the TinyPowellDogLegNewtonRaphsonSolver
class.TinyBroydenSolver
and TinyBroyden2Solver
classes. The first Broyden algorithm is coupled with a Powell’ dog leg step in the TinyPowellDogLegBroydenSolver
class.TinyLevenbergMarquardtSolver
class.Those classes implements the curiously recurring template pattern (CRTP) to avoid the use of virtual calls [3]. 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::
: 2u, double, NewtonRaphsonSolver> {
TinyNewtonRaphsonSolver<
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;
1 - x(0), 100 * (x(1) - x(0) * x(0))};
f = {1., 0., //
J = {-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:
zeros
.fzeros
.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.
TinyNonLinearSolverBase
classAll 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
andTinyNonLinearSolverBase
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.
TinyNonLinearSolverBase
classThe 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 valueiterMax
: 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:
processNewEstimate
must be called must be called one this initial estimate is known.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.
solveNonLinearSystem2
methodThis 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.
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.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 theStandardElastoViscoplasticity
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 thesetMaximumIncrementValuePerIteration
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.
solveNonLinearSystem
methodAs 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).
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;
this->additionalConvergenceChecks(converged, error);
if((mfront_internals_converged) && (!converged)){
this->is_delta_zeros_defined = false;
} }
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.↩︎