Merge branch 'main' into interpolation-gpu

LasNikas · LasNikas · commit a9b4f32eaaff · 2025-06-02T11:07:25.000+02:00
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.31.2
+        uses: crate-ci/typos@v1.32.0
diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl
@@ -32,7 +32,7 @@ function particle_shifting!(integrator)
     v_ode, u_ode = integrator.u.x
     dt = integrator.dt
     # Internal cache vector, which is safe to use as temporary array
-    u_cache = first(get_tmp_cache(integrator))
+    vu_cache = first(get_tmp_cache(integrator))
 
     # Update quantities that are stored in the systems. These quantities (e.g. pressure)
     # still have the values from the last stage of the previous step if not updated here.
@@ -41,7 +41,7 @@ function particle_shifting!(integrator)
     @trixi_timeit timer() "particle shifting" foreach_system(semi) do system
         u = wrap_u(u_ode, system, semi)
         v = wrap_v(v_ode, system, semi)
-        particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt)
+        particle_shifting!(u, v, system, v_ode, u_ode, semi, vu_cache, dt)
     end
 
     # Tell OrdinaryDiffEq that `u` has been modified
@@ -55,14 +55,18 @@ function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt)
 end
 
 function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi,
-                            u_cache, dt)
+                            vu_cache, dt)
     # Wrap the cache vector to an NDIMS x NPARTICLES matrix.
     # We need this buffer because we cannot safely update `u` while iterating over it.
+    _, u_cache = vu_cache.x
     delta_r = wrap_u(u_cache, system, semi)
     set_zero!(delta_r)
 
-    v_max = maximum(particle -> norm(current_velocity(v, system, particle)),
-                    eachparticle(system))
+    # This has similar performance to `maximum(..., eachparticle(system))`,
+    # but is GPU-compatible.
+    v_max = maximum(x -> sqrt(dot(x, x)),
+                    reinterpret(reshape, SVector{ndims(system), eltype(v)},
+                                current_velocity(v, system)))
 
     # TODO this needs to be adapted to multi-resolution.
     # Section 3.2 explains what else needs to be changed.
@@ -88,7 +92,7 @@ function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi,
             grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
 
             # According to p. 29 below Eq. 9
-            R = 0.2
+            R = 2 // 10
             n = 4
 
             # Eq. 7 in Sun et al. (2017).
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -247,6 +247,23 @@ end
 
 system_correction(system::WeaklyCompressibleSPHSystem) = system.correction
 
+@inline function current_velocity(v, system::WeaklyCompressibleSPHSystem)
+    return current_velocity(v, system.density_calculator, system)
+end
+
+@inline function current_velocity(v, ::SummationDensity,
+                                  system::WeaklyCompressibleSPHSystem)
+    # When using `SummationDensity`, `v` contains only the velocity
+    return v
+end
+
+@inline function current_velocity(v, ::ContinuityDensity,
+                                  system::WeaklyCompressibleSPHSystem)
+    # When using `ContinuityDensity`, the velocity is stored
+    # in the first `ndims(system)` rows of `v`.
+    return view(v, 1:ndims(system), :)
+end
+
 @inline function current_density(v, system::WeaklyCompressibleSPHSystem)
     return current_density(v, system.density_calculator, system)
 end
diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl
@@ -136,14 +136,15 @@ end
                                                parallelization_backend=Main.parallelization_backend)
 
             # Note that this simulation only takes 42 time steps on the CPU.
-            # Due to https://github.com/JuliaGPU/Metal.jl/issues/549, it doesn't work on Metal.
+            # TODO This takes 43 time steps on Metal.
+            # Maybe related to https://github.com/JuliaGPU/Metal.jl/issues/549
             trixi_include_changeprecision(Float32, @__MODULE__,
                                           joinpath(examples_dir(), "fluid",
                                                    "dam_break_3d.jl"),
                                           tspan=(0.0f0, 0.1f0),
                                           fluid_particle_spacing=0.1,
                                           semi=semi_fullgrid,
-                                          maxiters=42)
+                                          maxiters=43)
             @test sol.retcode == ReturnCode.Success
             backend = TrixiParticles.KernelAbstractions.get_backend(sol.u[end].x[1])
             @test backend == Main.parallelization_backend
diff --git a/test/validation/validation.jl b/test/validation/validation.jl
@@ -83,4 +83,33 @@
         # Verify number of plots
         @test length(axs_edac[1].scene.plots) >= 2
     end
+
+    @trixi_testset "TGV_2D" begin
+        @trixi_test_nowarn trixi_include(@__MODULE__,
+                                         joinpath(validation_dir(),
+                                                  "taylor_green_vortex_2d",
+                                                  "validation_taylor_green_vortex_2d.jl"),
+                                         tspan=(0.0, 0.01)) [
+            r"WARNING: Method definition pressure_function.*\n",
+            r"WARNING: Method definition initial_pressure_function.*\n",
+            r"WARNING: Method definition velocity_function.*\n",
+            r"WARNING: Method definition initial_velocity_function.*\n"
+        ]
+        @test sol.retcode == ReturnCode.Success
+        @test count_rhs_allocations(sol, semi) == 0
+    end
+
+    @trixi_testset "LDC_2D" begin
+        @trixi_test_nowarn trixi_include(@__MODULE__,
+                                         joinpath(validation_dir(),
+                                                  "lid_driven_cavity_2d",
+                                                  "validation_lid_driven_cavity_2d.jl"),
+                                         tspan=(0.0, 0.02), dt=0.01,
+                                         SENSOR_CAPTURE_TIME=0.01) [
+            r"WARNING: Method definition lid_movement_function.*\n",
+            r"WARNING: Method definition is_moving.*\n"
+        ]
+        @test sol.retcode == ReturnCode.Success
+        @test count_rhs_allocations(sol, semi) == 0
+    end
 end
diff --git a/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl b/validation/lid_driven_cavity_2d/validation_lid_driven_cavity_2d.jl
@@ -2,17 +2,21 @@ using TrixiParticles
 
 # ==========================================================================================
 # ==== Resolution
-particle_spacings = [0.02, 0.01, 0.005]
+# The paper provides reference data for particle spacings particle_spacings = [0.02, 0.01, 0.005]
+particle_spacing = 0.02
 
 # ==========================================================================================
 # ==== Experiment Setup
 tspan = (0.0, 25.0)
 reynolds_numbers = [100.0, 1000.0, 10_000.0]
 
-interpolated_velocity(v, u, t, system) = nothing
+const SENSOR_CAPTURE_TIME = 24.8
+const CAPTURE_STARTED = Ref(false)
 
-function interpolated_velocity(v, u, t, system::TrixiParticles.FluidSystem)
-    if t < 24.8
+interpolated_velocity(system, v, u, semi, t) = nothing
+
+function interpolated_velocity(system::TrixiParticles.FluidSystem, v, u, semi, t)
+    if t < SENSOR_CAPTURE_TIME
         return nothing
     end
 
@@ -32,7 +36,7 @@ function interpolated_velocity(v, u, t, system::TrixiParticles.FluidSystem)
 
     file = joinpath(output_directory, "interpolated_velocities.csv")
 
-    if isfile(file)
+    if CAPTURE_STARTED[]
         data = TrixiParticles.CSV.read(file, TrixiParticles.DataFrame)
         vy_y_ = (data.vy_y .+ vy_y)
         vy_x_ = (data.vy_x .+ vy_x)
@@ -50,12 +54,13 @@ function interpolated_velocity(v, u, t, system::TrixiParticles.FluidSystem)
                                       counter=1, vy_y=vy_y, vy_x=vy_x, vx_y=vx_y, vx_x=vx_x)
 
         TrixiParticles.CSV.write(output_directory * "/interpolated_velocities.csv", df)
+        CAPTURE_STARTED[] = true
     end
 
     return nothing
 end
 
-for particle_spacing in particle_spacings, reynolds_number in reynolds_numbers,
+for reynolds_number in reynolds_numbers,
     density_calculator in [SummationDensity(), ContinuityDensity()], wcsph in [false, true]
     n_particles_xy = round(Int, 1.0 / particle_spacing)
 
@@ -68,8 +73,9 @@ for particle_spacing in particle_spacings, reynolds_number in reynolds_numbers,
                                        name_density_calculator,
                                        "validation_run_lid_driven_cavity_2d_nparticles_$(n_particles_xy)x$(n_particles_xy)_Re_$Re")
 
-    saving_callback = SolutionSavingCallback(dt=0.1, output_directory=output_directory)
+    saving_callback = SolutionSavingCallback(dt=0.02, output_directory=output_directory)
 
+    CAPTURE_STARTED[] = false
     pp_callback = PostprocessCallback(; dt=0.02,
                                       interpolated_velocity=interpolated_velocity,
                                       filename="interpolated_velocities",
diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl
@@ -2,7 +2,8 @@ using TrixiParticles
 
 # ==========================================================================================
 # ==== Resolution
-particle_spacings = [0.02, 0.01, 0.005]
+# The paper provides reference data for particle spacings particle_spacings = [0.02, 0.01, 0.005]
+particle_spacing = 0.02
 
 # ==========================================================================================
 # ==== Experiment Setup
@@ -12,10 +13,13 @@ reynolds_number = 100.0
 density_calculators = [ContinuityDensity(), SummationDensity()]
 perturb_coordinates = [false, true]
 
-function compute_l1v_error(v, u, t, system)
+function compute_l1v_error(system, v_ode, u_ode, semi, t)
     v_analytical_avg = 0.0
     v_avg = 0.0
 
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+
     for particle in TrixiParticles.eachparticle(system)
         position = TrixiParticles.current_coords(u, system, particle)
 
@@ -32,11 +36,14 @@ function compute_l1v_error(v, u, t, system)
     return v_avg /= v_analytical_avg
 end
 
-function compute_l1p_error(v, u, t, system)
+function compute_l1p_error(system, v_ode, u_ode, semi, t)
     p_max_exact = 0.0
 
     L1p = 0.0
 
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+
     for particle in TrixiParticles.eachparticle(system)
         position = TrixiParticles.current_coords(u, system, particle)
 
@@ -57,14 +64,14 @@ end
 
 # The pressure plotted in the paper is the difference of the local pressure minus
 # the average of the pressure of all particles.
-function diff_p_loc_p_avg(v, u, t, system)
-    p_avg_tot = avg_pressure(v, u, t, system)
+function diff_p_loc_p_avg(system, v, u, semi, t)
+    p_avg_tot = avg_pressure(system, v, u, semi, t)
 
     return v[end, :] .- p_avg_tot
 end
 
 for density_calculator in density_calculators, perturbation in perturb_coordinates,
-    particle_spacing in particle_spacings, wcsph in [false, true]
+    wcsph in [false, true]
     n_particles_xy = round(Int, 1.0 / particle_spacing)
 
     name_density_calculator = density_calculator isa SummationDensity ?
