Tensors.jl — Tensor Computations in Julia

Tensors.jl is a Julia package that provides efficient computations with symmetric and non-symmetric tensors. The focus is on the kind of tensors commonly used in e.g. continuum mechanics and fluid dynamics. Exploiting Julia’s ability to overload Unicode infix operators and using Unicode in identifiers, implemented tensor expressions commonly look very similar to their mathematical writing. This possibly reduces the number of bugs in implementations. Operations on tensors are often compiled into the minimum assembly instructions required, and, when beneficial, SIMD-instructions are used. Computations involving symmetric tensors take symmetry into account to reduce computational cost. Automatic differentiation is supported, which means that most functions written in pure Julia can be efficiently differentiated without having to implement the derivative by hand. The package is useful in applications where efficient tensor operations are required, e.g. in the Finite Element Method.


Introduction
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 Voigt notation or Voigt format, described in many classic FEM textbooks [1,3,6]. In this format, second-order tensors are stored in vectors 1 of length 2 dim n or n dim (n dim + 1)/2 for the non-symmetric and symmetric case, respectively, where n dim denotes the number of spatial dimensions for which the problem is formulated. Examples of second-order tensors are the symmetric Cauchy stress σ and the nonsymmetric Kirchoff stress τ T symmetric, 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 C and a second-order tensor ε can be formulated as a matrix-vector multiplication, and the open product between two second-order tensors as a column vector times a row vector, viz.
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 σ and ε are symmetric second-order tensors, it is common to define where γ ij is usually denoted "engineering strains". 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 component of a rank-2 tensor σ would look like σ [6], or σ[idx [1,2]]] where idx is a lookup table matching the "Cartesian index" [1,2] to the "linear index" [6] which is the location of σ xy in the Voigt vector. • Some operations will silently give wrong result when naively applied to a tensor in Voigt notation. One example is the norm function ij ij σ σ = σ . Using the norm function for vectors will give the correct answer if the tensor is stored as an non-symmetric tensor. However, it will silently give the wrong result if the tensor is stored as a symmetric tensor (since the offdiagonal components will only be accounted for once).
• 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 strain-like tensors are commonly stored with a factor two on the off-diagonals. This is a frequent source of confusion because the same mathematical object is stored differently.
Tensors.jl was created in an attempt to overcome the many deficiencies of Voigt notation. It is written in the programming language Julia [2].
The main design goals have been: 1. Performance -Operations should compile to (close to) the bare minimum of assembly operations needed. Symmetry should be exploited for computational efficiency. SIMD-instructions should be used when computationally beneficial. 2. Generality -The same code should work regardless if the number types are double or single precision (denoted Float64 and Float32 in Julia) or even user defined numerical types, e.g. dual numbers used in forward mode automatic differentiation. It should also be dimension independent so that the same code can be used for one, two and three spatial dimensions.
3. Clarity -Implementation details, such as the particular way a tensor is stored, should not be visible to the user. Symmetric and non-symmetric tensors should behave the same, with the difference that operations on symmetric tensors should, if possible, be faster. Operations should be visually similar to mathematical writing. This includes using Unicode infix operators such as ⊗ for the open product and · for the scalar product.
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. We note that these operations all get compiled into very efficient machine code, and use SIMD whenever it is beneficial. In Table 1 we show a non-exhaustive summary of tensor operations that are currently implemented.

Illustrative usage example
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 F = I + ∇ ⊗ u, where u is the displacement from the reference to the current configuration, the right Cauchy-Green deformation tensor is defined by C = F T · F. The Second Piola-Kirchoff stress tensor S is derived from the Helmholtz free energy ψ by the relation . The free energy for a hyperelastic material can be defined as where Cˆ = det(C) -1/3 C and The analytic expression for the Second Piola-Kirchoff tensor is which can be implemented by the following function:

Automatic Differentiation
Automatic differentiation [4] (AD) is a numerical method for differentiating functions implemented in a programming language. AD has many advantages over other numerical methods for differentiation. Comparing it to numerical differentiation, where components are perturbed to compute a gradient, AD does not suffer from cancellation and can compute multiple partial derivatives in a single function call. Tensors.jl supports AD and as an example of its use, we here recall the function 2 ∂Ψ ∂ = C S , compute it using AD, and compare with the analytical result: julia> H = rand(Tensor{2,3}); F = one(H) + H; C = symmetric(F' · F); julia> 2 * gradient(C -> ψ (C, μ, Kb), C) 3×3 Tensors.SymmetricTensor{2, 3 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 dual numbers defined in ForwardDiff.jl [5].

Performance
In this section we compare the performance of a selected number of operations when using the Tensor types to the Voigt format, implemented with standard Array types. 4 The results are shown in Table 2, where u denotes a vector, A, A sym denote second-order non-symmetric and symmetric tensors and C, C sym denote fourth-order non-symmetric and symmetric tensors, respectively. All tensors presented are in three dimensions. The Table 1: Summary of implemented tensor operations. u, v denotes vectors, A, B denotes second-order symmetric or non-symmetric tensors, and C, D denotes fourthorder symmetric or non-symmetric tensors. We note that instead of using : for infix double contraction we use (written as \boxdot). This is because : does not have the same operator precedence as multiplication in Julia.

Julia code infix
Single contraction

Correctness and performance testing
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 Arrays. Code coverage for the package is currently at 95%. The examples in the documentation are written as "doctests" which means that the code in the examples are running as part of CI. This is to prevent examples in the documentation from becoming stale. A comprehensive set of benchmarks are implemented and used to check that regressions in performance are not introduced.

Implementation details
The elements of a tensor are stored internally in a Tuple. In Julia, a Tuple is an immutable container, where the length is part of the type. Consequently, when compiling a function, the Julia compiler can statically know the length of the tuple-container and use that information for optimizations. Furthermore, tuples can be stored on the stack which removes any use of expensive heap allocations, alleviating the need to preallocate output buffers.
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.

Conclusion
We have presented the Julia package Tensors.jl. It provides a framework for doing computations with nonsymmetric and symmetric tensors of rank-1, rank-2 and rank-4 with arbitrary number types. The implementation has many advantages over using the common Voigt format, such as, higher performance, dimension generality, and allows for a more direct mapping between the source code and the mathematical notation. Automatic differentiation for first and second-order derivatives is supported and is implemented efficiently.

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

Programming language
Julia 0.6
These dependencies are automatically installed upon installing via Pkg.add("Tensors").