From afd786ace912baf166e9a0ec670e2857442d6422 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 May 2025 11:27:18 +0200 Subject: [PATCH 1/6] fix interpolation --- src/general/interpolation.jl | 47 ++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index 15ffd24003..fd73bbaf73 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -512,9 +512,9 @@ 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) + shepard_weighting = zeros(ELTYPE, n_points) neighbor_count = zeros(Int, n_points) cache = create_cache_interpolation(ref_system, n_points) @@ -541,12 +541,11 @@ end W_a = 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_a * W_a + shepard_weighting[point] += W_a interpolate_system!(cache, v, neighbor_system, - point, neighbor, volume, W_a, clip_negative_pressure) + point, neighbor, W_a, clip_negative_pressure) else other_density[point] += m_a * W_a end @@ -556,10 +555,10 @@ 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] || + shepard_weighting[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 @@ -570,27 +569,26 @@ end end end else - # Normalize all quantities by the shepard coefficient - interpolated_density[point] /= shepard_coefficient[point] - + # Normalize all quantities by the shepard weighting foreach(cache) do field if field isa AbstractVector - field[point] /= shepard_coefficient[point] + field[point] /= shepard_weighting[point] else - field[:, point] ./= shepard_coefficient[point] + field[:, point] ./= shepard_weighting[point] end 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) @@ -604,17 +602,20 @@ end end @inline function interpolate_system!(cache, v, system::FluidSystem, - point, neighbor, volume, W_a, clip_negative_pressure) + point, neighbor, W_a, 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] * W_a 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 * W_a + + density = current_density(v, system, neighbor) + cache.density[point] += density * W_a return cache end @@ -623,15 +624,15 @@ end point, neighbor, volume, W_a, 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] * W_a 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)) * W_a + cache.von_mises_stress[point] += von_mises_stress(system) * W_a 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] * W_a end return cache From 58f8544300586c9c0778a26e31b80a75cc7c92df Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 May 2025 12:40:17 +0200 Subject: [PATCH 2/6] add weighting per mass --- src/general/interpolation.jl | 38 ++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index fd73bbaf73..dc86c119ee 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -514,7 +514,6 @@ end ELTYPE = eltype(point_coords) computed_density = zeros(ELTYPE, n_points) other_density = zeros(ELTYPE, n_points) - shepard_weighting = zeros(ELTYPE, n_points) neighbor_count = zeros(Int, n_points) cache = create_cache_interpolation(ref_system, n_points) @@ -537,17 +536,18 @@ 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) + # Weighting mass + m_b = hydrodynamic_mass(neighbor_system, neighbor) + W_ab = kernel(ref_smoothing_kernel, distance, smoothing_length) if system_id == ref_id - computed_density[point] += m_a * W_a - shepard_weighting[point] += W_a + # The computed density is the shepard weighting + computed_density[point] += m_b * W_ab interpolate_system!(cache, v, neighbor_system, - point, neighbor, W_a, clip_negative_pressure) + point, neighbor, m_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 @@ -556,7 +556,7 @@ end @threaded parallelization_backend for point in axes(point_coords, 2) if other_density[point] > computed_density[point] || - shepard_weighting[point] < eps() + computed_density[point] < eps() # Return NaN values that can be filtered out in ParaView computed_density[point] = NaN neighbor_count[point] = 0 @@ -572,9 +572,9 @@ end # Normalize all quantities by the shepard weighting foreach(cache) do field if field isa AbstractVector - field[point] /= shepard_weighting[point] + field[point] /= computed_density[point] else - field[:, point] ./= shepard_weighting[point] + field[:, point] ./= computed_density[point] end end end @@ -602,37 +602,37 @@ end end @inline function interpolate_system!(cache, v, system::FluidSystem, - point, neighbor, W_a, clip_negative_pressure) + point, neighbor, m_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] * W_a + cache.velocity[i, point] += velocity[i] * W_ab * m_b end pressure = current_pressure(v, system, neighbor) if clip_negative_pressure pressure = max(zero(eltype(pressure)), pressure) end - cache.pressure[point] += pressure * W_a + cache.pressure[point] += pressure * W_ab * m_b density = current_density(v, system, neighbor) - cache.density[point] += density * W_a + cache.density[point] += density * W_ab * m_b return cache end @inline function interpolate_system!(cache, v, system::SolidSystem, - point, neighbor, volume, W_a, clip_negative_pressure) + point, neighbor, m_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] * W_a + cache.velocity[i, point] += velocity[i] * W_ab * m_b end - cache.jacobian[point] += det(deformation_gradient(system, neighbor)) * W_a - cache.von_mises_stress[point] += von_mises_stress(system) * W_a + cache.jacobian[point] += det(deformation_gradient(system, neighbor)) * W_ab * m_b + cache.von_mises_stress[point] += von_mises_stress(system) * W_ab * m_b 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] * W_a + cache.cauchy_stress[i, j, point] += sigma[i, j, neighbor] * W_ab * m_b end return cache From 192713b34f561cfde4bba4b0197d325ffc8c1b09 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 May 2025 14:01:28 +0200 Subject: [PATCH 3/6] fix tests --- test/general/interpolation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/general/interpolation.jl b/test/general/interpolation.jl index 8624b18041..b7ee61a66b 100644 --- a/test/general/interpolation.jl +++ b/test/general/interpolation.jl @@ -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 * @@ -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) @@ -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], @@ -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 * From 027e36cce785e14b51af7ef7f5bd1d4e1d7ca3e8 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 May 2025 15:58:10 +0200 Subject: [PATCH 4/6] add volume contribution --- src/general/interpolation.jl | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index dc86c119ee..ad6c2e91c7 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -515,6 +515,7 @@ end computed_density = zeros(ELTYPE, n_points) other_density = zeros(ELTYPE, n_points) neighbor_count = zeros(Int, n_points) + shepard_coefficient = zeros(ELTYPE, n_points) cache = create_cache_interpolation(ref_system, n_points) @@ -538,14 +539,18 @@ end distance # Weighting mass 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 - # The computed density is the shepard weighting 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, m_b, W_ab, clip_negative_pressure) + point, neighbor, volume_b, W_ab, clip_negative_pressure) else other_density[point] += m_b * W_ab end @@ -572,9 +577,9 @@ end # Normalize all quantities by the shepard weighting foreach(cache) do field if field isa AbstractVector - field[point] /= computed_density[point] + field[point] /= shepard_coefficient[point] else - field[:, point] ./= computed_density[point] + field[:, point] ./= shepard_coefficient[point] end end end @@ -602,37 +607,39 @@ end end @inline function interpolate_system!(cache, v, system::FluidSystem, - point, neighbor, m_b, W_ab, 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] * W_ab * m_b + 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 * W_ab * m_b + cache.pressure[point] += pressure * volume_b * W_ab density = current_density(v, system, neighbor) - cache.density[point] += density * W_ab * m_b + cache.density[point] += density * volume_b * W_ab return cache end @inline function interpolate_system!(cache, v, system::SolidSystem, - point, neighbor, m_b, W_ab, 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] * W_ab * m_b + cache.velocity[i, point] += velocity[i] * volume_b * W_ab end - cache.jacobian[point] += det(deformation_gradient(system, neighbor)) * W_ab * m_b - cache.von_mises_stress[point] += von_mises_stress(system) * W_ab * m_b + 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] * W_ab * m_b + cache.cauchy_stress[i, j, point] += sigma[i, j, neighbor] * volume_b * W_ab end return cache From 42bb3dd2806811be16c56cd519796b055634765d Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 May 2025 16:21:47 +0200 Subject: [PATCH 5/6] fix doc test --- src/general/interpolation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index ad6c2e91c7..6a0c57a177 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -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, @@ -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, @@ -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; @@ -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 From a78ca66a37b96f2d43c9d45c805c572146a877c0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 15 May 2025 17:57:34 +0200 Subject: [PATCH 6/6] fix comments --- src/general/interpolation.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index 6a0c57a177..a17cddddba 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -537,7 +537,6 @@ end foreach_point_neighbor(point_coords, neighbor_coords, nhs; parallelization_backend) do point, neighbor, pos_diff, distance - # Weighting mass 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) @@ -574,7 +573,7 @@ end end end else - # Normalize all quantities by the shepard weighting + # Normalize all quantities by the shepard coefficient foreach(cache) do field if field isa AbstractVector field[point] /= shepard_coefficient[point]