Add wetted-area contact angle

NEWS.md

@@ -4,11 +4,39 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
-## Version 0.5.3
-
-### Features
-
-- Added the computation of boundary normals for `RectangularTank`s and `SphereShape`s.
+## Version 0.5.3
+
+### API Changes
+
+- Corrected `SurfaceTensionMorris` to apply its local CSF acceleration once per particle and
+  retain the required one-phase surface delta. Previous coefficients compensated implicitly
+  for a dimensionally incomplete force repeated once per fluid neighbor and must be recalibrated.
+- For `SurfaceTensionMorris` and `SurfaceTensionMomentumMorris` with
+  `ColorfieldSurfaceNormal`, `ideal_density_threshold` now denotes a fraction of the continuous
+  complete-support kernel moment instead of an integer neighbor-count fraction. The default zero
+  still disables interior filtering.
+
+### Features
+
+- Reworked `SurfaceTensionMomentumMorris` as a balanced continuum-surface-stress operator with a
+  symmetric support correction. It conserves pairwise linear momentum without a cached stress
+  tensor or global reduction.
+- Added opt-in activity-weighted Shepard normal smoothing for Morris CSF/CSS and expanded VTK
+  output with raw and capillary normals, surface delta, activity, support, force, and reconstructed
+  stress diagnostics.
+- Added `FreeSurfaceTangentialShifting`, an opt-in treatment that smoothly removes the
+  interface-normal component of Sun particle shifting while retaining full interior shifting.
+- Added `InterfaceAwareTensileInstabilityControl`, an opt-in pressure formulation that applies
+  tensile-instability control in fluid interiors and smoothly disables it at free surfaces.
+- Added optional per-particle `surface_measure` quadrature data to dummy-particle boundaries.
+- Added `WettedAreaContactAngle`, an opt-in Young wall-energy model for supported 3D Morris CSS
+  simulations with force- and torque-conserving fixed-wall and rigid-body reactions.
+- Added C1 interface activation for Morris CSF and CSS with `ColorfieldSurfaceNormal`.
+  Color-gradient and continuous support-moment indicators taper the physical surface delta without
+  another neighbor pass.
+- Added `CorrectedCSFSurfaceNormal`, an explicit free-surface implementation of the C-CSF
+  interface geometry from Vergnaud et al. (2022) for `SurfaceTensionMorris`.
+- Added the computation of boundary normals for `RectangularTank`s and `SphereShape`s.
 - Added `flush` keyword argument to `InfoCallback` to flush `stdout` after each output,
   useful for monitoring progress in real-time on clusters or batch systems (#1246).
 - Added the number of split integration time steps to the `InfoCallback` output

README.md

@@ -45,24 +45,24 @@ It offers intuitive configuration, robust pre- and post-processing, and vendor-a
 We provide several example simulation setups in the `examples` folder (which can be accessed from Julia via `examples_dir()`).
 
 <table align="center" border="0">
-  <tr>
-  </tr>
   <tr>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/683e9363-5705-49cc-9a5c-3b47d73ea4b8" style="width: 80% !important;"/><br><figcaption>2D Dam Break</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/683e9363-5705-49cc-9a5c-3b47d73ea4b8" style="width: 80% !important;"/><br />
+      <figcaption>2D Dam Break</figcaption>
     </td>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/c10faddf-0400-47c9-b225-f5d286a8ecb8" style="width: 80% !important;"/><br><figcaption>Moving Wall</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/c10faddf-0400-47c9-b225-f5d286a8ecb8" style="width: 80% !important;"/><br />
+      <figcaption>Moving Wall</figcaption>
     </td>
   </tr>
-  <tr>
-  </tr>
   <tr>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/e05ace63-e330-441a-a391-eda3d2764074" style="width: 80% !important;"/><br><figcaption>Oscillating Beam</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/e05ace63-e330-441a-a391-eda3d2764074" style="width: 80% !important;"/><br />
+      <figcaption>Oscillating Beam</figcaption>
     </td>
     <td align="center">
-      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/ada0d554-e0ba-44ed-923d-2b77ef252258" style="width: 80% !important;"/><br><figcaption>Dam Break with Elastic Plate</figcaption>
+      <img src="https://github.com/trixi-framework/TrixiParticles.jl/assets/10238714/ada0d554-e0ba-44ed-923d-2b77ef252258" style="width: 80% !important;"/><br />
+      <figcaption>Dam Break with Elastic Plate</figcaption>
     </td>
   </tr>
 </table>

docs/literate/src/tut_2d_geometry.jl

@@ -0,0 +1,155 @@
+# # [Setting up a 2D simulation from geometry files](@id tut_2d_geometry)
+
+# In this tutorial, we build two genuine 2D setups from geometry files:
+# 1. a curved pipe, where one geometry file defines the outer wall envelope and a second
+#    one defines the empty channel cut out of it,
+# 2. a dam-break basin with a coastline profile, where one geometry file defines the
+#    filled coastline wall together with the seawall on the right.
+#
+# For a real 2D setup, we use 2D geometry formats such as `.asc` or `.dxf`.
+# STL files are surface meshes and therefore naturally lead to thin 3D setups instead.
+
+# First, we import TrixiParticles.jl together with
+# `OrdinaryDiffEqLowStorageRK` of
+# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
+# and [Plots.jl](https://docs.juliaplots.org/stable/).
+using TrixiParticles
+using OrdinaryDiffEqLowStorageRK
+using Plots
+
+# ## Resolution
+
+# We use the same particle spacing for the fluid and for the wall geometries.
+particle_spacing = 0.03
+fluid_density = 1000.0
+gravity = 9.81
+sound_speed = 10.0
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7)
+nothing # hide
+
+# ## Loading 2D geometry files
+
+# The following helper loads a closed 2D geometry file and samples particles in its interior:
+# 1. load the polygon with [`load_geometry`](@ref),
+# 2. fill the polygon with [`ComplexShape`](@ref).
+#
+# This creates a true 2D solid region instead of a hollow shell around the polygon edges.
+function solid_from_geometry_file(file; particle_spacing, density)
+    geometry = load_geometry(file)
+    solid = ComplexShape(geometry; particle_spacing, density,
+                         grid_offset=0.5particle_spacing)
+
+    return (; geometry, solid)
+end
+
+# ## A curved pipe from two filled geometries
+
+# The pipe wall is a solid L-shaped region with a channel cut out of it:
+# 1. one geometry file describes the outer pipe envelope,
+# 2. one geometry file describes the empty channel,
+# 3. the `setdiff` operation subtracts the channel from the solid envelope.
+pipe_outer_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data",
+                         "curved_pipe_outer_2d.asc")
+pipe_channel_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data",
+                           "curved_pipe_channel_2d.asc")
+
+pipe_outer = solid_from_geometry_file(pipe_outer_file; particle_spacing,
+                                      density=fluid_density)
+pipe_channel = load_geometry(pipe_channel_file)
+
+pipe_setup = (; wall=setdiff(pipe_outer.solid, pipe_channel),
+              outer_geometry=pipe_outer.geometry,
+              channel_geometry=pipe_channel)
+
+# ## A dam-break basin with a coastline profile
+
+# In the second setup, a single 2D geometry file defines a filled coastline wall:
+# the beach profile on top, a finite wall thickness below it, and the seawall on the right.
+coast_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data",
+                    "coastline_profile_2d.asc")
+coast = solid_from_geometry_file(coast_file; particle_spacing, density=fluid_density)
+
+# The geometry file gives the coastline bed and the right wall as a solid region.
+# We add the left wall explicitly as a rectangular particle block and place a
+# 1.5x taller rectangular dam-break water column next to it.
+left_wall = RectangularShape(particle_spacing, (5, 50), (0.0, -0.12),
+                             density=fluid_density)
+reservoir = RectangularShape(particle_spacing, (28, 42), (0.15, 0.03),
+                             acceleration=(0.0, -gravity),
+                             state_equation=state_equation)
+coast_setup = (; geometry=coast.geometry,
+               wall=union(coast.solid, left_wall),
+               fluid=setdiff(reservoir, coast.geometry))
+
+p_pipe = plot(pipe_setup.wall, label="wall", title="Curved pipe",
+              markerstrokewidth=0, markersize=4)
+plot!(p_pipe, showaxis=false, aspect_ratio=:equal,
+      xlims=(-0.03, 1.23), ylims=(-0.03, 1.23))
+
+p_coast = plot(coast_setup.fluid, coast_setup.wall,
+               labels=["fluid" "wall"], title="Coastline dam break",
+               markerstrokewidth=0, markersize=3)
+plot!(p_coast, showaxis=false, aspect_ratio=:equal,
+      xlims=(0.0, 2.75), ylims=(-0.15, 1.35))
+
+plot(p_pipe, p_coast, layout=(1, 2), size=(900, 360))
+savefig("tut_2d_geometry_plot.png"); # hide
+# ![2D geometry based initial conditions](tut_2d_geometry_plot.png)
+
+# ## Building the simulation systems
+
+# To keep the example focused, we continue with the coastline setup.
+# From this point on, the simulation setup is the same as in other 2D simulation files.
+setup = coast_setup
+tspan = (0.0, 0.03)
+nothing # hide
+
+# We define the state equation, smoothing kernel, and viscosity for a
+# weakly compressible SPH simulation.
+smoothing_length = 1.2 * particle_spacing
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
+
+fluid_density_calculator = ContinuityDensity()
+density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
+
+fluid_system = WeaklyCompressibleSPHSystem(setup.fluid;
+                                           density_calculator=fluid_density_calculator,
+                                           state_equation, smoothing_kernel,
+                                           smoothing_length, viscosity=viscosity,
+                                           density_diffusion=density_diffusion,
+                                           acceleration=(0.0, -gravity))
+nothing # hide
+
+# For the wall, we reuse the combined solid wall particles created above. The high-level
+# constructor obtains the smoothing kernel, smoothing length, and state equation from the fluid.
+boundary_model = BoundaryModelDummyParticles(setup.wall; fluid_system)
+boundary_system = WallBoundarySystem(setup.wall, boundary_model)
+nothing # hide
+
+# ## Semidiscretization
+
+# With fluid and wall particles defined, we can build the
+# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization) exactly as in other tutorials.
+semi = Semidiscretization(fluid_system, boundary_system)
+ode = semidiscretize(semi, tspan)
+nothing # hide
+
+# ## Time integration
+
+# The setup is now complete.
+# To start the simulation, run for example
+# ```julia
+# callbacks = CallbackSet(InfoCallback(interval=10))
+# sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks)
+# ```
+# This is the same final step as in [the basic setup tutorial](@ref tut_setup).
+callbacks = CallbackSet(InfoCallback(interval=10))
+nothing # hide
+
+sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks) #!md
+
+# For more accurate body-fitted particles around sharper features, you can also
+# apply the [particle packing workflow](@ref tut_packing) to the 2D geometry files
+# before starting the simulation.

docs/literate/src/tut_packing.jl

@@ -75,8 +75,8 @@ plot!(right_margin=5Plots.mm) #hide
 # ## Creating an initial configuration of boundary particles
 
 # To create the initial configuration of the boundary particles,
-# we use the sampled points of the SDF whose signed distance lies between 0
-# and `boundary_thickness`.
+# we use the sampled points of the SDF whose signed distance lies between the
+# geometry offset implied by `place_on_shell` and `boundary_thickness`.
 # Here, we need to specify the `density` of the boundary particles.
 # As an example, we choose `1.0` for all particles.
 # This gives us an [`InitialCondition`](@ref InitialCondition) for the boundary particles.
@@ -125,7 +125,8 @@ plot!(geometry, linestyle=:dash, label=nothing, showaxis=false, color=:black,
 # ## Particle packing
 
 # In the following, we will essentially follow the same steps described in the fluid tutorials.
-# That means we will generate systems that are then passed to the [`Semidiscretization`](@ref).
+# That means we will generate systems that are then passed to the
+# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization).
 # The difference from a typical physical simulation is that we use [`ParticlePackingSystem`](@ref),
 # which does not represent any physical law. Instead, we only use the simulation framework to time-integrate
 # the packing process.
@@ -211,7 +212,7 @@ plot!(geometry, seriestype=:path, color=:black, label=nothing, linewidth=2)
 boundary_system = ParticlePackingSystem(boundary_sampled; is_boundary=true,
                                         smoothing_kernel, smoothing_length,
                                         boundary_compress_factor=0.7, signed_distance_field,
-                                        background_pressure)
+                                        boundary_thickness, background_pressure)
 
 # We can now couple the boundary system with the interior system:
 semi = Semidiscretization(packing_system, boundary_system)
@@ -251,10 +252,12 @@ fixed_system = ParticlePackingSystem(packed_ic; smoothing_kernel, smoothing_leng
 
 # Now we define a rectangular domain that we want to pack.
 # In practice, you could create any `InitialCondition` that encloses your complex geometry.
-tank_domain = RectangularTank(particle_spacing, (4, 4), (0, 0), min_coordinates=(-1, -2),
-                              density)
+domain_size = (4, 4)
+n_particles_per_dimension = round.(Int, domain_size ./ particle_spacing)
+tank_domain = RectangularShape(particle_spacing, n_particles_per_dimension, (-1, -2);
+                               density)
 
-sampled_outer_domain = setdiff(tank_domain.fluid, packed_ic)
+sampled_outer_domain = setdiff(tank_domain, packed_ic)
 
 # If we plot these two `InitialCondition`s, we can see
 # that the geometry interface is not properly represented yet.

docs/literate/src/tut_rigid_body_fsi.jl

@@ -165,12 +165,8 @@ nothing # hide
 # See [the docs on dummy particles](@ref boundary_models) for a definition for these terms.
 
 boundary_density_calculator = AdamiPressureExtrapolation()
-tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
-                                                  tank.boundary.mass,
-                                                  boundary_density_calculator,
-                                                  fluid_smoothing_kernel,
-                                                  fluid_smoothing_length;
-                                                  state_equation)
+tank_boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system,
+                                                  boundary_density_calculator)
 
 boundary_system = WallBoundarySystem(tank.boundary, tank_boundary_model)
 nothing # hide

docs/literate/src/tut_setup.jl

@@ -131,20 +131,18 @@ nothing # hide
 
 # To model the boundary, we use particle-based boundary conditions, in which particles
 # are sampled in the boundary that interact with the fluid particles to avoid penetration.
-# In order to define a boundary system, we first have to choose a boundary model,
-# which defines how the fluid interacts with boundary particles.
-# We will use the [`BoundaryModelDummyParticles`](@ref) with [`AdamiPressureExtrapolation`](@ref).
-# See [here](@ref boundary_models) for a comprehensive overview over boundary models.
-boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
-                                             AdamiPressureExtrapolation(),
-                                             smoothing_kernel, smoothing_length;
-                                             state_equation)
+# Here, we explicitly choose the dummy-particle boundary model and use its high-level
+# builder to infer kernel and equation-of-state-related settings from the adjacent
+# fluid system. See [here](@ref boundary_models) for a comprehensive overview over
+# boundary models.
+boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system)
 boundary_system = WallBoundarySystem(tank.boundary, boundary_model)
 nothing # hide
 
 # ## [Semidiscretization](@id tut_setup_semi)
 
-# The key component of every simulation is the [`Semidiscretization`](@ref),
+# The key component of every simulation is the
+# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization),
 # which couples all systems of the simulation.
 # All simulation methods in TrixiParticles.jl are semidiscretizations, which discretize
 # the equations in space to provide an ordinary differential equation that still

docs/make.jl

@@ -42,6 +42,8 @@ Literate.markdown(joinpath("docs", "literate", "src", "tut_rigid_body_fsi.jl"),
                   joinpath("docs", "src", "tutorials"))
 Literate.markdown(joinpath("docs", "literate", "src", "tut_packing.jl"),
                   joinpath("docs", "src", "tutorials"))
+Literate.markdown(joinpath("docs", "literate", "src", "tut_2d_geometry.jl"),
+                  joinpath("docs", "src", "tutorials"))
 
 copy_file("AUTHORS.md",
           "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")
@@ -94,6 +96,8 @@ makedocs(sitename="TrixiParticles.jl",
                                                                                  "tut_rigid_body_fsi.md")
                  ],
                  "Preprocessing" => [
+                     "Setting up a 2D simulation from geometry files" => joinpath("tutorials",
+                                                                                  "tut_2d_geometry.md"),
                      "Particle packing tutorial" => joinpath("tutorials",
                                                              "tut_packing.md")
                  ]