[Dune-devel] Tensor data structures and operations on tensors

[Simon Praetorius]

Dear all,

Recently I came across the problem that I wanted to compute the second 
derivative of a vector-valued function, or even higher-order 
derivatives. And I was wondering which data-structure to use as the 
range type of such an operation. Is it a 
FieldMatrix<FieldVector<T,m>,n,o> or better a 
FieldVector<FieldMatrix<T,m,n>,o> and what about the higher-order 
tensors. Is it even "allowed" to nest FieldMatrix/FieldVector in such a 
way? The value-type in these cases is not a field-type anymore, but a 
vector or matrix. And how to multiply these nested matrices with vectors 
or other matrices and make coordinate transformations and so on...?

This was an initial motivation to think about a nice abstraction of 
tensors. Tensors that maybe could replace the current FieldMatrix and 
FieldVector implementations and even others of our Vector and Matrix 
classes. We have a lot of code duplication, although we already use the 
CRTP base-class DenseVector/-Matrix. And recently some effort had been 
made to get products of matrices of mixed type work in various combinations.

So, I came up with an implementation that follows the c++ standard 
proposals for mdarray and mdspan in its interface but additionally 
provides numerical operations, like vector-space operations and some 
tensor contractions like in the numpy or tensorflow framework.

See https://gitlab.dune-project.org/simon.praetorius/dune-tensor

The README gives an overview about the data structures and algorithms 
provided. The implementation is very flexible and allows to define your 
own tensor types, in terms of storage backend or index mapping. Five 
examples are given: FieldTensor, DynamicTensor, FixedSizeSparseTensor 
(with CSF storage), DynamicSizeSparseTensor, and a DiagonalTensor. These 
are just experiments but the FieldTensor and DynamicTensor provide 
already a lot of the functionality of the DenseVector and DenseMatrix 
classes in dune-common.

One main difference there is, though. Instead of interpreting a 
FieldVector<T,1> and FieldMatrix<T,1,1> as a scalar, the rank-0 tensor 
could be interpreted as a scalar, while the tensors of higher-rank with 
size=1 still are non-scalar tensors. This could solve some issue we had, 
especially in the 1d simulations where it was unclear whether a vector 
of size=1 is a scalar or a vector.

If you are interested, it would be nice to get some feedback. The 
implementation is currently not cleaned up, but already includes some 
documentation. It also contains some benchmarking sources.

Best regards,
Simon

-- 
Dr. Simon Praetorius
Technische Universität Dresden
Institute of Scientific Computing
phone: +49 351 463-34432
mail: simon.praetorius@​tu-dresden.de