Interpolation on GPU

src/general/interpolation.jl

@@ -393,7 +393,13 @@ function interpolate_line(start, end_, n_points, semi, ref_system, v_ode, u_ode;
     # Convert to coordinate matrix
     points_coords_ = collect(reinterpret(reshape, eltype(start_svector), points_coords))
 
-    return interpolate_points(points_coords_, semi, ref_system, v_ode, u_ode;
+    if semi.parallelization_backend isa KernelAbstractions.Backend
+        points_coords_new = Adapt.adapt(semi.parallelization_backend, points_coords_)
+    else
+        points_coords_new = points_coords_
+    end
+
+    return interpolate_points(points_coords_new, semi, ref_system, v_ode, u_ode;
                               smoothing_length=smoothing_length,
                               cut_off_bnd=cut_off_bnd, clip_negative_pressure)
 end
@@ -570,11 +576,14 @@ end
             computed_density[point] = NaN
             neighbor_count[point] = 0
 
-            set_nan!(cache, point, ref_system)
+            foreach(Tuple(cache)) do field
+                set_nan!(field, point)
+            end
         else
             # Normalize all quantities by the shepard coefficient
-            divide_by_shepard_coefficient!(cache, point, shepard_coefficient,
-                                           ref_system)
+            foreach(Tuple(cache)) do field
+                divide_by_shepard_coefficient!(field, shepard_coefficient, point)
+            end
         end
     end
 
@@ -584,55 +593,6 @@ end
             neighbor_count=Array(neighbor_count), cache_...)
 end
 
-function set_nan!(cache, point, ref_system::FluidSystem)
-    for dim in 1:ndims(ref_system)
-        cache.velocity[dim, point] = NaN
-    end
-    cache.pressure[point] = NaN
-    cache.density[point] = NaN
-
-    return cache
-end
-
-function set_nan!(cache, point, ref_system::SolidSystem)
-    for dim in 1:ndims(ref_system)
-        cache.velocity[dim, point] = NaN
-    end
-    cache.jacobian[point] = NaN
-    cache.von_mises_stress[point] = NaN
-
-    for j in axes(cache.cauchy_stress, 2), i in axes(cache.cauchy_stress, 1)
-        cache.cauchy_stress[i, j, point] = NaN
-    end
-
-    return cache
-end
-
-function divide_by_shepard_coefficient!(cache, point, shepard_coefficient,
-                                        ref_system::FluidSystem)
-    for dim in 1:ndims(ref_system)
-        cache.velocity[dim, point] /= shepard_coefficient[point]
-    end
-    cache.pressure[point] /= shepard_coefficient[point]
-    cache.density[point] /= shepard_coefficient[point]
-
-    return cache
-end
-
-function divide_by_shepard_coefficient!(cache, point, shepard_coefficient,
-                                        ref_system::SolidSystem)
-    for dim in 1:ndims(ref_system)
-        cache.velocity[dim, point] /= shepard_coefficient[point]
-    end
-    cache.jacobian[point] /= shepard_coefficient[point]
-    cache.von_mises_stress[point] /= shepard_coefficient[point]
-
-    for j in axes(cache.cauchy_stress, 2), i in axes(cache.cauchy_stress, 1)
-        cache.cauchy_stress[i, j, point] /= shepard_coefficient[point]
-    end
-    return cache
-end
-
 @inline function create_cache_interpolation(ref_system::FluidSystem, n_points, semi)
     (; parallelization_backend) = semi
 
@@ -704,3 +664,49 @@ end
 
     return cache
 end
+
+function set_nan!(field::AbstractVector, point)
+    field[point] = NaN
+
+    return field
+end
+
+function set_nan!(field::AbstractArray{T, 2}, point) where {T}
+    for dim in axes(field, 1)
+        field[dim, point] = NaN
+    end
+
+    return field
+end
+
+function set_nan!(field::AbstractArray{T, 3}, point) where {T}
+    for j in axes(field, 2), i in axes(field, 1)
+        field[i, j, point] = NaN
+    end
+
+    return field
+end
+
+function divide_by_shepard_coefficient!(field::AbstractVector, shepard_coefficient, point)
+    field[point] /= shepard_coefficient[point]
+
+    return field
+end
+
+function divide_by_shepard_coefficient!(field::AbstractArray{T, 2}, shepard_coefficient,
+                                        point) where {T}
+    for dim in axes(field, 1)
+        field[dim, point] /= shepard_coefficient[point]
+    end
+
+    return field
+end
+
+function divide_by_shepard_coefficient!(field::AbstractArray{T, 3}, shepard_coefficient,
+                                        point) where {T}
+    for j in axes(field, 2), i in axes(field, 1)
+        field[i, j, point] /= shepard_coefficient[point]
+    end
+
+    return field
+end

test/examples/gpu.jl

@@ -392,3 +392,72 @@ end
         end
     end
 end
+
+@testset verbose=true "Postprocessing $TRIXIPARTICLES_TEST_" begin
+    @testset verbose=true "Interpolation" begin
+        # Import variables into scope
+        trixi_include_changeprecision(Float32, @__MODULE__,
+                                      joinpath(examples_dir(), "fluid",
+                                               "hydrostatic_water_column_2d.jl"),
+                                      sol=nothing, ode=nothing)
+
+        # Neighborhood search with `FullGridCellList` for GPU compatibility
+        min_corner = minimum(tank.boundary.coordinates, dims=2)
+        max_corner = maximum(tank.boundary.coordinates, dims=2)
+        cell_list = FullGridCellList(; min_corner, max_corner, max_points_per_cell=500)
+        semi_fullgrid = Semidiscretization(fluid_system, boundary_system,
+                                           neighborhood_search=GridNeighborhoodSearch{2}(;
+                                                                                         cell_list),
+                                           parallelization_backend=Main.parallelization_backend)
+
+        trixi_include_changeprecision(Float32, @__MODULE__,
+                                      joinpath(examples_dir(),
+                                               "fluid", "hydrostatic_water_column_2d.jl");
+                                      semi=semi_fullgrid, tspan=(0.0f0, 0.1f0))
+
+        semi_new = TrixiParticles.Adapt.adapt(semi_fullgrid.parallelization_backend,
+                                              sol.prob.p)
+
+        # Interpolation parameters
+        position_x = tank_size[1] / 2
+        n_interpolation_points = 10
+        start_point = [position_x, -fluid_particle_spacing]
+        end_point = [position_x, tank_size[2]]
+
+        result = interpolate_line(start_point, end_point, n_interpolation_points,
+                                  semi_new, semi_new.systems[1], sol; cut_off_bnd=false)
+
+        @test isapprox(result.computed_density[1:(end - 1)], # Exclude last NaN
+                       Float32[62.50176,
+                               1053.805,
+                               1061.2959,
+                               1055.8348,
+                               1043.9069,
+                               1038.2051,
+                               1033.1708,
+                               1014.2249,
+                               672.61566])
+
+        @test isapprox(result.density[1:(end - 1)], # Exclude last NaN
+                       Float32[1078.3738,
+                               1070.8535,
+                               1061.2003,
+                               1052.4126,
+                               1044.5074,
+                               1037.0444,
+                               1028.4813,
+                               1014.7941,
+                               1003.6117])
+
+        @test isapprox(result.pressure[1:(end - 1)], # Exclude last NaN
+                       Float32[9940.595,
+                               8791.842,
+                               7368.837,
+                               6143.6562,
+                               5093.711,
+                               4143.313,
+                               3106.1575,
+                               1552.1078,
+                               366.71414])
+    end
+end