Skip to content
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.

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

Generalize surface normal calc #539

Open
wants to merge 228 commits into
base: main
Choose a base branch
from

Conversation

svchb
Copy link
Collaborator

@svchb svchb commented May 29, 2024

Depends on #599
extracted from #584

NEWS.md Outdated Show resolved Hide resolved
src/schemes/boundary/dummy_particles/dummy_particles.jl Outdated Show resolved Hide resolved
@svchb svchb requested a review from efaulhaber December 2, 2024 11:36
docs/src/systems/fluid.md Outdated Show resolved Hide resolved
initial_condition :: IC
mass :: MA # Array{ELTYPE, 1}
pressure :: P # Array{ELTYPE, 1}
density_calculator :: DC
state_equation :: SE
smoothing_kernel :: K
smoothing_length :: ELTYPE
ideal_neighbor_count :: Int
color :: Int
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is never used, is it? Same for EDAC.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes it will be used in #584

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As always, please remove it from this PR and add it when it's actually used.

test/schemes/fluid/surface_normal_sph.jl Outdated Show resolved Hide resolved
test/schemes/fluid/surface_normal_sph.jl Outdated Show resolved Hide resolved

# After computation, check that surface normals have been computed
@test all(isfinite.(system.cache.surface_normal))
@test all(isfinite.(system.cache.neighbor_count))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't this zero by default? Shouldn't you then check that this is not zero anymore?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The inside will stay at 0.0. though.

test/schemes/fluid/surface_normal_sph.jl Outdated Show resolved Hide resolved
test/schemes/fluid/surface_normal_sph.jl Show resolved Hide resolved
test/schemes/fluid/surface_normal_sph.jl Outdated Show resolved Hide resolved
Comment on lines 144 to 146
for i in surface_particles
@test isapprox(computed_normals[:, i], expected_normals[:, i], atol=0.05)
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
for i in surface_particles
@test isapprox(computed_normals[:, i], expected_normals[:, i], atol=0.05)
end
@test isapprox(computed_normals, expected_normals, atol=0.05)

The loop will generate a lot of tests.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't work since the normals are not always exactly 0.0 on the inside.

Copy link
Member

@efaulhaber efaulhaber Dec 9, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@test isapprox(computed_normals[:, surface_particles], expected_normals[:, surface_particles], atol=0.05)

Comment on lines +148 to +151
# Optionally, check that normals for interior particles are zero
# for i in setdiff(1:nparticles, surface_particles)
# @test isapprox(norm(system.cache.surface_normal[:, i]), 0.0, atol=1e-4)
# end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is also part of #584

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought I used it there but this is rather tricky... Even with the special conditions that are implemented in #584

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And with a larger tolerance? There should be some test to ensure the inner values are not completely random.
And there should not be a commented out test in here.

@svchb svchb requested a review from efaulhaber December 4, 2024 12:55
@@ -96,18 +84,18 @@ end
center = (0.0, 0.0)
NDIMS = 2

# Create a SphereShape (which is a circle in 2D)
# Create a `SphereShape`, which is a circle in 2D
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A circle is the surface, a disk is a filled circle.

Comment on lines 144 to 146
for i in surface_particles
@test isapprox(computed_normals[:, i], expected_normals[:, i], atol=0.05)
end
Copy link
Member

@efaulhaber efaulhaber Dec 9, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@test isapprox(computed_normals[:, surface_particles], expected_normals[:, surface_particles], atol=0.05)

@@ -1,7 +1,7 @@
name = "TrixiParticles"
uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
authors = ["erik.faulhaber <[email protected]>"]
version = "0.2.4-dev"
version = "0.2.4"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it makes more sense to wait for #529 before we release.

Comment on lines +68 to +72
foreach_point_neighbor(system, neighbor_system,
system_coords, neighbor_system_coords,
nhs) do particle, neighbor, pos_diff, distance
colorfield[neighbor] += smoothing_kernel(system, distance)
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
foreach_point_neighbor(system, neighbor_system,
system_coords, neighbor_system_coords,
nhs) do particle, neighbor, pos_diff, distance
colorfield[neighbor] += smoothing_kernel(system, distance)
end
foreach_point_neighbor(neighbor_system, system,
neighbor_system_coords, system_coords,
nhs) do particle, neighbor, pos_diff, distance
colorfield[particle] += smoothing_kernel(system, distance)
end

Otherwise you will get race conditions on multiple threads.

end

@threaded neighbor_system for bnd_particle in eachparticle(neighbor_system)
colorfield[bnd_particle] = colorfield[bnd_particle] / (colorfield[bnd_particle] +
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is colorfield_bnd initialized elsewhere and not here where colorfield is computed?


nparticles = size(u, 2)

# check the threshold has been applied correctly
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# check the threshold has been applied correctly
# Check that the threshold has been applied correctly

#133

Comment on lines +60 to +66
if system.cache.neighbor_count[i] < threshold
@test all(system.cache.surface_normal[:, i] .== 0.0)
else
# For the linear arrangement, surface normals may still be zero
# Adjust the test to account for this possibility
@test true
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you please fix this here then?

Comment on lines +87 to +93
@inline function surface_normal_method(system::FluidSystem)
return system.surface_normal_method
end

@inline function surface_normal_method(system)
return nothing
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Never used.

Comment on lines +411 to +426
function initialize_colorfield!(system, ::BoundaryModelDummyParticles, neighborhood_search)
system_coords = system.coordinates
(; smoothing_kernel, smoothing_length) = system.boundary_model

foreach_point_neighbor(system, system,
system_coords, system_coords,
neighborhood_search,
points=eachparticle(system)) do particle, neighbor, pos_diff,
distance
system.boundary_model.cache.colorfield_bnd[particle] += kernel(smoothing_kernel,
distance,
smoothing_length)
system.boundary_model.cache.neighbor_count[particle] += 1
end
return system
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a test.

return floor(Int, pi * compact_support^2 / particle_spacing^2)
end

@inline @fastpow function ideal_neighbor_count(::Val{3}, particle_spacing, compact_support)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please check the codecov report. This is not covered, but the 2D version is.

# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids"
# Note: This is the simplest form of normal approximation commonly used in SPH and comes
# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries
function calc_normal_akinci!(system, neighbor_system::BoundarySystem, u_system, v,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please check Codecov report. This is missing tests.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants