WIP: MatrixFree and Assembled GPU operators

[termi-official]

This PR shows how to implement matrix-free operators on the GPU and provides the necessary infrastructure.

TODOs

• Fix bug when block size > 1
• Use Float32
• Remove problematic string interpolations in exceptions (breaks GPU compilation) - also check out Base.LazyString
• Performance comparison
• GPUCellValues/MFCellValues or similar
• Move adapt stuff into extension
• Documentation entries
• Document example
• Check if AMDGPU.jl works, too - maybe even move to KernelAbstractions.jl

Future work

• Constraint handler integration
• Move from mass to steady state heat problem
• Matrix-Free p-Multigrid example from https://arxiv.org/pdf/2204.01722.pdf with help from https://dl.acm.org/doi/pdf/10.1145/3322813

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 92.77%. Comparing base (1a84b77) to head (08256af).
Report is 11 commits behind head on master.

❗ Current head 08256af differs from pull request most recent head 2b7bb24. Consider uploading reports for the commit 2b7bb24 to get more accurate results

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #766      +/-   ##
==========================================
- Coverage   93.29%   92.77%   -0.53%     
==========================================
  Files          36       33       -3     
  Lines        5235     4952     -283     
==========================================
- Hits         4884     4594     -290     
- Misses        351      358       +7     

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

[KnutAM]

Driveby comment: I think we could change

Ferrite.jl/src/interpolations.jl

Line 521 in 3cff601

throw(ArgumentError("no shape function $i for interpolation $ip"))

to @boundscheck throw(ArgumentError("no shape function $i for interpolation $ip"))

(also in all other places) and then in this example use @inbounds?

[termi-official]

Well, this is rather hacky. Give me a bit more time to think about this.

[KnutAM]

Right, this is actually not type stable.
Probably should return eltype(ξ)(NaN) after the @boundscheck and update the operations to not multiply by floats.

But I wouldn't call that hacky though since it just allows @inbounds to be used to elude the bounds check?

[termi-official]

The first SpMV kernels run, but there is still quite a bit of features and optimizations to implement. Here some first benchmarks. cuSPARSE currently wins by a landslide and the full matrix-free kernel is about 3x slower thatn the CPU implementation for now.

julia> @benchmark CUDA.@sync mul!($u,$Agpu,$bgpu)
BenchmarkTools.Trial: 27 samples with 1 evaluation.
 Range (min … max):  185.123 ms … 186.871 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     186.140 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   186.102 ms ± 469.661 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

                      ▃ ▃            █          ▃      █         
  ▇▁▁▁▁▁▁▁▇▁▁▁▁▁▇▇▁▁▁▇█▁█▇▁▁▁▇▇▁▁▁▁▁▇█▁▁▁▁▇▇▁▁▁▁█▁▁▇▁▇▁█▁▁▁▁▇▁▇ ▁
  185 ms           Histogram: frequency by time          187 ms <

 Memory estimate: 46.61 KiB, allocs estimate: 475.

julia> @benchmark mul!($ucpu, $Kcpu, $bcpu)
BenchmarkTools.Trial: 63 samples with 1 evaluation.
 Range (min … max):  79.763 ms …  83.167 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     80.435 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   80.446 ms ± 512.121 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

              █                                                 
  ▆▁▅▆▁▆▅▆▅▁▆██▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅ ▁
  79.8 ms       Histogram: log(frequency) by time        83 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark mul!($ugpu,$Kgpu,$bgpu)
BenchmarkTools.Trial: 8654 samples with 1 evaluation.
 Range (min … max):   18.390 μs …   2.875 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     610.724 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   574.855 μs ± 144.657 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

                                                              █  
  ▃▂▁▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂█ ▂
  18.4 μs          Histogram: frequency by time          614 μs <

 Memory estimate: 1.44 KiB, allocs estimate: 60.

julia> @benchmark mul!($ugpu,$Aeagpu,$bgpu)
BenchmarkTools.Trial: 201 samples with 1 evaluation.
 Range (min … max):  24.722 ms …  30.899 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     24.823 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   24.871 ms ± 437.805 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▃▃█▅▄▁                                                  
  ▂▂▃▃▄███████▆▅▄▄▄▃▁▁▁▁▂▂▂▂▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▃▂ ▃
  24.7 ms         Histogram: frequency by time         25.4 ms <

 Memory estimate: 4.42 KiB, allocs estimate: 151.