Skip to content

Add variable smoothing length #736

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
svchb merged 28 commits into from
Apr 5, 2025
Merged

Add variable smoothing length #736

Show file tree
Hide file tree
Changes from 25 commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
67647b5
adaptive smoothing length
Mar 17, 2025
9e2ac47
fix unit tests
Mar 17, 2025
3df9b54
fix tests
Mar 17, 2025
7266c4a
fix tests again
Mar 17, 2025
9d58ec6
fix correction tests
Mar 26, 2025
87c5409
fix tests
Mar 30, 2025
aaa3163
fix
Mar 30, 2025
fa9c9ee
fix
Mar 30, 2025
a698d25
fix packing
Mar 30, 2025
d937bc6
fix
Mar 31, 2025
fa8d8a9
fix
Apr 1, 2025
6731151
comment `particle_spacing`
Apr 1, 2025
dd44a73
implement suggestions
Apr 1, 2025
dc028e7
Merge branch 'main' into prepare-multiresolution
Apr 1, 2025
d87382e
fix
Apr 1, 2025
bdaacfa
Merge branch 'main' into prepare-multiresolution
Apr 3, 2025
4dd22da
fix tests
Apr 3, 2025
ba6f08b
fix typo
Apr 3, 2025
fd61e4d
Merge branch 'main' into prepare-multiresolution
Apr 4, 2025
923fabc
add `particle` to pressure acceleration
Apr 4, 2025
efc8758
add suggestions
Apr 4, 2025
ddf4f96
fix
Apr 4, 2025
c34d234
rm
Apr 4, 2025
b1b8411
implement suggestions
Apr 4, 2025
4278290
apply formatter
Apr 4, 2025
be4e96d
Merge branch 'main' into prepare-multiresolution
Apr 4, 2025
76fea69
Merge branch 'main' into prepare-multiresolution
Apr 5, 2025
8cb5107
rm TODO
Apr 5, 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
2 changes: 1 addition & 1 deletion examples/solid/oscillating_beam_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ semi = Semidiscretization(solid_system,
neighborhood_search=PrecomputedNeighborhoodSearch{2}())
ode = semidiscretize(semi, tspan, data_type=nothing)

info_callback = InfoCallback(interval=100)
info_callback = InfoCallback(interval=1000)

# Track the position of the particle in the middle of the tip of the beam.
middle_particle_id = Int(n_particles_per_dimension[1] * (n_particles_per_dimension[2] + 1) /
Expand Down
29 changes: 18 additions & 11 deletions src/general/corrections.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,8 @@ function compute_shepard_coeff!(system, system_coords, v_ode, u_ode, semi,
volume = m_b / rho_b

kernel_correction_coefficient[particle] += volume *
smoothing_kernel(system, distance)
smoothing_kernel(system, distance,
particle)
end
end

Expand Down Expand Up @@ -217,10 +218,15 @@ function compute_correction_values!(system,
m_b = hydrodynamic_mass(neighbor_system, neighbor)
volume = m_b / rho_b

kernel_correction_coefficient[particle] += volume *
smoothing_kernel(system, distance)
# Use uncorrected kernel to compute correction coefficients
W = kernel(system_smoothing_kernel(system), distance,
smoothing_length(system, particle))

kernel_correction_coefficient[particle] += volume * W
if distance > sqrt(eps())
tmp = volume * smoothing_kernel_grad(system, pos_diff, distance)
grad_W = kernel_grad(system_smoothing_kernel(system), pos_diff, distance,
smoothing_length(system, particle))
tmp = volume * grad_W
for i in axes(dw_gamma, 1)
dw_gamma[i, particle] += tmp[i]
end
Expand Down Expand Up @@ -311,7 +317,7 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search,
pos_diff, distance
volume = mass[neighbor] / density_fun(neighbor)

grad_kernel = smoothing_kernel_grad(system, pos_diff, distance)
grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)

iszero(grad_kernel) && return

Expand All @@ -329,7 +335,7 @@ end

function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
coordinates, v_ode, u_ode, semi,
correction, smoothing_length, smoothing_kernel)
correction, smoothing_kernel)
set_zero!(corr_matrix)

# Loop over all pairs of particles and neighbors within the kernel cutoff
Expand All @@ -345,23 +351,24 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
pos_diff, distance
volume = hydrodynamic_mass(neighbor_system, neighbor) /
particle_density(v_neighbor_system, neighbor_system, neighbor)
smoothing_length_ = smoothing_length(system, particle)

function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance,
smoothing_length, system, particle)
return smoothing_kernel_grad(system, pos_diff, distance)
smoothing_length_, system, particle)
return smoothing_kernel_grad(system, pos_diff, distance, particle)
end

# Compute gradient of corrected kernel
function compute_grad_kernel(correction::MixedKernelGradientCorrection,
smoothing_kernel, pos_diff, distance,
smoothing_length, system, particle)
smoothing_length_, system, particle)
return corrected_kernel_grad(smoothing_kernel, pos_diff, distance,
smoothing_length, KernelCorrection(), system,
smoothing_length_, KernelCorrection(), system,
particle)
end

grad_kernel = compute_grad_kernel(correction, smoothing_kernel, pos_diff,
distance, smoothing_length, system, particle)
distance, smoothing_length_, system, particle)

iszero(grad_kernel) && return

Expand Down
2 changes: 1 addition & 1 deletion src/general/density_calculators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ function summation_density!(system, semi, u, u_ode, density;
points=particles) do particle, neighbor,
pos_diff, distance
mass = hydrodynamic_mass(neighbor_system, neighbor)
density[particle] += mass * smoothing_kernel(system, distance)
density[particle] += mass * smoothing_kernel(system, distance, particle)
end
end
end
51 changes: 27 additions & 24 deletions src/general/interpolation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ 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,
smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true,
clip_negative_pressure=false)

Interpolates properties along a plane in a TrixiParticles simulation.
Expand All @@ -24,7 +24,7 @@ See also: [`interpolate_plane_2d_vtk`](@ref), [`interpolate_plane_3d`](@ref),
- `sol`: The solution state from which the properties are interpolated.

# Keywords
- `smoothing_length=ref_system.smoothing_length`: The smoothing length used in the interpolation.
- `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation.
- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point
is "closer" to the boundary than to the fluid in a kernel-weighted sense.
Or, in more detail, when the boundary has more influence than the fluid
Expand Down Expand Up @@ -53,7 +53,7 @@ results = interpolate_plane_2d([0.0, 0.0], [1.0, 1.0], 0.2, semi, ref_system, so
"""
function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system,
sol::ODESolution;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
# Filter out particles without neighbors
filter_no_neighbors = true
Expand All @@ -69,7 +69,7 @@ end

function interpolate_plane_2d(min_corner, max_corner, resolution, semi, ref_system,
v_ode, u_ode;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
# Filter out particles without neighbors
filter_no_neighbors = true
Expand All @@ -84,7 +84,7 @@ end

@doc raw"""
interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system, sol;
smoothing_length=ref_system.smoothing_length, cut_off_bnd=true,
smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true,
clip_negative_pressure=false, output_directory="out", filename="plane")

Interpolates properties along a plane in a TrixiParticles simulation and exports the result
Expand All @@ -107,7 +107,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_3d`](@ref),
- `sol`: The solution state from which the properties are interpolated.

# Keywords
- `smoothing_length=ref_system.smoothing_length`: The smoothing length used in the interpolation.
- `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation.
- `output_directory="out"`: Directory to save the VTI file.
- `filename="plane"`: Name of the VTI file.
- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point
Expand All @@ -132,7 +132,7 @@ results = interpolate_plane_2d([0.0, 0.0], [1.0, 1.0], 0.2, semi, ref_system, so
"""
function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system,
sol::ODESolution; clip_negative_pressure=false,
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true,
output_directory="out", filename="plane")
v_ode, u_ode = sol.u[end].x
Expand All @@ -144,7 +144,7 @@ end

function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_system,
v_ode, u_ode;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false,
output_directory="out", filename="plane")
# Don't filter out particles without neighbors to keep 2D grid structure
Expand Down Expand Up @@ -206,7 +206,7 @@ end

@doc raw"""
interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system, sol;
smoothing_length=ref_system.smoothing_length, cut_off_bnd=true,
smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true,
clip_negative_pressure=false)

Interpolates properties along a plane in a 3D space in a TrixiParticles simulation.
Expand All @@ -229,7 +229,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_2d_vtk`](@ref),
- `sol`: The solution state from which the properties are interpolated.

# Keywords
- `smoothing_length=ref_system.smoothing_length`: The smoothing length used in the interpolation.
- `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation.
- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point
is "closer" to the boundary than to the fluid in a kernel-weighted sense.
Or, in more detail, when the boundary has more influence than the fluid
Expand Down Expand Up @@ -259,7 +259,7 @@ results = interpolate_plane_3d([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]
"""
function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system,
sol::ODESolution;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
v_ode, u_ode = sol.u[end].x

Expand All @@ -269,7 +269,8 @@ function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_syst
end

function interpolate_plane_3d(point1, point2, point3, resolution, semi, ref_system,
v_ode, u_ode; smoothing_length=ref_system.smoothing_length,
v_ode, u_ode;
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
if ndims(ref_system) != 3
throw(ArgumentError("`interpolate_plane_3d` requires a 3D simulation"))
Expand Down Expand Up @@ -299,7 +300,7 @@ end

@doc raw"""
interpolate_line(start, end_, n_points, semi, ref_system, sol; endpoint=true,
smoothing_length=ref_system.smoothing_length, cut_off_bnd=true,
smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true,
clip_negative_pressure=false)

Interpolates properties along a line in a TrixiParticles simulation.
Expand All @@ -321,7 +322,7 @@ See also: [`interpolate_point`](@ref), [`interpolate_plane_2d`](@ref),

# Keywords
- `endpoint=true`: A boolean to include (`true`) or exclude (`false`) the end point in the interpolation.
- `smoothing_length=ref_system.smoothing_length`: The smoothing length used in the interpolation.
- `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation.
- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point
is "closer" to the boundary than to the fluid in a kernel-weighted sense.
Or, in more detail, when the boundary has more influence than the fluid
Expand Down Expand Up @@ -351,15 +352,17 @@ 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::ODESolution;
endpoint=true, smoothing_length=ref_system.smoothing_length,
endpoint=true,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
v_ode, u_ode = sol.u[end].x

interpolate_line(start, end_, n_points, semi, ref_system, v_ode, u_ode;
endpoint, smoothing_length, cut_off_bnd, clip_negative_pressure)
end
function interpolate_line(start, end_, n_points, semi, ref_system, v_ode, u_ode;
endpoint=true, smoothing_length=ref_system.smoothing_length,
endpoint=true,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
start_svector = SVector{ndims(ref_system)}(start)
end_svector = SVector{ndims(ref_system)}(end_)
Expand All @@ -376,11 +379,11 @@ end

@doc raw"""
interpolate_point(points_coords::Array{Array{Float64,1},1}, semi, ref_system, sol;
smoothing_length=ref_system.smoothing_length, cut_off_bnd=true,
smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true,
clip_negative_pressure=false)

interpolate_point(point_coords, semi, ref_system, sol;
smoothing_length=ref_system.smoothing_length, cut_off_bnd=true,
smoothing_length=initial_smoothing_length(ref_system), cut_off_bnd=true,
clip_negative_pressure=false)

Performs interpolation of properties at specified points or an array of points in a TrixiParticles simulation.
Expand All @@ -400,7 +403,7 @@ See also: [`interpolate_line`](@ref), [`interpolate_plane_2d`](@ref),
- `sol`: The current solution state from which properties are interpolated.

# Keywords
- `smoothing_length=ref_system.smoothing_length`: The smoothing length used in the interpolation.
- `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation.
- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point
is "closer" to the boundary than to the fluid in a kernel-weighted sense.
Or, in more detail, when the boundary has more influence than the fluid
Expand Down Expand Up @@ -434,7 +437,7 @@ results = interpolate_point(points, semi, ref_system, sol)
accurate as a real surface reconstruction.
"""
@inline function interpolate_point(point_coords, semi, ref_system, sol::ODESolution;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
v_ode, u_ode = sol.u[end].x

Expand All @@ -444,7 +447,7 @@ end

@inline function interpolate_point(points_coords::AbstractArray{<:AbstractArray}, semi,
ref_system, v_ode, u_ode;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
num_points = length(points_coords)
coords = similar(points_coords)
Expand Down Expand Up @@ -472,7 +475,7 @@ end
end

function interpolate_point(point_coords, semi, ref_system, v_ode, u_ode;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
neighborhood_searches = process_neighborhood_searches(semi, u_ode, ref_system,
smoothing_length)
Expand All @@ -483,7 +486,7 @@ function interpolate_point(point_coords, semi, ref_system, v_ode, u_ode;
end

function process_neighborhood_searches(semi, u_ode, ref_system, smoothing_length)
if isapprox(smoothing_length, ref_system.smoothing_length)
if isapprox(smoothing_length, initial_smoothing_length(ref_system))
# Update existing NHS
update_nhs!(semi, u_ode)
neighborhood_searches = semi.neighborhood_searches[system_indices(ref_system, semi)]
Expand All @@ -507,7 +510,7 @@ end

@inline function interpolate_point(point_coords, semi, ref_system, v_ode, u_ode,
neighborhood_searches;
smoothing_length=ref_system.smoothing_length,
smoothing_length=initial_smoothing_length(ref_system),
cut_off_bnd=true, clip_negative_pressure=false)
interpolated_density = 0.0
other_density = 0.0
Expand Down
5 changes: 3 additions & 2 deletions src/general/semidiscretization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,8 +134,9 @@ function create_neighborhood_search(neighborhood_search, system, neighbor)
end

@inline function compact_support(system, neighbor)
(; smoothing_kernel, smoothing_length) = system
return compact_support(smoothing_kernel, smoothing_length)
(; smoothing_kernel) = system
# TODO: Variable search radius for NHS?
return compact_support(smoothing_kernel, initial_smoothing_length(system))
end

@inline function compact_support(system::OpenBoundarySPHSystem, neighbor)
Expand Down
Loading
Loading