From d25fca303c3b6c808432d1345781b4196c207873 Mon Sep 17 00:00:00 2001 From: miguelmaso Date: Tue, 22 Sep 2026 22:17:49 +0200 Subject: [PATCH] Zero-allocation tangent for ViscousPolyconvex MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add generated implementations of `-` (binary and unary) for `TensorValue`, mirroring the existing `+`. Without them, subtraction of fourth-order tensors fell back to the generic Gridap implementation, which exceeds Julia's inlining threshold for `TensorValue{9,9}` and heap-allocates the result. Those subtractions appear in `∂Sv∂C_Cᵥfix` and `∂Cv⁻¹∂C`, i.e. exactly on the path of the viscoelastic tangent operator, which is why only `∂∂Ψ∂FF` allocated while `Ψ` and `∂Ψ∂F` did not. `+` is also generalized from `TensorValue{D,D}` to `TensorValue{D1,D2}` so that non-square tensors are covered as well. No change in the mathematics. The three methods build the result with an explicit element type, which also covers empty tensors: Gridap builds `TensorValue{0,3}` values when it writes the vertex grid of a 3D model, and leaving the type implicit made them throw `UndefVarError: T not defined`, breaking `writevtk(model)` for any code that loads HyperFEM. The price of fixing the element type with `promote_type` is that adding two `Bool` tensors now errors instead of promoting to `Int` as Gridap does; the tensors involved in a finite element evaluation are floating point. Visco-polyconvex ∂∂Ψ∂FF: 187 allocs / 15552 B / 5367 ns -> 0 allocs / 0 B / 599 ns. Visco-elastic ∂∂Ψ∂FF: 1021 allocs / 49600 ns -> 501 allocs / 34200 ns. Also interpolate the arguments of the tensor algebra benchmarks. Without `$`, the benchmarked expressions read non-const globals, so the call is dynamically dispatched and the result is boxed: this is the source of the spurious "1 alloc of 0.656 kB" reported for `IIsym`, `push_forward_C_to_F` and `×ᵢ⁴`, which are already allocation-free. Co-Authored-By: Claude Opus 5 --- .../TensorAlgebraBenchmarks.jl | 10 +++---- src/TensorAlgebra/Operations.jl | 26 ++++++++++++++----- src/TensorAlgebra/TensorAlgebra.jl | 2 ++ test/TestTensorAlgebra/TensorAlgebraTests.jl | 9 +++++++ 4 files changed, 36 insertions(+), 11 deletions(-) diff --git a/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl b/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl index d7974b2..9334a52 100644 --- a/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl +++ b/benchmark/TensorAlgebraBenchmarks/TensorAlgebraBenchmarks.jl @@ -49,8 +49,8 @@ H = TensorValue(1.:81...) SUITE["Tensor algebra"]["δδ_μ_2d"] = @benchmarkable δᵢₖδⱼₗ2D + δᵢₗδⱼₖ2D SUITE["Tensor algebra"]["δδ_λ_2d"] = @benchmarkable 1.0 * δᵢⱼδₖₗ2D -SUITE["Tensor algebra"]["Cofactor"] = @benchmarkable cof(A) -SUITE["Tensor algebra"]["Det(A)Inv(A')"] = @benchmarkable det(A)*inv(A') -SUITE["Tensor algebra"]["×ᵢ⁴"] = @benchmarkable ×ᵢ⁴(A) -SUITE["Tensor algebra"]["IIsym"] = @benchmarkable IIsym(A) -SUITE["Tensor algebra"]["push_forward_C_to_F"] = @benchmarkable push_forward_C_to_F(F,H) +SUITE["Tensor algebra"]["Cofactor"] = @benchmarkable cof($A) +SUITE["Tensor algebra"]["Det(A)Inv(A')"] = @benchmarkable det($A)*inv($A') +SUITE["Tensor algebra"]["×ᵢ⁴"] = @benchmarkable ×ᵢ⁴($A) +SUITE["Tensor algebra"]["IIsym"] = @benchmarkable IIsym($A) +SUITE["Tensor algebra"]["push_forward_C_to_F"] = @benchmarkable push_forward_C_to_F($F,$H) diff --git a/src/TensorAlgebra/Operations.jl b/src/TensorAlgebra/Operations.jl index 24a974b..0f36f22 100644 --- a/src/TensorAlgebra/Operations.jl +++ b/src/TensorAlgebra/Operations.jl @@ -10,12 +10,26 @@ function (*)(Ten1::TensorValue, Ten2::TensorValue) end -@inline @generated function (+)(A::TensorValue{D,D}, B::TensorValue{D,D}) where {D} - str = "" - for i in 1:D*D - str *= "A.data[$i] + B.data[$i], " - end - Meta.parse("TensorValue{D,D}($str)") +# The element type is explicit so that empty tensors (e.g. `TensorValue{0,3}`, +# built by Gridap for vertex grids) are also supported. +@inline @generated function (+)(A::TensorValue{D1,D2}, B::TensorValue{D1,D2}) where {D1,D2} + T = promote_type(eltype(A), eltype(B)) + data = [:(A.data[$i] + B.data[$i]) for i in 1:D1*D2] + :(TensorValue{D1,D2,$T}($(Expr(:tuple, data...)))) +end + + +@inline @generated function (-)(A::TensorValue{D1,D2}, B::TensorValue{D1,D2}) where {D1,D2} + T = promote_type(eltype(A), eltype(B)) + data = [:(A.data[$i] - B.data[$i]) for i in 1:D1*D2] + :(TensorValue{D1,D2,$T}($(Expr(:tuple, data...)))) +end + + +@inline @generated function (-)(A::TensorValue{D1,D2}) where {D1,D2} + T = eltype(A) + data = [:(-A.data[$i]) for i in 1:D1*D2] + :(TensorValue{D1,D2,$T}($(Expr(:tuple, data...)))) end diff --git a/src/TensorAlgebra/TensorAlgebra.jl b/src/TensorAlgebra/TensorAlgebra.jl index dfc7839..76abf9a 100644 --- a/src/TensorAlgebra/TensorAlgebra.jl +++ b/src/TensorAlgebra/TensorAlgebra.jl @@ -6,9 +6,11 @@ using StaticArrays using LinearAlgebra import Base: * import Base: + +import Base: - export (*) export (+) +export (-) export (⊗₁₂³) export (⊗₁₃²) export (⊗₁²³) diff --git a/test/TestTensorAlgebra/TensorAlgebraTests.jl b/test/TestTensorAlgebra/TensorAlgebraTests.jl index 8b34bad..d99e96d 100644 --- a/test/TestTensorAlgebra/TensorAlgebraTests.jl +++ b/test/TestTensorAlgebra/TensorAlgebraTests.jl @@ -107,6 +107,15 @@ end B = TensorValue(4.1, 5.2, 6.3, 7.4, 8.5, 9.6, 1.7, 2.8, 3.9) @test A + B == TensorValue(5.1, 7.2, 9.3, 11.4, 13.5, 15.6, 8.7, 10.8, 12.9) @test norm(A + B) ≈ 32.842807431765024 + @test B - A ≈ TensorValue(3.1, 3.2, 3.3, 3.4, 3.5, 3.6, -5.3, -5.2, -5.1) + @test -A == TensorValue(-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0) + C = TensorValue{2,3}(1.0, 2.0, 3.0, 4.0, 5.0, 6.0) + @test C + C == TensorValue{2,3}(2.0, 4.0, 6.0, 8.0, 10.0, 12.0) + @test TensorValue{2,2}(1, 2, 3, 4) + TensorValue{2,2}(0.5, 0.5, 0.5, 0.5) isa TensorValue{2,2,Float64} + E = TensorValue{0,3,Float64}(()) # empty tensor, built by Gridap for vertex grids + @test E + E == E + @test E - E == E + @test -E == E end