Further improve TLSPH deformation gradient performance on GPUs in 2D

NEWS.md

@@ -4,7 +4,7 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
-## Version 0.4.4
+## Version 0.5
 
 ### API Changes
 
@@ -21,6 +21,13 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - Added new validation case hydrostatic water column (#724).
 - Added Carreau–Yasuda non-Newtonian viscosity model (#1010).
 
+### Performance
+
+- Greatly improved GPU performance of WCSPH and TLSPH
+  (#1128, #1117, #1124, #1125, #1130, #1116, #1139, #1149).
+  See [#1131](https://github.com/trixi-framework/TrixiParticles.jl/issues/1131)
+  for a detailed breakdown including benchmark results.
+
 ### Important Bugfixes
 
 - Fixed the periodic array of cylinders example file (#975).

src/schemes/structure/total_lagrangian_sph/rhs.jl

@@ -71,7 +71,8 @@ end
             F_b = @inbounds deformation_gradient(system, neighbor)
 
             current_pos_diff_ = current_coords_a - current_coords_b
-            # On GPUs, convert `Float64` coordinates to `Float32` after computing the difference
+            # In mixed-precision simulations, convert from `coordinates_eltype(system)`
+            # to `eltype(system)` immediately after computing the difference.
             current_pos_diff = convert.(eltype(system), current_pos_diff_)
             current_distance = norm(current_pos_diff)
 

src/schemes/structure/total_lagrangian_sph/system.jl

@@ -486,9 +486,6 @@ end
 @inline function calc_deformation_grad!(deformation_grad, system, semi)
     (; mass, material_density) = system
 
-    # Reset deformation gradient
-    set_zero!(deformation_grad)
-
     # For `distance == 0`, the analytical gradient is zero, but the unsafe gradient
     # and the density diffusion divide by zero.
     # To account for rounding errors, we check if `distance` is almost zero.
@@ -530,10 +527,15 @@ end
             grad_kernel = smoothing_kernel_grad_unsafe(system, initial_pos_diff,
                                                        initial_distance, particle)
 
-            volume = @inbounds mass[neighbor] / material_density[neighbor]
+            # Since this is one of the most performance critical functions, using fast
+            # divisions here gives a significant speedup on GPUs.
+            # See the docs page "Development" for more details on `div_fast`.
+            volume = @inbounds div_fast(mass[neighbor], material_density[neighbor])
             current_coords_b = @inbounds current_coords(system, neighbor)
+
             pos_diff_ = current_coords_a - current_coords_b
-            # On GPUs, convert `Float64` coordinates to `Float32` after computing the difference
+            # In mixed-precision simulations, convert from `coordinates_eltype(system)`
+            # to `eltype(system)` immediately after computing the difference.
             pos_diff = convert.(eltype(system), pos_diff_)
 
             # The tensor product pos_diff ⊗ (L_{0a} * ∇W) is equivalent to multiplication
@@ -542,7 +544,8 @@ end
         end
 
         for j in 1:ndims(system), i in 1:ndims(system)
-            @inbounds deformation_grad[i, j, particle] += result[][i, j]
+            # We overwrite every entry of `deformation_grad`, so no `set_zero!` is required.
+            @inbounds deformation_grad[i, j, particle] = result[][i, j]
         end
     end