Plane interpolation (#306)

* 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 * add func * add doc * add test * add another test * comment out example * remove pyplot again * fix * format * small update improve accuracy * update * up * fix * format * fix plots * fix plots.jl plot * cleanup * fix tests * cleanup * cleanup * fix docs * fix doc * remove unused func * remove unused code * fix * fix viscosity * working 3d example for plane interpolation * add 2d plot as well for 3d case * add docs * fix tests * fix * fix plot * fix comments * fix comment * update * format * add Plots * fix * adjust plot * reduce time * Update interpolation.jl * Update interpolation_plane.jl
trixi-framework · Jan 25, 2024 · 898d8d8 · 898d8d8
1 parent bcc33d7
commit 898d8d8
Show file tree

Hide file tree

Showing 9 changed files with 649 additions and 3 deletions.
diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl
@@ -36,6 +36,7 @@ smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
+viscosity_wall = ViscosityAdami(nu=1.0)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
@@ -44,11 +45,13 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
 
 # ==========================================================================================
 # ==== Boundary
+
 boundary_density_calculator = AdamiPressureExtrapolation()
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
                                              boundary_density_calculator,
-                                             smoothing_kernel, smoothing_length)
+                                             smoothing_kernel, smoothing_length,
+                                             viscosity=viscosity_wall)
 
 # Uncomment to use repulsive boundary particles by Monaghan & Kajtar.
 # Also change spacing ratio and boundary layers (see comment above).

diff --git a/examples/fluid/rectangular_tank_3d.jl b/examples/fluid/rectangular_tank_3d.jl
@@ -0,0 +1,8 @@
+using TrixiParticles
+
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_2d.jl"),
+              fluid_particle_spacing=0.05, initial_fluid_size=(1.0, 1.0, 0.9),
+              tank_size=(1.0, 1.0, 1.2), acceleration=(0.0, 0.0, -9.81),
+              smoothing_kernel=SchoenbergCubicSplineKernel{3}(), tspan=(0.0, 1.0),
+              maxiters=10^5, fluid_density_calculator=ContinuityDensity(),
+              clip_negative_pressure=false)
diff --git a/examples/postprocessing/interpolation_plane.jl b/examples/postprocessing/interpolation_plane.jl
@@ -0,0 +1,91 @@
+# Example for using interpolation
+#######################################################################################
+using TrixiParticles
+using Plots
+using Plots.PlotMeasures
+
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_2d.jl"),
+              tspan=(0.0, 0.1))
+
+# Interpolation parameters
+interpolation_start = [0.0, 0.0]
+interpolation_end = [1.0, 1.0]
+resolution = 0.005
+
+# We can interpolate a plane by providing the lower left and top right coordinates of a plane.
+# Per default the same `smoothing_length` will be used as during the simulation.
+original_plane = interpolate_plane_2d(interpolation_start, interpolation_end, resolution,
+                                      semi, fluid_system, sol)
+original_x = [point[1] for point in original_plane.coord]
+original_y = [point[2] for point in original_plane.coord]
+original_density = original_plane.density
+
+# Plane with double smoothing length
+# Utilizing a higher `smoothing_length` in SPH interpolation increases the amount of smoothing,
+# thereby reducing the visibility of disturbances. It also increases the distance
+# from free surfaces where the fluid is cut off. This adjustment in `smoothing_length`
+# can affect both the accuracy and smoothness of the interpolated results.
+double_smoothing_plane = interpolate_plane_2d(interpolation_start, interpolation_end,
+                                              resolution, semi, fluid_system, sol,
+                                              smoothing_length=2.0 * smoothing_length)
+double_x = [point[1] for point in double_smoothing_plane.coord]
+double_y = [point[2] for point in double_smoothing_plane.coord]
+double_density = double_smoothing_plane.density
+
+# Plane with half smoothing length
+# Employing a lower `smoothing_length` in SPH interpolation reduces the amount of smoothing,
+# consequently increasing the visibility of disturbances. Simultaneously, it allows for a more
+# precise cutoff of the fluid at free surfaces. This change in `smoothing_length` can impact the
+# balance between the detail of disturbances captured and the precision of fluid representation near surfaces.
+half_smoothing_plane = interpolate_plane_2d(interpolation_start, interpolation_end,
+                                            resolution,
+                                            semi, fluid_system, sol,
+                                            smoothing_length=0.5 * smoothing_length)
+half_x = [point[1] for point in half_smoothing_plane.coord]
+half_y = [point[2] for point in half_smoothing_plane.coord]
+half_density = half_smoothing_plane.density
+
+scatter1 = scatter(original_x, original_y, zcolor=original_density, marker=:circle,
+                   markersize=2, markercolor=:viridis, markerstrokewidth=0)
+scatter2 = scatter(double_x, double_y, zcolor=double_density, marker=:circle, markersize=2,
+                   markercolor=:viridis, markerstrokewidth=0)
+scatter3 = scatter(half_x, half_y, zcolor=half_density, marker=:circle, markersize=2,
+                   markercolor=:viridis, markerstrokewidth=0)
+
+plot1 = plot(scatter1, xlabel="X Coordinate", ylabel="Y Coordinate",
+             title="Density Distribution", colorbar_title="Density", ylim=(0.0, 1.0),
+             legend=false, clim=(1000, 1010), colorbar=true)
+plot2 = plot(scatter2, xlabel="X Coordinate", ylabel="Y Coordinate",
+             title="Density with 2x Smoothing Length", colorbar_title="Density",
+             ylim=(0.0, 1.0), legend=false, clim=(1000, 1010), colorbar=true)
+plot3 = plot(scatter3, xlabel="X Coordinate", ylabel="Y Coordinate",
+             title="Density with 0.5x Smoothing Length", colorbar_title="Density",
+             ylim=(0.0, 1.0), legend=false, clim=(1000, 1010), colorbar=true)
+
+trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_3d.jl"),
+              tspan=(0.0, 0.01), initial_fluid_size=(2.0, 1.0, 0.9),
+              tank_size=(2.0, 1.0, 1.2))
+
+# Interpolation parameters
+p1 = [0.0, 0.0, 0.0]
+p2 = [1.0, 1.0, 0.1]
+p3 = [1.0, 0.5, 0.2]
+resolution = 0.025
+
+# We can also interpolate a 3D plane but in this case we have to provide 3 points instead!
+original_plane = interpolate_plane_3d(p1, p2, p3, resolution, semi,
+                                      fluid_system, sol)
+original_x = [point[1] for point in original_plane.coord]
+original_y = [point[2] for point in original_plane.coord]
+original_z = [point[3] for point in original_plane.coord]
+original_density = original_plane.density
+
+scatter_3d = scatter3d(original_x, original_y, original_z, marker_z=original_density,
+                       color=:viridis, markerstrokewidth=0)
+
+plot_3d = plot(scatter_3d, xlabel="X", ylabel="Y", zlabel="Z",
+               title="3D Scatter Plot with Density Coloring", legend=false,
+               clim=(1000, 1010), colorbar=false)
+
+combined_plot = plot(plot1, plot2, plot3, plot_3d, layout=(2, 2),
+                     size=(1000, 1500), margin=3mm)
diff --git a/examples/postprocessing/interpolation.jl → ...ostprocessing/interpolation_point_line.jl b/examples/postprocessing/interpolation.jl → ...ostprocessing/interpolation_point_line.jl
@@ -1,6 +1,7 @@
 # Example for using interpolation
 #######################################################################################
 using TrixiParticles
+# this needs to be commented out to use PythonPlot
 using Plots
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "rectangular_tank_2d.jl"))

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
@@ -58,6 +58,6 @@ export VoxelSphere, RoundSphere, reset_wall!
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection
 export nparticles
-export interpolate_line, interpolate_point
+export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d
 
 end # module
diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl
@@ -1,8 +1,175 @@
+using LinearAlgebra
+
+@doc raw"""
+    interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, sol;
+                         smoothing_length=ref_system.smoothing_length, cut_off_bnd=true)
+
+Interpolates properties along a plane in a TrixiParticles simulation.
+The region for interpolation is defined by its lower left and top right corners,
+with a specified resolution determining the density of the interpolation points.
+
+The function generates a grid of points within the defined region,
+spaced uniformly according to the given resolution.
+
+# Arguments
+- `min_corner`: The lower left corner of the interpolation region.
+- `max_corner`: The top right corner of the interpolation region.
+- `resolution`: The distance between adjacent interpolation points in the grid.
+- `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`.
+- `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 within the plane.
+
+!!! note
+    - 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 across a plane from [0.0, 0.0] to [1.0, 1.0] with a resolution of 0.2
+results = interpolate_plane_2d([0.0, 0.0], [1.0, 1.0], 0.2, semi, ref_system, sol)
+```
+"""
+function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system, sol;
+                              smoothing_length=ref_system.smoothing_length,
+                              cut_off_bnd=true)
+    dims = length(min_corner)
+    if dims != 2 || length(max_corner) != 2
+        throw(ArgumentError("function is intended for 2D coordinates only"))
+    end
+
+    if any(min_corner .> max_corner)
+        throw(ArgumentError("`min_corner` should be smaller than `max_corner` in every dimension"))
+    end
+
+    # Calculate the number of points in each dimension based on the resolution
+    no_points_x = ceil(Int, (max_corner[1] - min_corner[1]) / resolution) + 1
+    no_points_y = ceil(Int, (max_corner[2] - min_corner[2]) / resolution) + 1
+
+    x_range = range(min_corner[1], max_corner[1], length=no_points_x)
+    y_range = range(min_corner[2], max_corner[2], length=no_points_y)
+
+    # Generate points within the plane
+    points_coords = [SVector(x, y) for x in x_range, y in y_range]
+
+    results = interpolate_point(points_coords, semi, ref_system, sol,
+                                smoothing_length=smoothing_length,
+                                cut_off_bnd=cut_off_bnd)
+
+    # Find indices where neighbor_count > 0
+    indices = findall(x -> x > 0, results.neighbor_count)
+
+    # Filter all arrays in the named tuple using these indices
+    filtered_results = map(x -> x[indices], results)
+
+    return filtered_results
+end
+
+@doc raw"""
+    interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, sol;
+                         smoothing_length=ref_system.smoothing_length, cut_off_bnd=true)
+
+Interpolates properties along a plane in a 3D space in a TrixiParticles simulation.
+The plane for interpolation is defined by three points in 3D space,
+with a specified resolution determining the density of the interpolation points.
+
+The function generates a grid of points on a parallelogram within the plane defined by the three points, spaced uniformly according to the given resolution.
+
+# Arguments
+- `point1`: The first point defining the plane.
+- `point2`: The second point defining the plane.
+- `point3`: The third point defining the plane. The points must not be collinear.
+- `resolution`: The distance between adjacent interpolation points in the grid.
+- `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`: 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`.
+- `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 within the plane.
+
+!!! note
+    - 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 across a plane defined by points [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], and [0.0, 1.0, 0.0]
+# with a resolution of 0.1
+results = interpolate_plane_3d([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], 0.1, semi, ref_system, sol)
+```
+"""
+function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, sol;
+                              smoothing_length=ref_system.smoothing_length,
+                              cut_off_bnd=true)
+    # Verify that points are in 3D space
+    if length(point1) != 3 || length(point2) != 3 || length(point3) != 3
+        throw(ArgumentError("all points must be 3D coordinates"))
+    end
+
+    point1_ = SVector{3}(point1)
+    point2_ = SVector{3}(point2)
+    point3_ = SVector{3}(point3)
+
+    # Vectors defining the edges of the parallelogram
+    edge1 = point2_ - point1_
+    edge2 = point3_ - point1_
+
+    # Check if the points are collinear
+    if norm(cross(edge1, edge2)) == 0
+        throw(ArgumentError("the points must not be collinear"))
+    end
+
+    # Determine the number of points along each edge
+    num_points_edge1 = ceil(Int, norm(edge1) / resolution)
+    num_points_edge2 = ceil(Int, norm(edge2) / resolution)
+
+    # Create a set of points on the plane
+    points_coords = Vector{SVector{3, Float64}}(undef,
+                                                (num_points_edge1 + 1) *
+                                                (num_points_edge2 + 1))
+    index = 1
+    for i in 0:num_points_edge1
+        for j in 0:num_points_edge2
+            point_on_plane = point1 + (i / num_points_edge1) * edge1 +
+                             (j / num_points_edge2) * edge2
+            points_coords[index] = point_on_plane
+            index += 1
+        end
+    end
+
+    # Interpolate using the generated points
+    results = interpolate_point(points_coords, semi, ref_system, sol,
+                                smoothing_length=smoothing_length,
+                                cut_off_bnd=cut_off_bnd)
+
+    # Filter results
+    indices = findall(x -> x > 0, results.neighbor_count)
+    filtered_results = map(x -> x[indices], results)
+
+    return filtered_results
+end
+
 @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.
+Interpolates properties along a line in a 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,

diff --git a/test/Project.toml b/test/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"

diff --git a/test/examples/examples.jl b/test/examples/examples.jl
@@ -18,6 +18,22 @@
             @test sol.retcode == ReturnCode.Success
         end
 
+        @trixi_testset "fluid/rectangular_tank_3d.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid",
+                                                "rectangular_tank_3d.jl"), tspan=(0.0, 0.1))
+            @test sol.retcode == ReturnCode.Success
+        end
+
+        @trixi_testset "fluid/rectangular_tank_3d.jl with SummationDensity" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "fluid",
+                                                "rectangular_tank_3d.jl"), tspan=(0.0, 0.1),
+                                       fluid_density_calculator=SummationDensity(),
+                                       clip_negative_pressure=true)
+            @test sol.retcode == ReturnCode.Success
+        end
+
         @trixi_testset "fluid/rectangular_tank_edac_2d.jl" begin
             @test_nowarn trixi_include(@__MODULE__,
                                        joinpath(examples_dir(), "fluid",
@@ -136,4 +152,20 @@
                                                 "n_body_benchmark_reference_faster.jl"))
         end
     end
+
+    @testset verbose=true "Postprocessing" begin
+        @trixi_testset "postprocessing/interpolation_plane.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "postprocessing",
+                                                "interpolation_plane.jl"),
+                                       tspan=(0.0, 0.01))
+            @test sol.retcode == ReturnCode.Success
+        end
+        @trixi_testset "postprocessing/interpolation_point_line.jl" begin
+            @test_nowarn trixi_include(@__MODULE__,
+                                       joinpath(examples_dir(), "postprocessing",
+                                                "interpolation_point_line.jl"))
+            @test sol.retcode == ReturnCode.Success
+        end
+    end
 end