[Proof of Concept] Performance comparison with DualSPHysics

[efaulhaber]

Based on #1147, #1211 and #1264. Needs to be used together with trixi-framework/PointNeighbors.jl#167.

This PR adds

• a new example file to replicate the 3D dam break setup from DualSPHysics.
• a new time integration scheme SymplecticPositionVerletWithSorting that is the same as SymplecticPositionVerlet, but sorts particles before every single RHS call. This is necessary to use the algorithm of DualSPHysics. Note that I did not implement the fast GPU sorting that they have, so that we can benchmark the interactions kernel but not the total runtime, as the sorting will now use most runtime.
• a new fluid-fluid interact! function that implements the algorithm used by DualSPHysics to benchmark in a hardcoded custom kernel to benchmark the fluid-fluid interaction with the same algorithm.
• another new fluid-fluid interact! function that also implements the same algorithm, but wrapped in our user-friendly API. This version trades a bit of raw performance for better readability.

Run the benchmark with

trixi_include_changeprecision(Float32, "examples/fluid/dam_break_3d_dualsphysics.jl", fluid_particle_spacing=0.002f0, fluid_dt=5f-5, parallelization_backend=CUDABackend(), tspan=(0.0, 0.1f0), saving_callback=nothing);

Here is a timer output on an RTX A4500 (with a different dt and tspan though):

──────────────────────────────────────────────────────────────────────────────────────────────────
               TrixiParticles.jl                         Time                    Allocations
                                                ───────────────────────   ────────────────────────
               Tot / % measured:                      918s /  99.0%            301GiB / 100.0%

Section                                 ncalls     time    %tot     avg     alloc    %tot      avg
──────────────────────────────────────────────────────────────────────────────────────────────────
sort                                     2.00k     686s   75.5%   343ms    300GiB   99.8%   154MiB
  sort systems                           2.00k     486s   53.5%   243ms    150GiB   49.9%  76.7MiB
  sort arrays                            2.00k     159s   17.5%  79.5ms    220MiB    0.1%   112KiB
  permutation                            2.00k    32.4s    3.6%  16.2ms    150GiB   49.8%  76.6MiB
  deactivate out of bounds particle      2.00k    8.38s    0.9%  4.18ms   48.6MiB    0.0%  24.8KiB
  ~sort~                                 2.00k   7.04ms    0.0%  3.51μs    628KiB    0.0%     321B
kick!                                    2.00k     221s   24.3%   110ms    413MiB    0.1%   211KiB
  system interaction                     2.00k     161s   17.7%  80.2ms    140MiB    0.0%  71.8KiB
    fluid1-*                             2.00k     156s   17.2%  78.1ms   43.5MiB    0.0%  22.2KiB
    boundary2-fluid1                     2.00k    4.24s    0.5%  2.12ms   88.0MiB    0.0%  45.0KiB
    ~system interaction~                 2.00k   23.0ms    0.0%  11.5μs   8.93MiB    0.0%  4.56KiB
    boundary2-boundary2                  2.00k    103μs    0.0%  51.6ns     0.00B    0.0%    0.00B
  update systems and nhs                 2.00k    58.7s    6.4%  29.3ms    241MiB    0.1%   123KiB
    update nhs                           2.00k    30.9s    3.4%  15.4ms    143MiB    0.0%  73.1KiB
    deactivate out of bounds particle    2.00k    26.1s    2.9%  13.1ms   48.7MiB    0.0%  24.9KiB
    ~update systems and nhs~             2.00k    1.65s    0.2%   825μs   49.3MiB    0.0%  25.2KiB
    update density diffusion             2.00k    146μs    0.0%  72.8ns     0.00B    0.0%    0.00B
  source terms                           2.00k    1.48s    0.2%   738μs   27.7MiB    0.0%  14.2KiB
  reset ∂v/∂t                            2.00k   23.6ms    0.0%  11.8μs   2.48MiB    0.0%  1.27KiB
  ~kick!~                                2.00k   4.85ms    0.0%  2.42μs   1.84MiB    0.0%     961B
drift!                                   2.00k    1.79s    0.2%   892μs   17.6MiB    0.0%  8.98KiB
sorting callback                             2    490ms    0.1%   245ms    307MiB    0.1%   153MiB
──────────────────────────────────────────────────────────────────────────────────────────────────

and DualSPHysics on the same card:

Particles of simulation (initial): 10,743,136
DTs adjusted to DtMin............: 0
Excluded particles...............: 0
Total Runtime....................: 191.181488 sec.
Simulation Runtime...............: 190.139053 sec.
Runtime per physical second......: 18994.910156 sec.
Steps per second.................: 5.264568
Steps of simulation..............: 1,001
PART files.......................: 1
Maximum number of particles......: 10,743,136
Maximum number of cells..........: 261,900
CPU Memory.......................: 966,939,840 (922.15 MiB)
GPU Memory.......................: 2,017,608,392 (1924.14 MiB)
GPU Memory (for cells)...........: 31,190,584 (29.75 MiB)
 
[GPU Timers]
VA-Init..........................: 1.042425 sec.
NL-Limits........................: 5.859418 sec.
NL-PreSort.......................: 0.003581 sec.
NL-RadixSort.....................: 3.150588 sec.
NL-CellBegin.....................: 0.250372 sec.
NL-SortData......................: 11.852192 sec.
NL-OutCheck......................: 0 sec.
CF-PreForces.....................: 1.866821 sec.
CF-PreMDBC.......................: 0 sec.
CF-Forces........................: 160.079714 sec.
SU-Shifting......................: 0 sec.
SU-ComputeStep...................: 7.074752 sec.
SU-Floating......................: 0 sec.
SU-Motion........................: 0 sec.
SU-Periodic......................: 0 sec.
SU-ResizeNp......................: 0 sec.
SU-DownData......................: 0 sec.
SU-SavePart......................: 0 sec.
SU-Chrono........................: 0 sec.
SU-Moorings......................: 0 sec.
SU-InOut.........................: 0 sec.
SU-Gauges........................: 0 sec.
SU-FlexStruc.....................: 0 sec.

So that is 161s for forces/interactions in TrixiParticles and 160s in DualSPHysics.

Here are some more average interaction runtimes (lower is better!):

Method	NVIDIA RTX A4500	NVIDIA H100
DualSPHysics	Baseline — 81.5 ms	Baseline — 29.2 ms
TrixiParticles, hardcoded (combined) kernel	1.03× — 83.9 ms	0.96× — 28 ms
TrixiParticles, their neighborhood search algorithm	1.29× — 105 ms	1.22× — 35.5 ms
TrixiParticles, standard code, our neighborhood search	1.60× — 130 ms	1.42× — 41.6 ms

And total runtimes:

Method	NVIDIA RTX A4500	NVIDIA H100
DualSPHysics	Baseline — 387 s	Baseline — 132 s
TrixiParticles, standard code, our neighborhood search	1.57× — 609 s	1.36× — 180 s
TrixiParticles, Carpenter-Kennedy, our NHS	1.00× — 387 s	0.86× — 114 s