Optimization bnd pressure computation

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -324,13 +324,32 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation,
         v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)
         u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
 
-        nhs = get_neighborhood_search(system, neighbor_system, semi)
-
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
-        adami_pressure_extrapolation!(boundary_model, system, neighbor_system,
-                                      system_coords, neighbor_coords,
-                                      v_neighbor_system, nhs)
+        # This is an optimization for simulations with large and complex boundaries.
+        # Especially, in 3D simulations with large and/or complex structures outside
+        # of areas with permanent flow.
+        # Note: The version iterating neighbors first is not thread parallizable.
+        # The factor is based on the achievable speed-up of the thread parallizable version.
+        if nparticles(system) >
+           ceil(Int, 0.5 * Threads.nthreads()) * nparticles(neighbor_system)
+            nhs = get_neighborhood_search(neighbor_system, system, semi)
+
+            adami_pressure_extrapolation_neighbor!(boundary_model, system, neighbor_system,
+                                                   system_coords, neighbor_coords,
+                                                   v_neighbor_system, nhs)
+        else
+            nhs = get_neighborhood_search(system, neighbor_system, semi)
+
+            adami_pressure_extrapolation!(boundary_model, system, neighbor_system,
+                                          system_coords, neighbor_coords,
+                                          v_neighbor_system, nhs)
+        end
+        @simd for particle in eachparticle(system)
+            # Limit pressure to be non-negative to avoid attractive forces between fluid and
+            # boundary particles at free surfaces (sticking artifacts).
+            pressure[particle] = max(pressure[particle], 0.0)
+        end
     end
 
     @trixi_timeit timer() "inverse state equation" @threaded for particle in eachparticle(system)
@@ -365,6 +384,35 @@ function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZe
     return boundary_model
 end
 
+@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system,
+                                                        neighbor_system::FluidSystem,
+                                                        system_coords, neighbor_coords,
+                                                        v_neighbor_system,
+                                                        neighborhood_search)
+    (; pressure, cache, viscosity) = boundary_model
+
+    for_particle_neighbor(neighbor_system, system,
+                          neighbor_coords, system_coords,
+                          neighborhood_search;
+                          particles=eachparticle(neighbor_system),
+                          parallel=false) do neighbor, particle,
+                                             pos_diff, distance
+        # since neighbor and particle are switched
+        pos_diff = -pos_diff
+        adami_pressure_inner!(boundary_model, system, neighbor_system::FluidSystem,
+                              v_neighbor_system, particle, neighbor, pos_diff,
+                              distance, viscosity, cache, pressure)
+    end
+end
+
+@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system,
+                                                        neighbor_system,
+                                                        system_coords, neighbor_coords,
+                                                        v_neighbor_system,
+                                                        neighborhood_search)
+    return boundary_model
+end
+
 @inline function adami_pressure_extrapolation!(boundary_model, system,
                                                neighbor_system::FluidSystem,
                                                system_coords, neighbor_coords,
@@ -377,28 +425,9 @@ end
                           neighborhood_search;
                           particles=eachparticle(system)) do particle, neighbor,
                                                              pos_diff, distance
-        density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor)
-
-        resulting_acc = neighbor_system.acceleration -
-                        current_acceleration(system, particle)
-
-        kernel_weight = smoothing_kernel(boundary_model, distance)
-
-        pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system,
-                                                 neighbor) +
-                               dot(resulting_acc, density_neighbor * pos_diff)) *
-                              kernel_weight
-
-        cache.volume[particle] += kernel_weight
-
-        compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
-                                   kernel_weight, particle, neighbor)
-    end
-
-    for particle in eachparticle(system)
-        # Limit pressure to be non-negative to avoid attractive forces between fluid and
-        # boundary particles at free surfaces (sticking artifacts).
-        pressure[particle] = max(pressure[particle], 0.0)
+        adami_pressure_inner!(boundary_model, system, neighbor_system::FluidSystem,
+                              v_neighbor_system, particle, neighbor, pos_diff,
+                              distance, viscosity, cache, pressure)
     end
 end
 
@@ -408,6 +437,28 @@ end
     return boundary_model
 end
 
+@inline function adami_pressure_inner!(boundary_model, system,
+                                       neighbor_system::FluidSystem,
+                                       v_neighbor_system, particle, neighbor, pos_diff,
+                                       distance, viscosity, cache, pressure)
+    density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor)
+
+    resulting_acc = neighbor_system.acceleration -
+                    current_acceleration(system, particle)
+
+    kernel_weight = smoothing_kernel(boundary_model, distance)
+
+    pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system,
+                                             neighbor) +
+                           dot(resulting_acc, density_neighbor * pos_diff)) *
+                          kernel_weight
+
+    cache.volume[particle] += kernel_weight
+
+    compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
+                               kernel_weight, particle, neighbor)
+end
+
 function compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
                                     kernel_weight, particle, neighbor)
     return cache

test/validation/validation.jl

@@ -46,14 +46,16 @@
         ]
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0
-        @test isapprox(error_edac_P1, 0, atol=1e-9)
-        @test isapprox(error_edac_P2, 0, atol=1e-11)
 
         if VERSION >= v"1.10"
+            @test isapprox(error_edac_P1, 0, atol=1e-9)
+            @test isapprox(error_edac_P2, 0, atol=6e-12)
             @test isapprox(error_wcsph_P1, 0, atol=eps())
             @test isapprox(error_wcsph_P2, 0, atol=eps())
         else
             # 1.9 causes a large difference in the solution
+            @test isapprox(error_edac_P1, 0, atol=1e-8)
+            @test isapprox(error_edac_P2, 0, atol=1e-10)
             @test isapprox(error_wcsph_P1, 0, atol=10)
             @test isapprox(error_wcsph_P2, 0, atol=1e-3)
         end