Skip to content

Fix interpolation edge cases #1183

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
svchb wants to merge 1 commit into
base: main
Choose a base branch
from
+70 −23
Open

Fix interpolation edge cases #1183

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 37 additions & 23 deletions src/general/interpolation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,10 +198,9 @@ function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_syst
x_range = range(min_corner[1], max_corner[1], length=n_points_per_dimension[1])
y_range = range(min_corner[2], max_corner[2], length=n_points_per_dimension[2])

# Generate points within the plane. Use `place_on_shell=true` to generate points
# on the shell of the geometry.
point_coords = rectangular_shape_coords(resolution, n_points_per_dimension, min_corner,
place_on_shell=true)
# Generate points from the exact ranges used for VTK output. This keeps interpolation
# points inside the requested box even when `resolution` does not divide its side lengths.
point_coords = plane_point_coords(x_range, y_range)

results = interpolate_points(point_coords, semi, ref_system, v_ode, u_ode;
smoothing_length, cut_off_bnd, include_wall_velocity,
Expand All @@ -211,14 +210,8 @@ function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_syst
# Find indices where neighbor_count > 0
indices = findall(x -> x > 0, results.neighbor_count)

# Filter all arrays in the named tuple using these indices
results = map(results) do x
if isa(x, AbstractVector)
return x[indices]
else
return x[:, indices]
end
end
# Filter all arrays in the named tuple using these indices.
results = filter_interpolation_results(results, indices)
end

return results, x_range, y_range
Expand Down Expand Up @@ -309,13 +302,7 @@ function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_syst

# Filter results
indices = findall(x -> x > 0, results.neighbor_count)
filtered_results = map(results) do x
if isa(x, AbstractVector)
return x[indices]
else
return x[:, indices]
end
end
filtered_results = filter_interpolation_results(results, indices)

return filtered_results
end
Expand Down Expand Up @@ -513,7 +500,8 @@ function process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length
old_nhs = get_neighborhood_search(ref_system, system, semi)
nhs = PointNeighbors.copy_neighborhood_search(old_nhs, search_radius,
nparticles(system))
PointNeighbors.initialize!(nhs, point_coords, system_coords)
PointNeighbors.initialize!(nhs, point_coords, system_coords,
eachindex_y=each_active_particle(system))

return nhs
end
Expand Down Expand Up @@ -555,8 +543,9 @@ end
ref_id = system_indices(ref_system, semi)
ref_smoothing_kernel = ref_system.smoothing_kernel

# If we don't cut at the boundary, we only need to iterate over the reference system
systems = cut_off_bnd ? semi : (ref_system,)
# If we don't cut at the boundary and don't include wall velocity, we only need to
# iterate over the reference system.
systems = (cut_off_bnd || include_wall_velocity) ? semi : (ref_system,)

foreach_system(systems) do neighbor_system
system_id = system_indices(neighbor_system, semi)
Expand Down Expand Up @@ -590,7 +579,9 @@ end
interpolate_system!(cache, v, neighbor_system,
point, neighbor, volume_b, W_ab, clip_negative_pressure)
else
other_density[point] += m_b * W_ab
if cut_off_bnd
other_density[point] += m_b * W_ab
end

if include_wall_velocity
velocity_neighbor_ = current_velocity(v, neighbor_system, neighbor)
Expand Down Expand Up @@ -640,6 +631,29 @@ end
return (; computed_density, point_coords, neighbor_count, cache...)
end

function plane_point_coords(x_range, y_range)
ELTYPE = promote_type(eltype(x_range), eltype(y_range))
point_coords = Array{ELTYPE, 2}(undef, 2, length(x_range) * length(y_range))

point = 1
for y in y_range, x in x_range
point_coords[1, point] = x
point_coords[2, point] = y
point += 1
end

return point_coords
end

function filter_interpolation_results(results, indices)
return map(field -> filter_interpolation_field(field, indices), results)
end

function filter_interpolation_field(field::AbstractArray, indices)
selectors = ntuple(dim -> dim == ndims(field) ? indices : Colon(), ndims(field))
return field[selectors...]
end

@inline function create_cache_interpolation(ref_system::AbstractFluidSystem, n_points, semi)
(; parallelization_backend) = semi

Expand Down
33 changes: 33 additions & 0 deletions test/general/interpolation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,31 @@
@testset verbose=true "SPH Interpolation" begin
@testset "Plane point coordinates" begin
x_range = range(0.0, 1.0, length=5)
y_range = range(2.0, 3.0, length=3)
point_coords = TrixiParticles.plane_point_coords(x_range, y_range)

@test point_coords[:, 1] == [0.0, 2.0]
@test point_coords[:, 5] == [1.0, 2.0]
@test point_coords[:, end] == [1.0, 3.0]
@test extrema(point_coords[1, :]) == (0.0, 1.0)
@test extrema(point_coords[2, :]) == (2.0, 3.0)
end

@testset "Filter tensor interpolation results" begin
results = (density=[10.0, 20.0, 30.0],
velocity=[1.0 2.0 3.0
4.0 5.0 6.0],
cauchy_stress=reshape(1:12, 2, 2, 3))
filtered = TrixiParticles.filter_interpolation_results(results, [1, 3])

@test filtered.density == [10.0, 30.0]
@test filtered.velocity == [1.0 3.0
4.0 6.0]
@test size(filtered.cauchy_stress) == (2, 2, 2)
@test filtered.cauchy_stress[:, :, 1] == results.cauchy_stress[:, :, 1]
@test filtered.cauchy_stress[:, :, 2] == results.cauchy_stress[:, :, 3]
end

function compare_interpolation_result(actual, expected; tolerance=5e-4)
@test length(actual.density) == length(expected.density)
for i in 1:length(expected.density)
Expand Down Expand Up @@ -578,8 +605,14 @@
v_no_wall_velocity = interpolate_points(points_coords, semi_boundary,
include_wall_velocity=false,
fluid_system, v_ode, u_ode).velocity
v_wall_velocity_without_cutoff = interpolate_points(points_coords, semi_boundary,
include_wall_velocity=true,
cut_off_bnd=false,
fluid_system, v_ode,
u_ode).velocity

@test isapprox(v_wall_velocity[2, 1], 0.0; atol=eps())
@test isapprox(v_wall_velocity_without_cutoff[2, 1], 0.0; atol=eps())
@test isapprox(v_no_wall_velocity[2, 1], 0.1; atol=eps())
@test any(isapprox.(v_wall_velocity[:, 3:end], v_no_wall_velocity[:, 3:end],
atol=eps()))
Expand Down
Loading