Skip to content

Convert to voigt and mandel as SVector/SMatrix#194

Closed
KnutAM wants to merge 10 commits intoFerrite-FEM:masterfrom
KnutAM:kam/static_voigt
Closed

Convert to voigt and mandel as SVector/SMatrix#194
KnutAM wants to merge 10 commits intoFerrite-FEM:masterfrom
KnutAM:kam/static_voigt

Conversation

@KnutAM
Copy link
Copy Markdown
Member

@KnutAM KnutAM commented Feb 22, 2023

For local residuals (in e.g. plasticity models) defined as custom types this gives very good speedup when converting to static vectors for iterations:

struct Residual{T}
    a::T
    b::SymmetricTensor{2,3,T,6}
end
function direct(x::Residual{T}) where T
    bs = tosmandel(x.b)
    return vcat(x.a, bs)
end
function via_marray(x::Residual{T}) where T
    y = MVector{7,T}(undef)
    y[1] = x.a
    tomandel!(y, x.b; offset=1)
    return SVector{7,T}(y)
end
function via_array(x::Residual{T}) where T
    y = Vector{T}(undef,7)
    y[1] = x.a
    tomandel!(y, x.b; offset=1)
    return SVector{7,T}(y)
end
x = Residual(rand(), rand(SymmetricTensor{2,3}));
@btime direct($x); # 1.600 ns (0 allocations: 0 bytes)
@btime via_marray($x); # 20.261 ns (1 allocation: 64 bytes)
@btime via_array($x); # 35.650 ns (1 allocation: 112 bytes)
# Edit with mutating version
function unsafe_as_array!(y::Vector, x::Residual)
    @inbounds y[1] = x.a
    @inbounds tomandel!(y, x.b; offset=1)
    return y
end
y = zeros(7);
@btime unsafe_as_array!($y, $x); # 16.032 ns (0 allocations: 0 bytes)

and we can also use this to speed up internal functions, such as the inverse where currently conversions to regular arrays are used (not included now, but can be in this or a later pr):

A=rand(SymmetricTensor{4,3});
@btime inv($A); # 960.000 ns (5 allocations: 3.95 KiB)
sinv(S::TT) where {TT<:FourthOrderTensor} = frommandel(TT, inv(tosmandel(S)));
@btime sinv($A); # 333.628 ns (0 allocations: 0 bytes)
A=rand(SymmetricTensor{4,2}); 
@btime inv($A); # 472.449 ns (5 allocations: 1.98 KiB)
@btime sinv($A); # 12.012 ns (0 allocations: 0 bytes)

Also xref: #182 for speeding up eigenvalue calculation of 4th order (2-dim) tensors.

@KristofferC
Copy link
Copy Markdown
Collaborator

Makes sense. Is there any reason not to just change tovoigt to return this? A user can always just call e.g. collect on it if they want a normal array.

@KnutAM
Copy link
Copy Markdown
Member Author

KnutAM commented Feb 22, 2023

The static version doesn't support custom ordering.

@KnutAM
Copy link
Copy Markdown
Member Author

KnutAM commented Feb 23, 2023

The static version doesn't support custom ordering.

Unless doing something bad of course :) (Or implementing an ordering type)

using Tensors
v64=Tensor{2,2}((11.0, 21, 12, 22));
v32=Tensor{2,2,Float32}(v64);
julia> tosvoigt(v64)
4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 11.0
 22.0
 12.0
 21.0

julia> Tensors.DEFAULT_VOIGT_ORDER[2]
2×2 Matrix{Int64}:
 1  3
 4  2

julia> Tensors.DEFAULT_VOIGT_ORDER[2][1]=2
2

julia> Tensors.DEFAULT_VOIGT_ORDER[2][4]=1
1

julia> Tensors.DEFAULT_VOIGT_ORDER[2]
2×2 Matrix{Int64}:
 2  3
 4  1

julia> tosvoigt(v64)
4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 11.0
 22.0
 12.0
 21.0

julia> tosvoigt(v32)
4-element StaticArraysCore.SVector{4, Float32} with indices SOneTo(4):
 22.0
 11.0
 12.0
 21.0

(Although it should never happen unless users do something very strange, would it make sense to change DEFAULT_VOIGT_ORDER[dim] to an SMatrix?)

@KristofferC
Copy link
Copy Markdown
Collaborator

For local residuals (in e.g. plasticity models) defined as custom types this gives very good speedup when converting to static vectors for iterations:

Question, you only need the voigt format when inserting into the residual and jacobian (which I guess tend to be normal arrays). And you can do that without allocations. So where does the speedup come from?

@KnutAM
Copy link
Copy Markdown
Member Author

KnutAM commented Feb 23, 2023

Using static arrays for the local equation system can be ~2x faster for reasonable equation system sizes, see some benchmarks: https://knutam.github.io/Newton.jl/dev/#Benchmarks (And these are with using RecursiveFactorization for the mutating version, so should be as good as it gets (if not let me know and I can update it:))
Edit 1: And the convenience of not having to pass a cache around which gets slightly cumbersome for ForwardDiff caching
Edit 2: But I haven't tested the pipeline converting a cache'd arrays into static arrays before solving the equation system
Edit 3: Added benchmark of mutating version, faster than using marrays, but slower than the static version.

@codecov-commenter
Copy link
Copy Markdown

codecov-commenter commented Aug 3, 2023

Codecov Report

Patch coverage: 84.25% and project coverage change: -1.45% ⚠️

Comparison is base (800fa66) 98.23% compared to head (ead8221) 96.78%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #194      +/-   ##
==========================================
- Coverage   98.23%   96.78%   -1.45%     
==========================================
  Files          17       17              
  Lines        1244     1338      +94     
==========================================
+ Hits         1222     1295      +73     
- Misses         22       43      +21     
Files Changed Coverage Δ
src/voigt.jl 85.10% <84.25%> (-14.90%) ⬇️

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@kimauth
Copy link
Copy Markdown
Member

kimauth commented Aug 3, 2023

I wanted to have this as well! After a bit of discussion with @KnutAM and @fredrikekre, the API now looks a bit different

julia> A = rand(SymmetricTensor{2,3})                                                                                                                                                                                                                           
3×3 SymmetricTensor{2, 3, Float64, 6}:                                                                                                                                                                                                                          
 0.761033  0.89233   0.280383                                                                                                                                                                                                                                   
 0.89233   0.022293  0.366955                                                                                                                                                                                                                                   
 0.280383  0.366955  0.597853                                                                                                                                                                                                                                   
                                                                                                                                                                                                                                                                
julia> v = tovoigt(A)  # old syntax, PR is not breaking                                                                                                                                                                                                                                         
6-element Vector{Float64}:                                                                                                                                                                                                                                      
 0.761033056378681                                                                                                                                                                                                                                              
 0.02229297696504351                                                                                                                                                                                                                                            
 0.5978529160509549                                                                                                                                                                                                                                             
 0.3669551182206563                                                                                                                                                                                                                                             
 0.28038330882003726                                                                                                                                                                                                                                            
 0.8923299861624863     

julia> v = tovoigt(SVector, A) # tovoigt(SArray, A) also allowed so that we can have dimension agnostic code                                                                                                                                                                                                                               
6-element SVector{6, Float64} with indices SOneTo(6):                                                                                                                                                                                                           
 0.761033056378681                                                                                                                                                                                                                                              
 0.02229297696504351                                                                                                                                                                                                                                            
 0.5978529160509549                                                                                                                                                                                                                                             
 0.3669551182206563                                                                                                                                                                                                                                             
 0.28038330882003726                                                                                                                                                                                                                                            
 0.8923299861624863 

julia> v = tovoigt(SMatrix, A) # not allowed, 2nd order Tensor gives a vector after all
ERROR: MethodError: no method matching tovoigt(::Type{SMatrix}, ::SymmetricTensor{2, 3, Float64, 6})

Closest candidates are:
  tovoigt(::Type{<:SMatrix}, ::Tensor{4, dim, T, M}; order) where {dim, T, M}
   @ Tensors ~/.julia/dev/Tensors/src/voigt.jl:308
  tovoigt(::Type{<:SVector}, ::SymmetricTensor{2, dim, T, M}; offdiagscale, order) where {dim, T, M}
   @ Tensors ~/.julia/dev/Tensors/src/voigt.jl:311
  tovoigt(::Type{<:SMatrix}, ::SymmetricTensor{4, dim, T}; offdiagscale, order) where {dim, T}
   @ Tensors ~/.julia/dev/Tensors/src/voigt.jl:314
  ...

Stacktrace:
 [1] top-level scope
   @ REPL[117]:1

Custom ordering for StaticArrays outputs is now allowed as well:

julia> tovoigt(SVector, A; order = Tensors.DEFAULT_VOIGT_ORDER[3])
6-element SVector{6, Float64} with indices SOneTo(6):
 0.761033056378681
 0.02229297696504351
 0.5978529160509549
 0.3669551182206563
 0.28038330882003726
 0.8923299861624863

@kimauth
Copy link
Copy Markdown
Member

kimauth commented Aug 3, 2023

It turned out that @KnutAM generated functions are blazingly fast, so I figured we should reuse that for tovoigt! too. As this is only possible for the default voigt order, custom ordering falls back to the old implementations (this is true when returning SArrays too).

Here are a couple of benchmarking results. Converting to SArrays with the default voigt order is pretty fast now 🙂

 Row │ tensor type                        number of entries  tovoigt!(v, A)  tovoigt!(v, A; order=custom)  tovoigt(SArray, A)  tovoigt(SArray, A; order=custom) 
     │ DataType                           Int64              Float64         Float64                       Float64             Float64                          
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Tensor{2, 1, Float64, 1}                           1          1.458                         1.458               1.5                              1.458
   2 │ Tensor{4, 1, Float64, 1}                           1          1.5                           1.459               1.458                            1.458
   3 │ SymmetricTensor{2, 1, Float64, 1}                  1          1.459                         1.5                 1.5                              1.459
   4 │ SymmetricTensor{4, 1, Float64, 1}                  1          1.458                         1.5                 1.5                              1.458
   5 │ SymmetricTensor{2, 2, Float64, 3}                  3          1.458                         2.083               1.459                            3.083
   6 │ Tensor{2, 2, Float64, 4}                           4          1.458                         1.75                1.5                              2.875
   7 │ SymmetricTensor{2, 3, Float64, 6}                  6          1.5                           3.458               1.458                            5.084
   8 │ SymmetricTensor{4, 2, Float64, 9}                  9          1.791                        14.4449              1.5                             16.992
   9 │ Tensor{2, 3, Float64, 9}                           9          1.459                         2.458               1.458                            5.5
  10 │ Tensor{4, 2, Float64, 16}                         16          7.666                         5.041               2.125                            8.05005
  11 │ SymmetricTensor{4, 3, Float64, 3…                 36         20.3109                       67.6817             10.8859                          83.5473
  12 │ Tensor{4, 3, Float64, 81}                         81         60.8452                       26.9408             25.7279                          38.2228
Benchmarking code
using Tensors, BenchmarkTools, StaticArrays, DataFrames

to_static_voigt(A; order=nothing) = tovoigt(SArray, A; order)

number_entries = Int[]
tensors = Any[]
mutate_default_order = Float64[]
mutate_custom_order = Float64[]
static_default_order = Float64[]
static_custom_order = Float64[]
for dim in 1:3, TT in (Tensor, SymmetricTensor), order in (2,4)
    println(TT, "{", order, ", ", dim, "}")
    A = rand(TT{order,dim})
    v = tovoigt(A)
    
    push!(tensors, A)
    push!(number_entries, reduce(*, size(v)))

    push!(mutate_default_order, @belapsed @inbounds tovoigt!($v, $A))
    push!(mutate_custom_order, @belapsed @inbounds tovoigt!($v, $A; order=$(Tensors.DEFAULT_VOIGT_ORDER[dim])))
    
    push!(static_default_order, @belapsed to_static_voigt($A))
    push!(static_custom_order, @belapsed to_static_voigt($A; order=$(Tensors.DEFAULT_VOIGT_ORDER[dim])))
end

df = DataFrame("tensor type"=>typeof.(tensors),
              "number of entries" => number_entries,
              "tovoigt!(v, A)" => mutate_default_order .* 1e9,
              "tovoigt!(v, A; order=custom)" => mutate_custom_order .* 1e9,
             "tovoigt(SArray, A)" => static_default_order .* 1e9,
            "tovoigt(SArray, A; order=custom)" => static_custom_order .*1e9)
sort!(df, :("number of entries"))

Comment thread src/voigt.jl Outdated
Copy link
Copy Markdown
Member Author

@KnutAM KnutAM left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for pushing this over the line @kimauth !
I will see if I have time today or tomorrow to try my suggestions to minimize a bit of code redundancy (or someone else feel free to do that)

Comment thread src/voigt.jl
Comment thread src/voigt.jl
@KnutAM
Copy link
Copy Markdown
Member Author

KnutAM commented Aug 3, 2023

Could squeeze a bit more performance for SymmetricTensor{4, 3} and Tensor{4,3}.

PR:

 Row │ tensor type                        number of entries  tovoigt!(v, A)  tovoigt!(v, A; order=custom)  tovoigt(SArray, A)  tovoigt(SArray, A; order=custom) 
     │ DataType                           Int64              Float64         Float64                       Float64             Float64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Tensor{2, 1, Float64, 1}                           1          1.8                           1.4                 1.4                              1.6
   2 │ Tensor{4, 1, Float64, 1}                           1          1.6                           1.6                 1.4                              1.7
   3 │ SymmetricTensor{2, 1, Float64, 1}                  1          1.4                           1.4                 1.6                              1.6
   4 │ SymmetricTensor{4, 1, Float64, 1}                  1          1.5                           1.9                 1.8                              2.0
   5 │ SymmetricTensor{2, 2, Float64, 3}                  3          1.7                           2.9                 2.0                              7.0
   6 │ Tensor{2, 2, Float64, 4}                           4          1.8                           2.4                 2.1                              7.90791
   7 │ SymmetricTensor{2, 3, Float64, 6}                  6          1.7                           4.9                 2.0                              8.60861
   8 │ SymmetricTensor{4, 2, Float64, 9}                  9          2.0                          19.2578              1.9                             16.9339
   9 │ Tensor{2, 3, Float64, 9}                           9          2.0                           3.4                 2.2                              7.3
  10 │ Tensor{4, 2, Float64, 16}                         16          9.1                           7.8                 1.9                              8.20821
  11 │ SymmetricTensor{4, 3, Float64, 3…                 36         23.2932                       97.7754              4.1                             95.7939
  12 │ Tensor{4, 3, Float64, 81}                         81         35.241                        34.5382             21.0421                          36.7573

master

 Row │ tensor type                        number of entries  tovoigt!(v, A)  tovoigt!(v, A; order=custom)  tovoigt(SArray, A)  tovoigt(SArray, A; order=custom) 
    │ DataType                           Int64              Float64         Float64                       Float64             Float64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  1 │ Tensor{2, 1, Float64, 1}                           1         1.4                            1.4                    NaN                               NaN
  2 │ Tensor{4, 1, Float64, 1}                           1         2.2                            2.2                    NaN                               NaN
  3 │ SymmetricTensor{2, 1, Float64, 1}                  1         1.9                            2.3                    NaN                               NaN
  4 │ SymmetricTensor{4, 1, Float64, 1}                  1         2.5                            2.3                    NaN                               NaN
  5 │ SymmetricTensor{2, 2, Float64, 3}                  3         3.3                            3.1                    NaN                               NaN
  6 │ Tensor{2, 2, Float64, 4}                           4         2.8                            2.5                    NaN                               NaN
  7 │ SymmetricTensor{2, 3, Float64, 6}                  6         5.0                            5.0                    NaN                               NaN
  8 │ SymmetricTensor{4, 2, Float64, 9}                  9        21.7653                        20.3611                 NaN                               NaN
  9 │ Tensor{2, 3, Float64, 9}                           9         4.7                            3.9                    NaN                               NaN
 10 │ Tensor{4, 2, Float64, 16}                         16         8.90891                        8.2                    NaN                               NaN
 11 │ SymmetricTensor{4, 3, Float64, 3…                 36       101.801                        108.091                  NaN                               NaN
 12 │ Tensor{4, 3, Float64, 81}                         81        39.0515                        35.9879                 NaN                               NaN

@termi-official
Copy link
Copy Markdown
Member

Oopsie?

@KnutAM
Copy link
Copy Markdown
Member Author

KnutAM commented Aug 4, 2023

Oopsie?

Merged directly via commit

@KnutAM KnutAM deleted the kam/static_voigt branch August 4, 2023 10:54
@fredrikekre
Copy link
Copy Markdown
Member

Yea, it was easier like that. It is a shame Github doesn't understand the link between the commit and the branch/PR.

@KnutAM KnutAM mentioned this pull request Aug 4, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants