Performance of the structured types

TensND is symbolic and numerical. The structured types (TensISO, TensTI, TensOrtho) exist so that products and inverses are computed on a handful of stored coefficients rather than on   components. This tutorial measures what that buys, and where it stops applying.

Theory: Isotropic tensors, The Walpole basis, Orthotropy.

using TensND
using LinearAlgebra
using BenchmarkTools
using Printf

arr(x) = get_array(x)

arr (generic function with 1 method)

What each type stores

The whole point in one table: the number of scalars held, against the 81 components of a dense order-4 tensor in 3-D.

ℬ = CanonicalBasis{3, Float64}()
n = (0.0, 0.0, 1.0)

iso = TensISO{3}(3 * 20.0, 2 * 8.0)
ti = TensTI{4}(12.0, 13.0, 3.5, 7.0, 4.0, n)
ort = TensOrtho(10.0, 8.0, 9.0, 3.0, 2.0, 4.0, 2.5, 3.0, 1.5, ℬ)
dense = Tens(arr(ort))

for (name, t) in ("TensISO" => iso, "TensTI{4,5}" => ti, "TensOrtho" => ort)
    @printf "  %-12s stores %2d scalars   (dense: 81 components)\n" name length(get_data(t))
end

  TensISO      stores  2 scalars   (dense: 81 components)
  TensTI{4,5}  stores  5 scalars   (dense: 81 components)
  TensOrtho    stores  9 scalars   (dense: 81 components)

Double contraction

Each structured type has a closed-form ⊡ that stays in the stored representation; the dense route expands all 81 components through an einsum.

function bench_dcontract()
    rows = Any[]
    for (name, a, b) in (
            ("TensISO ⊡ TensISO", iso, TensISO{3}(3 * 30.0, 2 * 12.0)),
            ("TensTI ⊡ TensTI", ti, TensTI{4}(10.0, 9.0, 2.0, 5.0, 3.0, n)),
            ("TensOrtho ⊡ TensOrtho", ort, TensOrtho(6.0, 5.0, 7.0, 2.0, 1.0, 3.0, 1.2, 2.2, 0.8, ℬ)),
        )
        ts = @belapsed $a ⊡ $b
        aa, bb = arr(a), arr(b)
        td = @belapsed $aa ⊡ $bb
        push!(rows, (name, ts, td))
    end
    return rows
end

println("  operation                     structured        dense        speed-up")
for (name, ts, td) in bench_dcontract()
    @printf "  %-26s %8.1f ns   %10.1f ns   ×%6.1f\n" name 1.0e9 * ts 1.0e9 * td td / ts
end

  operation                     structured        dense        speed-up
  TensISO ⊡ TensISO               2.8 ns       5593.7 ns   ×2008.5
  TensTI ⊡ TensTI                11.7 ns       5582.0 ns   × 477.4
  TensOrtho ⊡ TensOrtho          87.5 ns       5578.7 ns   ×  63.8

The isotropic case is two scalar products; the Walpole case a   matrix product plus two scalar products; the orthotropic case one   product plus three scalar products.

Inversion

The same story, and a starker one: inverting a dense minor-symmetric order-4 tensor means inverting a   matrix.

function bench_inv()
    rows = Any[]
    for (name, t) in ("TensISO" => iso, "TensTI{4,5}" => ti, "TensOrtho" => ort)
        ts = @belapsed inv($t)
        km = Matrix(KM(t))
        td = @belapsed inv($km)
        push!(rows, (name, ts, td))
    end
    return rows
end

println("  type            structured        6×6 inverse     speed-up")
for (name, ts, td) in bench_inv()
    @printf "  %-14s %8.1f ns   %12.1f ns   ×%6.1f\n" name 1.0e9 * ts 1.0e9 * td td / ts
end

  type            structured        6×6 inverse     speed-up
  TensISO             3.1 ns          908.2 ns   × 294.4
  TensTI{4,5}         3.4 ns          925.2 ns   × 271.6
  TensOrtho          10.0 ns          926.7 ns   ×  92.6

Allocations

A structured operation allocates a single small object — the tuple of stored coefficients — where the dense route allocates intermediate arrays of 81 components. The figure below is the whole cost of the result, and it scales with the number of stored scalars, not with .

for (name, f) in (
        "TensISO ⊡ TensISO" => () -> iso ⊡ iso,
        "inv(TensISO)" => () -> inv(iso),
        "TensTI ⊡ TensTI" => () -> ti ⊡ ti,
        "inv(TensTI)" => () -> inv(ti),
        "inv(TensOrtho)" => () -> inv(ort),
    )
    b = @benchmark $f()
    @printf "  %-22s %6d bytes in %d allocations\n" name b.memory b.allocs
end

  TensISO ⊡ TensISO          32 bytes in 1 allocations
  inv(TensISO)               32 bytes in 1 allocations
  TensTI ⊡ TensTI            80 bytes in 1 allocations
  inv(TensTI)                80 bytes in 1 allocations
  inv(TensOrtho)             80 bytes in 1 allocations

Element access

A structured type has no component array to index into, so getindex reconstructs the component from the stored coefficients. That is cheap, but repeated indexing is not the way to use these types: if a whole array is needed, ask for it once with get_array.

@printf "  getindex on TensOrtho : %6.1f ns\n" 1.0e9 * @belapsed $ort[1, 2, 1, 2]
@printf "  getindex on a dense   : %6.1f ns\n" 1.0e9 * @belapsed $(arr(ort))[1, 2, 1, 2]
@printf "  get_array(TensOrtho)  : %6.1f ns  (once, then index freely)\n" 1.0e9 * @belapsed get_array($ort)

  getindex on TensOrtho :    3.4 ns
  getindex on a dense   :    3.4 ns
  get_array(TensOrtho)  :  471.0 ns  (once, then index freely)

Where the structure stops paying

The closed forms apply only while the result stays in the class. Two cases leave it, and both fall back to the generic route:

• a product of two TensOrtho expressed in different material frames is generally fully anisotropic;

• a product of two orthotropic tensors in the same frame is orthotropic but not major-symmetric — twelve constants where TensOrtho stores nine — so the result is returned as a generic tensor (Orthotropy).

ort2 = TensOrtho(6.0, 5.0, 7.0, 2.0, 1.0, 3.0, 1.2, 2.2, 0.8, ℬ)
ort_rot = TensOrtho(6.0, 5.0, 7.0, 2.0, 1.0, 3.0, 1.2, 2.2, 0.8, Basis(0.3, 0.7, 0.2))

println("same frame      : ", typeof(ort ⊡ ort2))
println("different frames: ", typeof(ort ⊡ ort_rot))

same frame      : TensND.TensCanonical{4, 3, Float64, Tensors.SymmetricTensor{4, 3, Float64, 36}}
different frames: TensND.TensCanonical{4, 3, Float64, Tensors.SymmetricTensor{4, 3, Float64, 36}}

Both agree with the dense computation to machine precision — the fallback is exact, not approximate:

(
    norm(arr(ort ⊡ ort2) - arr(ort) ⊡ arr(ort2)),
    norm(arr(ort ⊡ ort_rot) - arr(ort) ⊡ arr(ort_rot)),
)

(7.377764055609925e-15, 3.595820250601821e-15)

Mixed-class arithmetic

Operations between different classes promote to the least constrained of the two, in closed form where possible (iso_to_ortho, walpole_to_ortho):

for (name, x, y) in (
        ("TensISO ⊡ TensTI", iso, ti),
        ("TensISO ⊡ TensOrtho", iso, ort),
        ("TensTI ⊡ TensISO", ti, iso),
    )
    @printf "  %-22s → %s\n" name string(typeof(x ⊡ y))
end

  TensISO ⊡ TensTI       → TensTI{4, Float64, 6}
  TensISO ⊡ TensOrtho    → TensND.TensCanonical{4, 3, Float64, Tensors.SymmetricTensor{4, 3, Float64, 36}}
  TensTI ⊡ TensISO       → TensTI{4, Float64, 6}

Summary

	stored scalars	⊡	inv
TensISO	2	two scalar products	two reciprocals
TensTI {4,5}	5	product + 2 scalars	inverse + 2 reciprocals
TensOrtho	9	product + 3 scalars	inverse + 3 reciprocals
dense	81	einsum over 81 components	inverse

Use a structured type whenever the physics guarantees the symmetry. Convert to a dense array once, with get_array, when component-by-component access is genuinely needed.

---

This page was generated using Literate.jl.