Merge branch 'wetting_model' into improved_normal_calc_methods

trixi-framework · Dec 4, 2024 · 6860ebc · 6860ebc
2 parents ab9f486 + a05adfb
commit 6860ebc
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 6 deletions.
diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl
@@ -248,6 +248,7 @@ function remove_invalid_normals!(system::FluidSystem,
     for particle in each_moving_particle(system)
 
         # heuristic condition if there is no gas phase to find the free surface
+        # We remove normals for particles which have alot of support e.g. they are in the inside
         if ideal_density_threshold > 0 &&
            ideal_density_threshold * ideal_neighbor_count < cache.neighbor_count[particle]
             if surface_tension.contact_model isa HuberContactModel

diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl
@@ -1,5 +1,6 @@
 function create_fluid_system(coordinates, velocity, mass, density, particle_spacing,
                              surface_tension;
+                             surface_normal_method=ColorfieldSurfaceNormal(),
                              NDIMS=2, smoothing_length=1.0)
     tspan = (0.0, 0.01)
 
@@ -13,7 +14,7 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac
     system = WeaklyCompressibleSPHSystem(fluid, SummationDensity(), state_equation,
                                          SchoenbergCubicSplineKernel{NDIMS}(),
                                          smoothing_length;
-                                         surface_normal_method=ColorfieldSurfaceNormal(),
+                                         surface_normal_method=surface_normal_method,
                                          reference_particle_spacing=particle_spacing,
                                          surface_tension=surface_tension)
 
@@ -25,6 +26,7 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac
     return system, semi, ode
 end
 
+# TODO: change to rect
 function compute_and_test_surface_normals(system, semi, ode; NDIMS=2)
     v0_ode, u0_ode = ode.u0.x
     v = TrixiParticles.wrap_v(v0_ode, system, semi)
@@ -34,6 +36,9 @@ function compute_and_test_surface_normals(system, semi, ode; NDIMS=2)
     TrixiParticles.compute_surface_normal!(system, system.surface_normal_method, v, u,
                                            v0_ode, u0_ode, semi, 0.0)
 
+    TrixiParticles.remove_invalid_normals!(system, SurfaceTensionMorris(),
+                                           system.surface_normal_method)
+
     # After computation, check that surface normals have been computed and are not NaN or Inf
     @test all(isfinite.(system.cache.surface_normal))
     @test all(isfinite.(system.cache.neighbor_count))
@@ -73,7 +78,7 @@ end
 
     system, semi, ode = create_fluid_system(coordinates, velocity, mass, density,
                                             particle_spacing, nothing;
-                                            buffer_size=0, NDIMS=NDIMS)
+                                            NDIMS=NDIMS)
 
     compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS)
 end
@@ -96,8 +101,10 @@ end
     # To get somewhat accurate normals we increase the smoothing length unrealistically
     system, semi, ode = create_fluid_system(coordinates, velocity, mass, density,
                                             particle_spacing, nothing;
-                                            buffer_size=0, NDIMS=NDIMS,
-                                            smoothing_length=3.0 * particle_spacing)
+                                            NDIMS=NDIMS,
+                                            smoothing_length=3.0 * particle_spacing,
+                                            surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1,
+                                                                                          ideal_density_threshold=0.9))
 
     compute_and_test_surface_normals(system, semi, ode; NDIMS=NDIMS)
 
@@ -130,8 +137,10 @@ end
         end
     end
 
-    # Compare computed normals to expected normals for surface particles
-    @test isapprox(computed_normals, expected_normals, atol=0.05)
+    for i in surface_particles
+        @test isapprox(computed_normals[:, i][:, i], expected_normals[:, i][:, i],
+                       atol=0.04)
+    end
 
     # Optionally, check that normals for interior particles are zero
     # for i in setdiff(1:nparticles, surface_particles)