Add interpolation functions (#304)

* add interpolation

* basic working

* support periodic

* allow to define different smoothing_length

* more precisely identify system regions

* to improve accuracy calc shepard

* add support for multiple points

* add line interpolation

* add doc

* format

* add test

* small change to doc

* revert accidental change

* fix

* another fix

* format

* format

* also interpolate pressure and velocity

* lower tolerance

* add example plot

* missing dependency

* format

* comment out example

* remove pyplot again

* fix

* small update improving accuracy

* fix tests

* remove line

* fix

* format

* fix

* merge

* optimization

* fix result format

* update test

* fix

* fix docs

* rename to copy_neighborhood_search

* review fixes

* format

* bs

* tests

* more tests

* fix all tests

* fix

* fix tests

* increase error since it passes locally but not on github

* fix

* fix

* review comments

* fix docs

* format

* add interpolation.jl to doc gen

* review

* update

* add Plots version

* improve plot

* remove other dispatch

* fix

* fix docs

* merge

* remove background pressure

* fix comment

docs/make.jl

@@ -43,6 +43,7 @@ makedocs(sitename="TrixiParticles.jl",
                      "Semidiscretization" => joinpath("general", "semidiscretization.md"),
                      "Initial Condition and Setups" => joinpath("general",
                                                                 "initial_condition.md"),
+                     "Interpolation" => joinpath("general", "interpolation.md"),
                      "Density Calculators" => joinpath("general", "density_calculators.md"),
                      "Smoothing Kernels" => joinpath("general", "smoothing_kernels.md"),
                      "Neighborhood Search" => joinpath("general", "neighborhood_search.md"),

docs/src/general/interpolation.md

@@ -0,0 +1,6 @@
+# Interpolation
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("general", "interpolation.jl")]
+```

examples/postprocessing/interpolation.jl

@@ -0,0 +1,61 @@
+# Example for using interpolation
+#######################################################################################
+using TrixiParticles
+using Plots
+
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_2d.jl"))
+
+# `interpolate_point` can be used to interpolate the properties of the `fluid_system` with the original kernel and `smoothing_length`
+println(interpolate_point([1.0, 0.01], semi, fluid_system, sol))
+# Or with an increased `smoothing_length` smoothing the result
+println(interpolate_point([1.0, 0.01], semi, fluid_system, sol,
+                          smoothing_length=2.0 * smoothing_length))
+
+# A point outside of the domain will result in properties with value 0
+# On the boundary a result can still be obtained
+println(interpolate_point([1.0, 0.0], semi, fluid_system, sol))
+# Slightly outside of the fluid domain the result is 0
+println(interpolate_point([1.0, -0.01], semi, fluid_system, sol))
+
+# Multiple points can be interpolated by providing an array
+println(interpolate_point([
+                              [1.0, 0.01],
+                              [1.0, 0.1],
+                              [1.0, 0.0],
+                              [1.0, -0.01],
+                              [1.0, -0.05],
+                          ], semi, fluid_system, sol))
+
+# It is also possible to interpolate along a line
+result = interpolate_line([1.0, -0.05], [1.0, 1.0], 10, semi, fluid_system, sol)
+result_endpoint = interpolate_line([1.0, -0.05], [1.0, 1.0], 10, semi, fluid_system, sol,
+                                   endpoint=false)
+
+# Extracting wall distance for the standard and endpoint cases
+walldistance = [coord[2] for coord in result.coord]
+walldistance_endpoint = [coord[2] for coord in result_endpoint.coord]
+
+# Instead of using Plots.jl one can also use PythonPlot which uses matplotlib
+# using PythonPlot
+
+# figure()
+# plot(walldistance, result.density, marker="o", linestyle="-", label="With Endpoint")
+# plot(walldistance_endpoint, result_endpoint.density, marker="x", linestyle="--",
+#      label="Without Endpoint")
+
+# xlabel("Wall distance")
+# ylabel("Density")
+# title("Density Interpolation Along a Line")
+# legend()
+
+# plotshow()
+
+p = plot(walldistance, result.density, marker=:circle, color=:blue, markerstrokecolor=:blue,
+         linewidth=2, label="With Endpoint")
+
+plot!(p, walldistance_endpoint, result_endpoint.density, marker=:xcross, linewidth=2,
+      linestyle=:dash, label="Without Endpoint", color=:orange)
+
+plot!(p, framestyle=:box, legend=:best, xlabel="Wall distance",
+      ylabel="Density", title="Density Interpolation Along a Line", size=(800, 600),
+      dpi=300)

src/TrixiParticles.jl

@@ -58,5 +58,6 @@ export VoxelSphere, RoundSphere, reset_wall!
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection
 export nparticles
+export interpolate_line, interpolate_point
 
 end # module

src/general/corrections.jl

@@ -192,7 +192,7 @@ function compute_shepard_coeff!(system, v, u, v_ode, u_ode, semi,
         system_coords = current_coordinates(u, system)
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
-        neighborhood_search = neighborhood_searches(system, neighbor_system, semi)
+        neighborhood_search = get_neighborhood_search(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff
         for_particle_neighbor(system, neighbor_system, system_coords,
@@ -259,7 +259,7 @@ function compute_correction_values!(system,
         system_coords = current_coordinates(u, system)
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
-        neighborhood_search = neighborhood_searches(system, neighbor_system, semi)
+        neighborhood_search = get_neighborhood_search(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff
         for_particle_neighbor(system, neighbor_system, system_coords,
@@ -402,7 +402,7 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
         v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)
 
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
-        neighborhood_search = neighborhood_searches(system, neighbor_system, semi)
+        neighborhood_search = get_neighborhood_search(system, neighbor_system, semi)
 
         for_particle_neighbor(system, neighbor_system, coordinates, neighbor_coords,
                               neighborhood_search;

src/general/density_calculators.jl

@@ -46,7 +46,7 @@ function summation_density!(system, semi, u, u_ode, density;
         system_coords = current_coordinates(u, system)
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
-        nhs = neighborhood_searches(system, neighbor_system, semi)
+        nhs = get_neighborhood_search(system, neighbor_system, semi)
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff.
         for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, nhs,

src/general/general.jl

@@ -22,3 +22,4 @@ include("corrections.jl")
 include("smoothing_kernels.jl")
 include("initial_condition.jl")
 include("system.jl")
+include("interpolation.jl")

src/general/interpolation.jl

@@ -0,0 +1,237 @@
+@doc raw"""
+    interpolate_line(start, end_, no_points, semi, ref_system, sol; endpoint=true,
+                     smoothing_length=ref_system.smoothing_length, cut_of_bnd=true)
+
+Interpolates properties along a line in an TrixiParticles simulation.
+The line interpolation is accomplished by generating a series of
+evenly spaced points between `start` and `end_`.
+If `endpoint` is `false`, the line is interpolated between the start and end points,
+but does not include these points.
+
+# Arguments
+- `start`: The starting point of the line.
+- `end_`: The ending point of the line.
+- `n_points`: The number of points to interpolate along the line.
+- `semi`: The semidiscretization used for the simulation.
+- `ref_system`: The reference system for the interpolation.
+- `sol`: The solution state from which the properties are interpolated.
+
+# Keywords
+- `cut_off_bnd`: `cut_off_bnd`: Boolean to indicate if quantities should be set to zero when a
+                  point is "closer" to a boundary than to the fluid system
+                  (see an explanation for "closer" below). Defaults to `true`.
+- `endpoint`: A boolean to include (`true`) or exclude (`false`) the end point in the interpolation. Default is `true`.
+- `smoothing_length`: The smoothing length used in the interpolation. Default is `ref_system.smoothing_length`.
+
+# Returns
+- A `NamedTuple` of arrays containing interpolated properties at each point along the line.
+
+!!! note
+    - This function is particularly useful for analyzing gradients or creating visualizations
+      along a specified line in the SPH simulation domain.
+    - The interpolation accuracy is subject to the density of particles and the chosen smoothing length.
+    - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as
+      accurate as a real surface reconstruction.
+
+# Examples
+```julia
+# Interpolating along a line from [1.0, 0.0] to [1.0, 1.0] with 5 points
+results = interpolate_line([1.0, 0.0], [1.0, 1.0], 5, semi, ref_system, sol)
+```
+"""
+function interpolate_line(start, end_, n_points, semi, ref_system, sol; endpoint=true,
+                          smoothing_length=ref_system.smoothing_length,
+                          cut_off_bnd=true)
+    start_svector = SVector{ndims(ref_system)}(start)
+    end_svector = SVector{ndims(ref_system)}(end_)
+    points_coords = range(start_svector, end_svector, length=n_points)
+
+    if !endpoint
+        points_coords = points_coords[2:(end - 1)]
+    end
+
+    return interpolate_point(points_coords, semi, ref_system, sol,
+                             smoothing_length=smoothing_length,
+                             cut_off_bnd=cut_off_bnd)
+end
+
+@doc raw"""
+    interpolate_point(points_coords::Array{Array{Float64,1},1}, semi, ref_system, sol;
+                      smoothing_length=ref_system.smoothing_length, cut_of_bnd=true)
+
+    interpolate_point(point_coords, semi, ref_system, sol;
+                      smoothing_length=ref_system.smoothing_length, cut_of_bnd=true)
+
+Performs interpolation of properties at specified points or an array of points in a TrixiParticles simulation.
+
+When given an array of points (`points_coords`), it iterates over each point and applies interpolation individually.
+For a single point (`point_coords`), it performs the interpolation at that specific location.
+The interpolation utilizes the same kernel function of the SPH simulation to weigh contributions from nearby particles.
+
+# Arguments
+- `points_coords`: An array of point coordinates, for which to interpolate properties.
+- `point_coords`: The coordinates of a single point for interpolation.
+- `semi`: The semidiscretization used in the SPH simulation.
+- `ref_system`: The reference system defining the properties of the SPH particles.
+- `sol`: The current solution state from which properties are interpolated.
+
+# Keywords
+- `cut_off_bnd`: Cut-off at the boundary.
+- `smoothing_length`: The smoothing length used in the kernel function. Defaults to `ref_system.smoothing_length`.
+
+# Returns:
+- For multiple points:  A `NamedTuple` of arrays containing interpolated properties at each point.
+- For a single point: A `NamedTuple` of interpolated properties at the point.
+
+# Examples:
+```julia
+# For a single point
+result = interpolate_point([1.0, 0.5], semi, ref_system, sol)
+
+# For multiple points
+points = [[1.0, 0.5], [1.0, 0.6], [1.0, 0.7]]
+results = interpolate_point(points, semi, ref_system, sol)
+```
+!!! note
+    - This function is particularly useful for analyzing gradients or creating visualizations
+      along a specified line in the SPH simulation domain.
+    - The interpolation accuracy is subject to the density of particles and the chosen smoothing length.
+    - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as
+    accurate as a real surface reconstruction.
+"""
+function interpolate_point(points_coords::AbstractArray{<:AbstractArray}, semi, ref_system,
+                           sol; smoothing_length=ref_system.smoothing_length,
+                           cut_off_bnd=true)
+    num_points = length(points_coords)
+    coords = similar(points_coords)
+    velocities = similar(points_coords)
+    densities = Vector{Float64}(undef, num_points)
+    pressures = Vector{Float64}(undef, num_points)
+    neighbor_counts = Vector{Int}(undef, num_points)
+
+    neighborhood_searches = process_neighborhood_searches(semi, sol, ref_system,
+                                                          smoothing_length)
+
+    for (i, point) in enumerate(points_coords)
+        result = interpolate_point(SVector{ndims(ref_system)}(point), semi, ref_system, sol,
+                                   neighborhood_searches, smoothing_length=smoothing_length,
+                                   cut_off_bnd=cut_off_bnd)
+        densities[i] = result.density
+        neighbor_counts[i] = result.neighbor_count
+        coords[i] = result.coord
+        velocities[i] = result.velocity
+        pressures[i] = result.pressure
+    end
+
+    return (density=densities, neighbor_count=neighbor_counts, coord=coords,
+            velocity=velocities, pressure=pressures)
+end
+
+function interpolate_point(point_coords, semi, ref_system, sol;
+                           smoothing_length=ref_system.smoothing_length,
+                           cut_off_bnd=true)
+    neighborhood_searches = process_neighborhood_searches(semi, sol, ref_system,
+                                                          smoothing_length)
+
+    return interpolate_point(SVector{ndims(ref_system)}(point_coords), semi, ref_system,
+                             sol, neighborhood_searches, smoothing_length=smoothing_length,
+                             cut_off_bnd=cut_off_bnd)
+end
+
+function process_neighborhood_searches(semi, sol, ref_system, smoothing_length)
+    if isapprox(smoothing_length, ref_system.smoothing_length)
+        # Update existing NHS
+        update_nhs(sol.u[end].x[2], semi)
+        neighborhood_searches = semi.neighborhood_searches[system_indices(ref_system, semi)]
+    else
+        ref_smoothing_kernel = ref_system.smoothing_kernel
+        search_radius = compact_support(ref_smoothing_kernel, smoothing_length)
+        neighborhood_searches = map(semi.systems) do system
+            u = wrap_u(sol.u[end].x[2], system, semi)
+            system_coords = current_coordinates(u, system)
+            old_nhs = get_neighborhood_search(ref_system, system, semi)
+            nhs = copy_neighborhood_search(old_nhs, search_radius, system_coords)
+            return nhs
+        end
+    end
+
+    return neighborhood_searches
+end
+
+@inline function interpolate_point(point_coords, semi, ref_system, sol,
+                                   neighborhood_searches;
+                                   smoothing_length=ref_system.smoothing_length,
+                                   cut_off_bnd=true)
+    interpolated_density = 0.0
+    interpolated_velocity = zero(SVector{ndims(ref_system)})
+    interpolated_pressure = 0.0
+
+    shepard_coefficient = 0.0
+    ref_id = system_indices(ref_system, semi)
+    neighbor_count = 0
+    other_density = 0.0
+    ref_smoothing_kernel = ref_system.smoothing_kernel
+
+    if cut_off_bnd
+        systems = semi
+    else
+        # Don't loop over other systems
+        systems = (ref_system,)
+    end
+
+    foreach_system(systems) do system
+        system_id = system_indices(system, semi)
+        nhs = neighborhood_searches[system_id]
+        (; search_radius, periodic_box) = nhs
+
+        v = wrap_v(sol.u[end].x[1], system, semi)
+        u = wrap_u(sol.u[end].x[2], system, semi)
+
+        system_coords = current_coordinates(u, system)
+
+        # This is basically `for_particle_neighbor` unrolled
+        for particle in eachneighbor(point_coords, nhs)
+            coords = extract_svector(system_coords, Val(ndims(system)), particle)
+
+            pos_diff = point_coords - coords
+            distance2 = dot(pos_diff, pos_diff)
+            pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2,
+                                                            search_radius, periodic_box)
+            if distance2 > search_radius^2
+                continue
+            end
+
+            distance = sqrt(distance2)
+            mass = hydrodynamic_mass(system, particle)
+            kernel_value = kernel(ref_smoothing_kernel, distance, smoothing_length)
+            m_W = mass * kernel_value
+
+            if system_id == ref_id
+                interpolated_density += m_W
+
+                volume = mass / particle_density(v, system, particle)
+                particle_velocity = current_velocity(v, system, particle)
+                interpolated_velocity += particle_velocity * (volume * kernel_value)
+
+                pressure = particle_pressure(v, system, particle)
+                interpolated_pressure += pressure * (volume * kernel_value)
+                shepard_coefficient += volume * kernel_value
+            else
+                other_density += m_W
+            end
+
+            neighbor_count += 1
+        end
+    end
+
+    # Point is not within the ref_system
+    if other_density > interpolated_density || shepard_coefficient < eps()
+        return (density=0.0, neighbor_count=0, coord=point_coords,
+                velocity=zero(SVector{ndims(ref_system)}), pressure=0.0)
+    end
+
+    return (density=interpolated_density / shepard_coefficient,
+            neighbor_count=neighbor_count,
+            coord=point_coords, velocity=interpolated_velocity / shepard_coefficient,
+            pressure=interpolated_pressure / shepard_coefficient)
+end