Generalize surface normal calc #539

svchb · 2024-05-29T10:49:04Z

Depends on #599
extracted from #584

Conflicts: src/schemes/fluid/weakly_compressible_sph/system.jl

…vchb/TrixiParticles.jlOpen into generalize_surface_normal_calc

Conflicts: examples/fluid/dam_break_oil_film_2d.jl src/schemes/fluid/weakly_compressible_sph/state_equations.jl

NEWS.md

src/schemes/boundary/dummy_particles/dummy_particles.jl

src/schemes/fluid/surface_normal_sph.jl

src/schemes/fluid/entropically_damped_sph/system.jl

…vchb/TrixiParticles.jlOpen into generalize_surface_normal_calc

docs/src/systems/fluid.md

efaulhaber · 2024-12-03T08:06:35Z

src/schemes/fluid/weakly_compressible_sph/system.jl

    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


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

Yes it will be used in #584

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

Great unnecessary conflicts in 5 Branches...

Rebasing should really eliminate all your problems.

test/schemes/fluid/surface_normal_sph.jl

efaulhaber · 2024-12-03T08:50:38Z

test/schemes/fluid/surface_normal_sph.jl

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


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

The inside will stay at 0.0. though.

test/schemes/fluid/surface_normal_sph.jl

efaulhaber · 2024-12-03T09:09:38Z

test/schemes/fluid/surface_normal_sph.jl

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


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.

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

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

That also doesn't work

Expression: isapprox(computed_normals[:, surface_particles], expected_normals[:, surface_particles], atol = 0.04) Evaluated: isapprox([0.3480639311878508 -0.3480639311878508 … 0.34806393118785095 -0.3480639311878508; 0.9374707994418061 0.9374707994418061 … -0.937470799441806 -0.937470799441806], [0.31622776601683794 -0.31622776601683794 … 0.31622776601683794 -0.31622776601683794; 0.9486832980505138 0.9486832980505138 … -0.9486832980505138 -0.9486832980505138]; atol = 0.04)

According to the docs of isapprox:

The norm keyword defaults to abs for numeric (x,y) and to LinearAlgebra.norm for arrays (where an alternative norm choice is sometimes useful).

In the loop, you checked the Linf norm of the difference, here you check the L2 norm, which becomes larger with the size of the vector when the Linf norm stays the same:

julia> norm(ones(100_000) - zeros(100_000)) 316.22776601683796

For isapprox, this means:

julia> isapprox(ones(100_000), zeros(100_000), atol=300) false julia> isapprox(ones(100_000), zeros(100_000), atol=400) true julia> isapprox(ones(100_000), zeros(100_000), atol=1.0, norm=x -> norm(x, Inf)) true julia> isapprox(ones(100_000), zeros(100_000), atol=0.9, norm=x -> norm(x, Inf)) false

efaulhaber · 2024-12-03T09:09:43Z

test/schemes/fluid/surface_normal_sph.jl

+    # 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


This is also part of #584

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

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.

The problem is that not only the outer layer has the surface normal set but this goes some layers deep so the only possiblity is to eliminate several of the outer layers.

Than the inside is 0

What happens between the surface and the inner layers where it's 0? Is this just a shorter normal then? Can you test that the normal either points in the correct direction or is almost zero?

…vchb/TrixiParticles.jlOpen into generalize_surface_normal_calc

test/schemes/fluid/surface_normal_sph.jl

efaulhaber · 2024-12-09T07:59:00Z

test/schemes/fluid/surface_normal_sph.jl

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


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

efaulhaber · 2024-12-09T08:01:37Z

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"]
-version = "0.2.4-dev"
+version = "0.2.4"


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

Depends on how long that still will take...

efaulhaber · 2024-12-09T08:13:31Z

src/schemes/fluid/surface_normal_sph.jl

+    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


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.

Doesn't work there will be no neighbors

Selfnote: I still need to revert this since this has been only changed at the moment from #666 and up

foreach_point_neighbor(neighbor_system, system, ..., particles=eachparticle(neighbor_system)), or the loop will go over all moving particles, which are none for a boundary system.

src/schemes/fluid/surface_normal_sph.jl

test/schemes/fluid/surface_normal_sph.jl

efaulhaber · 2024-12-09T08:16:45Z

test/schemes/fluid/surface_normal_sph.jl

+        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


Can you please fix this here then?

efaulhaber · 2024-12-09T08:18:21Z

src/schemes/fluid/fluid.jl

+@inline function surface_normal_method(system::FluidSystem)
+    return system.surface_normal_method
+end
+
+@inline function surface_normal_method(system)
+    return nothing
+end


Never used.

This is used in #584

Then add it there.

efaulhaber · 2024-12-09T08:19:03Z

src/schemes/boundary/system.jl

+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


Please add a test.

I have added a test but this test will change in all subsequent PRs

I think that is better than reducing coverage here and potentially forgetting about it.

efaulhaber · 2024-12-09T08:21:59Z

src/general/system.jl

+    return floor(Int, pi * compact_support^2 / particle_spacing^2)
+end
+
+@inline @fastpow function ideal_neighbor_count(::Val{3}, particle_spacing, compact_support)


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

A test for this exists in #584

Can you add it here or does it require other code from #584?

efaulhaber · 2024-12-09T08:24:06Z

src/schemes/fluid/surface_normal_sph.jl

+# 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,


Please check Codecov report. This is missing tests.

Testing this doesn't make much sense. Since there is no analytical value. And it will change in all subsequent PRs.

Let me restate that there is a analytical value but the results will not be close. Just closer than without it.

There must be some useful way to test this. Maybe with a high resolution and very large smoothing length?

svchb added 30 commits September 6, 2023 01:01

Merge remote-tracking branch 'origin/main' into main

f90f017

Merge branch 'trixi-framework:main' into main

372d67c

Merge branch 'main' of github.com:svchb/TrixiParticles.jl

f564478

Merge branch 'trixi-framework:main' into main

9c9e182

Merge branch 'main' of github.com:svchb/TrixiParticles.jl

90897f6

Merge branch 'main' of github.com:svchb/TrixiParticles.jl

1df6912

Merge branch 'main' of github.com:svchb/TrixiParticles.jl

6449b5f

merge by hand

439210f

Merge branch 'HEAD' into ipf_surface_tension

b31640f

Merge branch 'HEAD' into ipf_surface_tension

675c030

Merge branch 'HEAD' into ipf_surface_tension

2b0b09c

Merge branch 'HEAD' into ipf_surface_tension

b812869

Merge branch 'HEAD' into ipf_surface_tension

c3884d7

Merge branch 'HEAD' into ipf_surface_tension

45938eb

Merge tag 'HEAD' into ipf_surface_tension

3bdfc20

Merge branch 'conflict1'

de4f0e8

Conflicts: src/schemes/fluid/weakly_compressible_sph/system.jl

fix merge problems

c49dd87

fix merge

780cee3

fix merge

7550c3d

fix merge

87c21b9

fix merge

9cab05f

fix merge

896faef

more merge fixes

b567c55

format

3de9a75

fix

06a1223

fix

31d4404

fix

a95a54a

fix

ce74643

fix

17da7ef

more fixes

f002d2b

svchb added 8 commits November 20, 2024 16:11

Merge branch 'generalize_surface_normal_calc' of https://github.com/s…

01bddb0

…vchb/TrixiParticles.jlOpen into generalize_surface_normal_calc

Update NEWS.md

dc6fa70

Merge branch 'main'

ffc83d3

Conflicts: examples/fluid/dam_break_oil_film_2d.jl src/schemes/fluid/weakly_compressible_sph/state_equations.jl

Merge branch 'main' into generalize_surface_normal_calc

e424cfa

fix

041861e

fix

21f08b2

fix

0b63bca

fix

b7f448c

efaulhaber reviewed Nov 25, 2024

View reviewed changes

NEWS.md Outdated Show resolved Hide resolved

src/schemes/boundary/dummy_particles/dummy_particles.jl Outdated Show resolved Hide resolved

efaulhaber reviewed Nov 25, 2024

View reviewed changes

src/schemes/boundary/dummy_particles/dummy_particles.jl Outdated Show resolved Hide resolved

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

src/schemes/fluid/entropically_damped_sph/system.jl Show resolved Hide resolved

efaulhaber mentioned this pull request Nov 25, 2024

Move fields that are only sometimes needed to the cache? #676

Open

svchb added 3 commits December 2, 2024 12:35

Merge branch 'main' into generalize_surface_normal_calc

f4fbbb3

review fixes

589d416

Merge branch 'generalize_surface_normal_calc' of https://github.com/s…

1263e1a

…vchb/TrixiParticles.jlOpen into generalize_surface_normal_calc

svchb requested a review from efaulhaber December 2, 2024 11:36

efaulhaber requested changes Dec 3, 2024

View reviewed changes

svchb requested a review from efaulhaber December 4, 2024 12:55

svchb added 4 commits December 4, 2024 13:56

review comments

dd56f47

Merge branch 'generalize_surface_normal_calc' of https://github.com/s…

306f96d

…vchb/TrixiParticles.jlOpen into generalize_surface_normal_calc

fix test

93630aa

typo

9f758be

efaulhaber reviewed Dec 9, 2024

View reviewed changes

efaulhaber requested changes Dec 9, 2024

View reviewed changes

svchb added 7 commits December 12, 2024 16:15

Merge branch 'main' into generalize_surface_normal_calc

cd26885

review updates

80e51ba

review update

4721aec

add boundary system

98a5f6c

update

e05edea

update

36741d7

format

ba8765f

Generalize surface normal calc #539

Are you sure you want to change the base?

Generalize surface normal calc #539

Conversation

svchb commented May 29, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

efaulhaber Dec 9, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

efaulhaber Dec 9, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

svchb commented May 29, 2024 •

edited

Loading

efaulhaber Dec 9, 2024 •

edited

Loading

efaulhaber Dec 9, 2024 •

edited

Loading