Skip to content

Separate computed and interpolated density in interpolation #808

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

Merged
efaulhaber merged 8 commits into from
May 18, 2025
Merged

Separate computed and interpolated density in interpolation #808

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
afd786a
fix interpolation
May 15, 2025
58f8544
add weighting per mass
May 15, 2025
192713b
fix tests
May 15, 2025
027e36c
add volume contribution
May 15, 2025
42bb3dd
fix doc test
May 15, 2025
a78ca66
fix comments
May 15, 2025
a2a4a19
Merge branch 'main' into fix-interpolation
May 15, 2025
53a53ec
Merge branch 'main' into fix-interpolation
LasNikas May 15, 2025
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
63 changes: 35 additions & 28 deletions src/general/interpolation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ See also: [`interpolate_plane_2d_vtk`](@ref), [`interpolate_plane_3d`](@ref),
results = interpolate_plane_2d([0.0, 0.0], [1.0, 1.0], 0.2, semi, ref_system, sol)

# output
(density = ...)
(computed_density = ...)
```
"""
function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system,
Expand Down Expand Up @@ -267,7 +267,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_2d_vtk`](@ref),
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)

# output
(density = ...)
(computed_density = ...)
```
"""
function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system,
Expand Down Expand Up @@ -365,7 +365,7 @@ See also: [`interpolate_points`](@ref), [`interpolate_plane_2d`](@ref),
results = interpolate_line([1.0, 0.0], [1.0, 1.0], 5, semi, ref_system, sol)

# output
(density = ...)
(computed_density = ...)
```
"""
function interpolate_line(start, end_, n_points, semi, ref_system, sol::ODESolution;
Expand Down Expand Up @@ -441,7 +441,7 @@ points = [1.0 1.0 1.0; 0.5 0.6 0.7]
results = interpolate_points(points, semi, ref_system, sol)

# output
(density = ...)
(computed_density = ...)
```
!!! note
- This function is particularly useful for analyzing gradients or creating visualizations
Expand Down Expand Up @@ -512,10 +512,10 @@ end

n_points = size(point_coords, 2)
ELTYPE = eltype(point_coords)
interpolated_density = zeros(ELTYPE, n_points)
computed_density = zeros(ELTYPE, n_points)
other_density = zeros(ELTYPE, n_points)
shepard_coefficient = zeros(ELTYPE, n_points)
neighbor_count = zeros(Int, n_points)
shepard_coefficient = zeros(ELTYPE, n_points)

cache = create_cache_interpolation(ref_system, n_points)

Expand All @@ -537,29 +537,32 @@ end
foreach_point_neighbor(point_coords, neighbor_coords, nhs;
parallelization_backend) do point, neighbor, pos_diff,
distance
m_a = hydrodynamic_mass(neighbor_system, neighbor)
W_a = kernel(ref_smoothing_kernel, distance, smoothing_length)
m_b = hydrodynamic_mass(neighbor_system, neighbor)
volume_b = m_b / current_density(v, neighbor_system, neighbor)
W_ab = kernel(ref_smoothing_kernel, distance, smoothing_length)

if system_id == ref_id
interpolated_density[point] += m_a * W_a
volume = m_a / current_density(v, neighbor_system, neighbor)
shepard_coefficient[point] += volume * W_a
computed_density[point] += m_b * W_ab
shepard_coefficient[point] += volume_b * W_ab

# According to:
# u(r_a) = (∑_b u(r_b) ⋅ V_b ⋅ W(r_a-r_b)) / (∑_b V_b ⋅ W(r_a-r_b)),
# where V_b = m_b / ρ_b.
interpolate_system!(cache, v, neighbor_system,
point, neighbor, volume, W_a, clip_negative_pressure)
point, neighbor, volume_b, W_ab, clip_negative_pressure)
else
other_density[point] += m_a * W_a
other_density[point] += m_b * W_ab
end

neighbor_count[point] += 1
end
end

@threaded parallelization_backend for point in axes(point_coords, 2)
if other_density[point] > interpolated_density[point] ||
shepard_coefficient[point] < eps()
if other_density[point] > computed_density[point] ||
computed_density[point] < eps()
# Return NaN values that can be filtered out in ParaView
interpolated_density[point] = NaN
computed_density[point] = NaN
neighbor_count[point] = 0

foreach(cache) do field
Expand All @@ -571,8 +574,6 @@ end
end
else
# Normalize all quantities by the shepard coefficient
interpolated_density[point] /= shepard_coefficient[point]

foreach(cache) do field
if field isa AbstractVector
field[point] /= shepard_coefficient[point]
Expand All @@ -583,14 +584,15 @@ end
end
end

return (; density=interpolated_density, neighbor_count, point_coords, cache...)
return (; computed_density=computed_density, neighbor_count, point_coords, cache...)
end

@inline function create_cache_interpolation(ref_system::FluidSystem, n_points)
velocity = zeros(eltype(ref_system), ndims(ref_system), n_points)
pressure = zeros(eltype(ref_system), n_points)
density = zeros(eltype(ref_system), n_points)

return (; velocity, pressure)
return (; velocity, pressure, density)
end

@inline function create_cache_interpolation(ref_system::SolidSystem, n_points)
Expand All @@ -604,34 +606,39 @@ end
end

@inline function interpolate_system!(cache, v, system::FluidSystem,
point, neighbor, volume, W_a, clip_negative_pressure)
point, neighbor, volume_b, W_ab,
clip_negative_pressure)
velocity = current_velocity(v, system, neighbor)
for i in axes(cache.velocity, 1)
cache.velocity[i, point] += velocity[i] * (volume * W_a)
cache.velocity[i, point] += velocity[i] * volume_b * W_ab
end

pressure = current_pressure(v, system, neighbor)
if clip_negative_pressure
pressure = max(zero(eltype(pressure)), pressure)
end
cache.pressure[point] += pressure * (volume * W_a)
cache.pressure[point] += pressure * volume_b * W_ab

density = current_density(v, system, neighbor)
cache.density[point] += density * volume_b * W_ab

return cache
end

@inline function interpolate_system!(cache, v, system::SolidSystem,
point, neighbor, volume, W_a, clip_negative_pressure)
point, neighbor, volume_b, W_ab,
clip_negative_pressure)
velocity = current_velocity(v, system, neighbor)
for i in axes(cache.velocity, 1)
cache.velocity[i, point] += velocity[i] * (volume * W_a)
cache.velocity[i, point] += velocity[i] * volume_b * W_ab
end

cache.jacobian[point] += det(deformation_gradient(system, neighbor)) * (volume * W_a)
cache.von_mises_stress[point] += von_mises_stress(system) * (volume * W_a)
cache.jacobian[point] += det(deformation_gradient(system, neighbor)) * volume_b * W_ab
cache.von_mises_stress[point] += von_mises_stress(system) * volume_b * W_ab

sigma = cauchy_stress(system)
for j in axes(cache.cauchy_stress, 2), i in axes(cache.cauchy_stress, 1)
cache.cauchy_stress[i, j, point] += sigma[i, j, neighbor] * (volume * W_a)
cache.cauchy_stress[i, j, point] += sigma[i, j, neighbor] * volume_b * W_ab
end

return cache
Expand Down
8 changes: 4 additions & 4 deletions test/general/interpolation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -412,7 +412,7 @@
if cut_off_bnd
@test isapprox(distance_bottom_outside, 0, atol=1e-12)
else
@test isapprox(distance_bottom_outside, 0.11817145975510357, atol=1e-14)
@test isapprox(distance_bottom_outside, 0.11817392366574496, atol=1e-14)
end

result_bottom_outside = interpolation_walldistance(-0.5 *
Expand Down Expand Up @@ -664,7 +664,7 @@
(ny + 2) * particle_spacing,
interpolation_walldistance,
tolerance=1e-8)
@test isapprox(distance_top_outside, 0.09390581846237112, atol=1e-14)
@test isapprox(distance_top_outside, 0.09390704035758901, atol=1e-14)

result_zero = interpolation_walldistance(ny * particle_spacing +
distance_top_outside)
Expand Down Expand Up @@ -705,7 +705,7 @@
distance_bottom_outside = binary_search_outside(0.0, -2 * particle_spacing,
interpolation_walldistance,
tolerance=1e-12)
@test isapprox(distance_bottom_outside, 0.09390581390689477, atol=1e-14)
@test isapprox(distance_bottom_outside, 0.0939070362292114, atol=1e-14)

exp_res = (density=[666.0],
neighbor_count=[4],
Expand Down Expand Up @@ -805,7 +805,7 @@
if cut_off_bnd
@test isapprox(distance_bottom_outside, 0, atol=1e-12)
else
@test isapprox(distance_bottom_outside, 0.09390581390689477, atol=1e-14)
@test isapprox(distance_bottom_outside, 0.0939070362292114, atol=1e-14)
end

result_bottom_outside = interpolation_walldistance(-0.5 *
Expand Down
Loading