Partial Differential Equations (PDEs) describing natural phenomena are modelled using tensors of different order. Two commonly studied problems are heat transfer, which include temperature and heat flux (rank-0 and rank-1 tensor, respectively), and continuum mechanics, which include stress and strain (rank-2 tensors) and the so-called tangent stiffness (rank-4 tensor). These tensors are often symmetric but may also be non-symmetric, for example in finite deformation continuum mechanics.

In the implementation of numerical solution schemes for such PDEs, a large number of expressions involving tensors typically have to be computed. As an example, in the Finite Element Method (FEM), the weak form of the PDE to be solved is evaluated multiple times in every element of the mesh. The total number of elements in a model can easily exceed millions, and it is therefore desirable to have access to a library that can perform these tensor operations efficiently. Another aspect to consider is the level of difficulty to implement tensor expressions as computer source code based on their mathematical form. A close correspondence between the source code and the mathematical writing will likely lead to a quicker, less error prone, implementation process.

The classical way of treating tensors in computational mechanics is to use what is commonly called _{dim} (_{dim} + 1)/2 for the non-symmetric and symmetric case, respectively, where _{dim} denotes the number of spatial dimensions for which the problem is formulated. Examples of second-order tensors are the symmetric Cauchy stress

which are here represented in Voigt notation. Similarly, fourth-order tensors are stored in square matrices where the number of columns equals the length of the vector for the second-order tensors in Voigt notation.

The purpose of Voigt notation is that some standard linear algebra operators can be used for common tensor operations such as open products, dot products and double contractions. As an example, the double contraction between a fourth-order tensor

The left column is presented using the Einstein summation convention and the right column is presented with the corresponding, standard matrix operations.

In order to preserve the relation _{ij}ε_{ij}^{T}

Where γ_{ij}

While Voigt notation is simple to adopt when using programming languages or libraries that provide linear algebra functionality, it does, however, suffer from a few drawbacks regarding both performance and clarity:

Some operations, such as the scalar product between a rank-2 and rank-1 tensor, become difficult since the rank-2 tensor is not stored as a matrix.

Indexing into a tensor stored in Voigt format is not straightforward. Fetching the _{xy}

Some operations will silently give wrong result when naively applied to a tensor in Voigt notation. One example is the norm function

Since tensors in Voigt format are usually stored in arrays that allows arbitrary number of elements, a compiler cannot know the size of the array at compile time. This prohibits optimal code to be generated.

In order for certain operations to work as “expected” for symmetric tensors e.g. the double contraction _{ij}ε_{ij}^{T}

The main design goals have been:

The main purpose of the software is to be used for solving PDEs modelling physical phenomena (such as heat flux and stress equilibrium). In particular this excludes rectangular tensors, and tensors of dimension higher than 3, which are commonly used in e.g. machine learning. Currently only rank-1, rank-2 and rank-4 tensors, in up to 3 dimensions, are implemented. We note that support for rank-3 tensors fall within the scope of the package, but is not yet implemented.

The foundation of the package are the two types

Operations on these tensors can now be performed. When there is a corresponding mathematical symbol for the operation, it can often be used directly:

We note that these operations all get compiled into very efficient machine code, and use SIMD whenever it is beneficial. In Table

Summary of implemented tensor operations.

Operation | Julia code | infix |
---|---|---|

Single contraction | ||

_{i}v_{i} |
u · v | |

_{ij}v_{j} |
A · v | |

_{ij}B_{jk} |
A · B | |

Double contraction | ||

A ⊡ B | ||

_{ijkl}B_{kl} |
C ⊡ B | |

_{ijkl}D_{klmn} |
C ⊡ D | |

Outer product | ||

_{i}v_{j} |
u ⊗ v | |

_{ij}B_{kl} |
A ⊗ B | |

Other operations | ||

det( |
||

inv( |
||

norm( |
||

^{T} |
||

As an illustrative example of using the package, we here give the mathematical formulation of an energy potential and stress in large deformation continuum mechanics and its implementation as a function in Julia.

For a deformation gradient ^{T} ·

The free energy for a hyperelastic material can be defined as

where ^{–1/3}_{b}

The analytic expression for the Second Piola-Kirchoff tensor is

which can be implemented by the following function

These functions work in 1, 2 and 3 dimensions and have good performance.

Automatic Differentiation [

The slowdown from using AD instead of the analytic version is about a factor of 3. It is also possible to compute second-order derivatives exactly in an analogous manner as the example for first order derivatives above.

The AD-functionality is built upon the

In this section we compare the performance of a selected number of operations when using the ^{sym} denote second-order non-symmetric and symmetric tensors and ^{sym} denote fourth-order non-symmetric and symmetric tensors, respectively. All tensors presented are in three dimensions. The benchmarks were performed using the benchmarking tool

Comparison of performance for some tensor operations using

Operation | Speed-up | ||
---|---|---|---|

Single contraction | |||

1.241 ns | 9.795 ns | ×7.9 | |

2.161 ns | 58.769 ns | ×27.2 | |

3.117 ns | 44.395 ns | ×14.2 | |

^{sym} ^{sym} |
5.125 ns | 44.498 ns | ×8.7 |

Double contraction | |||

1.927 ns | 12.189 ns | ×6.3 | |

^{sym} ^{sym} |
1.927 ns | 12.187 ns | ×6.3 |

6.087 ns | 78.554 ns | ×12.9 | |

60.820 ns | 280.502 ns | ×4.6 | |

^{sym} ^{sym} |
22.104 ns | 281.003 ns | ×12.7 |

^{sym} ^{sym} : ^{sym} |
9.466 ns | 89.747 ns | ×9.5 |

Outer product | |||

2.167 ns | 32.447 ns | ×15.0 | |

9.801 ns | 6.568 ns | ×8.8 | |

Other operations | |||

det( |
1.924 ns | 177.134 ns | ×92.1 |

inv(^{sym}) |
4.587 ns | 635.858 ns | ×138.6 |

norm( |
1.990 ns | 16.752 ns | ×8.4 |

The package is tested for correctness using Continuous Integration (CI) on macOS, Linux and Windows versions of Julia. An extensive test suite based on unit testing is used. The testing includes tensor identities, and tests of different tensor operations. The results are compared with the matrix/vector representation of the tensors using standard Julia

The elements of a tensor are stored internally in a

For symmetric tensors only the “lower half” of the tensor is stored. A non-symmetric fourth-order tensors in 3 dimensions has 81 elements, while the symmetric version only need to store 36 elements.

Operations on tensors are frequently implemented using, what in Julia are called, “generated functions”. This allows the programmer to hook into the compilation process, at the time where the types of the arguments in the called function are known. Based on these types, the programmer can generate arbitrary Julia code which get compiled just like a normal function. This is, for example, used to generate the SIMD code for different tensor sizes.

Julia will also generate specialized code even when a user defined numeric type is used as the elements of the tensors. This is the case when doing automatic differentiation where a dual-number is used as the element type. This means that well performing code is generated even for automatic differentiation.

Since the number of dimensions and ranks for the tensors that this package support is limited, excessive compilation and code generation is prevented.

We have presented the Julia package

Any OS supported by Julia which include FreeBSD, Windows, macOS, Linux, Raspberry Pi and other ARM systems.

Julia 0.6

The following packages are required to use

These dependencies are automatically installed upon installing via

Kristoffer Carlsson, Fredrik Ekre, Keita Nakamura

English

Since tensor computations is of such fundamental use when modeling a large class of physical systems, a high quality implementation has high reuse potential in other packages modeling these physical systems, where

The term

For technical reasons the total number of elements stored is also a parameter of the type, but this is rarely of interest to a user of the package.

Unicode symbols like σ can be entered in the Julia REPL by entering

Benchmarks were performed with an Intel(R) Core(TM) i5-7600K CPU @ 3.80 GHz CPU.

We would like to thank the Julia community and, especially, Jarett Revels for creating

The authors have no competing interests to declare.