Generalize surface normal calc

NEWS.md

@@ -3,51 +3,85 @@
 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.2.4
+
+### Features
+
+- Support for surface tension was added to EDAC (#539)
+- A method to prevent penetration of fast moving particles with solids was added (#498)
+- Added the callback `SteadyStateReachedCallback` to detect converging simulations (#601)
+- Added Ideal Gas State Equation (#607)
+
+### Documentation
+
+- GPU Support Documentation was added (#660)
+- A new user tutorial was added (#514)
+
+### Fixes
+
+- Diverse Doc fixes (#663, #659, #637, #658, #664)
+- Simulations can be run with `Float32` (#662)
+
+### Refactored
+
+- Surface normal calculation was moved from surface_tension.jl to surface_normal_sph.jl (#539)
+
 ## Version 0.2.3
 
 ### Highlights
+
 Transport Velocity Formulation (TVF) based on the work of Ramachandran et al. "Entropically damped artificial compressibility for SPH" (2019) was added.
 
 ## Version 0.2.2
 
 ### Highlights
+
 Hotfix for threaded sampling of complex geometries.
 
 ## Version 0.2.1
 
 ### Highlights
+
 Particle sampling of complex geometries from `.stl` and `.asc` files.
 
 ## Version 0.2.0
 
 ### Removed
+
 Use of the internal neighborhood search has been removed and replaced with PointNeighbors.jl.
 
 ## Development Cycle 0.1
 
 ### Highlights
 
 #### Discrete Element Method
+
 A basic implementation of the discrete element method was added.
 
 #### Surface Tension and Adhesion Model
+
 A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH.
 
 #### Support for Open Boundaries
+
 Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC.
 
 ## Pre Initial Release (v0.1.0)
+
 This section summarizes the initial features that TrixiParticles.jl was released with.
 
 ### Highlights
 #### EDAC
+
 An implementation of EDAC (Entropically Damped Artificial Compressibility) was added,
 which allows for more stable simulations compared to basic WCSPH and reduces spurious pressure oscillations.
 
 #### WCSPH
+
 An implementation of WCSPH (Weakly Compressible Smoothed Particle Hydrodynamics), which is the classical SPH approach.
 
 Features:
+
 - Correction schemes (Shepard (0. Order) ... MixedKernelGradient (1. Order))
 - Density reinitialization
 - Kernel summation and Continuity equation density formulations
@@ -57,4 +91,5 @@ Features:
 
 
 #### TLSPH
+
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).

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"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/make.jl

@@ -139,8 +139,8 @@ makedocs(sitename="TrixiParticles.jl",
                      "Util" => joinpath("general", "util.md")
                  ],
                  "Systems" => [
-                     "Discrete Element Method (Solid)" => joinpath("systems",
-                                                                   "dem.md"),
+                     "Fluid Models" => joinpath("systems", "fluid.md"),
+                     "Discrete Element Method (Solid)" => joinpath("systems", "dem.md"),
                      "Weakly Compressible SPH (Fluid)" => joinpath("systems",
                                                                    "weakly_compressible_sph.md"),
                      "Entropically Damped Artificial Compressibility for SPH (Fluid)" => joinpath("systems",

docs/src/systems/fluid.md

@@ -0,0 +1,71 @@
+# [Fluid Models](@id fluid_models)
+In the following are common models for fluid effects that are used by both EDAC and WCSPH.
+
+
+## [Surface Normals](@id surface_normal)
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "surface_normal_sph.jl")]
+```
+
+## [Surface Tension](@id surface_tension)
+
+### Akinci-based intra-particle force surface tension and wall adhesion model
+The work by Akinci proposes three forces:
+- a cohesion force
+- a surface area minimization force
+- a wall adhesion force
+
+The classical model is composed of the curvature minimization and cohesion force.
+
+#### Cohesion force
+The model calculates the cohesion force based on the distance between particles and the support radius ``h_c``.
+This force is determined using two distinct regimes within the support radius:
+- For particles closer than half the support radius,
+  a repulsive force is calculated to prevent particle clustering too tightly,
+  enhancing the simulation's stability and realism.
+- Beyond half the support radius and within the full support radius,
+  an attractive force is computed, simulating the effects of surface tension that draw particles together.
+The cohesion force, ``F_{\text{cohesion}}``, for a pair of particles is given by:
+```math
+F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert},
+```
+where:
+- ``\sigma`` represents the surface tension coefficient, adjusting the overall strength of the cohesion effect.
+- ``C`` is a scalar function of the distance between particles.
+
+The cohesion kernel ``C`` is defined as
+```math
+C(r)=\frac{32}{\pi h_c^9}
+\begin{cases}
+(h_c-r)^3 r^3, & \text{if } 2r > h_c \\
+2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c \\
+0, & \text{otherwise}
+\end{cases}
+```
+
+#### Surface area minimization force
+To model the minimization of the surface area and curvature of the fluid, a curvature force is used, which is calculated as
+```math
+F_{\text{curvature}} = -\sigma (n_a - n_b)
+```
+
+#### Wall adhesion force
+The wall adhesion model proposed by Akinci et al. is based on a kernel function which is 0 from 0.0 to 0.5 support radiia with a maximum at 0.75.
+With the force calculated with an adhesion coefficient ``\beta`` as
+```math
+F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert},
+```
+with ``A`` being the adhesion kernel defined as
+```math
+A(r)= \frac{0.007}{h_c^{3.25}}
+\begin{cases}
+\sqrt[4]{- \frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c \\
+0, & \text{otherwise.}
+\end{cases}
+```
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "surface_tension.jl")]
+```

docs/src/systems/weakly_compressible_sph.md

@@ -126,64 +126,3 @@ Modules = [TrixiParticles]
 Pages = [joinpath("general", "corrections.jl")]
 ```
 
-## [Surface Tension](@id surface_tension)
-
-### Akinci-based intra-particle force surface tension and wall adhesion model
-The work by Akinci proposes three forces:
-- a cohesion force
-- a surface area minimization force
-- a wall adhesion force
-
-The classical model is composed of the curvature minimization and cohesion force.
-
-#### Cohesion force
-The model calculates the cohesion force based on the distance between particles and the support radius ``h_c``.
-This force is determined using two distinct regimes within the support radius:
-- For particles closer than half the support radius,
-  a repulsive force is calculated to prevent particle clustering too tightly,
-  enhancing the simulation's stability and realism.
-- Beyond half the support radius and within the full support radius,
-  an attractive force is computed, simulating the effects of surface tension that draw particles together.
-The cohesion force, ``F_{\text{cohesion}}``, for a pair of particles is given by:
-```math
-F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert},
-```
-where:
-- ``\sigma`` represents the surface tension coefficient, adjusting the overall strength of the cohesion effect.
-- ``C`` is a scalar function of the distance between particles.
-
-The cohesion kernel ``C`` is defined as
-```math
-C(r)=\frac{32}{\pi h_c^9}
-\begin{cases}
-(h_c-r)^3 r^3, & \text{if } 2r > h_c \\
-2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c \\
-0, & \text{otherwise}
-\end{cases}
-```
-
-#### Surface area minimization force
-To model the minimization of the surface area and curvature of the fluid, a curvature force is used, which is calculated as
-```math
-F_{\text{curvature}} = -\sigma (n_a - n_b)
-```
-
-#### Wall adhesion force
-The wall adhesion model proposed by Akinci et al. is based on a kernel function which is 0 from 0.0 to 0.5 support radiia with a maximum at 0.75.
-With the force calculated with an adhesion coefficient ``\beta`` as
-```math
-F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert},
-```
-with ``A`` being the adhesion kernel defined as
-```math
-A(r)= \frac{0.007}{h_c^{3.25}}
-\begin{cases}
-\sqrt[4]{- \frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c \\
-0, & \text{otherwise.}
-\end{cases}
-```
-
-```@autodocs
-Modules = [TrixiParticles]
-Pages = [joinpath("schemes", "fluid", "surface_tension.jl")]
-```

examples/fluid/dam_break_2d.jl

@@ -61,7 +61,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            smoothing_length, viscosity=viscosity,
                                            density_diffusion=density_diffusion,
                                            acceleration=(0.0, -gravity), correction=nothing,
-                                           surface_tension=nothing)
+                                           surface_tension=nothing,
+                                           reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Boundary

examples/fluid/dam_break_oil_film_2d.jl

@@ -30,13 +30,16 @@ oil_viscosity = ViscosityMorris(nu=nu_sim_oil)
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               sol=nothing, fluid_particle_spacing=fluid_particle_spacing,
               viscosity=ViscosityMorris(nu=nu_sim_water), smoothing_length=smoothing_length,
-              gravity=gravity, density_diffusion=nothing, sound_speed=sound_speed)
+              gravity=gravity, tspan=tspan, density_diffusion=nothing,
+              sound_speed=sound_speed, prefix="")
 
 # ==========================================================================================
 # ==== Setup oil layer
 
 oil_size = (W, 0.1 * H)
 oil_density = 700.0
+oil_eos = StateEquationCole(; sound_speed, reference_density=oil_density, exponent=1,
+                            clip_negative_pressure=false)
 
 oil = RectangularShape(fluid_particle_spacing,
                        round.(Int, oil_size ./ fluid_particle_spacing),
@@ -50,12 +53,12 @@ end
 oil_state_equation = StateEquationCole(; sound_speed, reference_density=oil_density,
                                        exponent=1, clip_negative_pressure=false)
 oil_system = WeaklyCompressibleSPHSystem(oil, fluid_density_calculator,
-                                         oil_state_equation, smoothing_kernel,
+                                         oil_eos, smoothing_kernel,
                                          smoothing_length, viscosity=oil_viscosity,
-                                         density_diffusion=density_diffusion,
                                          acceleration=(0.0, -gravity),
-                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.02),
-                                         correction=AkinciFreeSurfaceCorrection(oil_density))
+                                         surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.01),
+                                         correction=AkinciFreeSurfaceCorrection(oil_density),
+                                         reference_particle_spacing=fluid_particle_spacing)
 
 # ==========================================================================================
 # ==== Simulation

examples/fluid/falling_water_spheres_2d.jl

@@ -40,7 +40,7 @@ sphere2 = SphereShape(fluid_particle_spacing, sphere_radius, sphere2_center,
 
 # ==========================================================================================
 # ==== Fluid
-fluid_smoothing_length = 1.0 * fluid_particle_spacing
+fluid_smoothing_length = 1.0 * fluid_particle_spacing - eps()
 fluid_smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
@@ -50,13 +50,13 @@ alpha = 8 * nu / (fluid_smoothing_length * sound_speed)
 viscosity = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
 density_diffusion = DensityDiffusionAntuono(sphere2, delta=0.1)
 
-sphere_surface_tension = WeaklyCompressibleSPHSystem(sphere1, fluid_density_calculator,
-                                                     state_equation, fluid_smoothing_kernel,
+sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1, fluid_smoothing_kernel,
                                                      fluid_smoothing_length,
-                                                     viscosity=viscosity,
+                                                     sound_speed, viscosity=viscosity,
+                                                     density_calculator=ContinuityDensity(),
                                                      acceleration=(0.0, -gravity),
                                                      surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05),
-                                                     correction=AkinciFreeSurfaceCorrection(fluid_density))
+                                                     reference_particle_spacing=fluid_particle_spacing)
 
 sphere = WeaklyCompressibleSPHSystem(sphere2, fluid_density_calculator,
                                      state_equation, fluid_smoothing_kernel,
@@ -75,16 +75,16 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar
                                              viscosity=ViscosityAdami(nu=wall_viscosity))
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
-                                    adhesion_coefficient=0.001)
+                                    adhesion_coefficient=1.0)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(boundary_system, sphere_surface_tension, sphere)
+semi = Semidiscretization(sphere_surface_tension, sphere, boundary_system)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)
-saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out", prefix="",
-                                         write_meta_data=true)
+saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out",
+                                         prefix="", write_meta_data=true)
 
 callbacks = CallbackSet(info_callback, saving_callback)
 

examples/fluid/falling_water_spheres_3d.jl

@@ -3,35 +3,35 @@ using OrdinaryDiffEq
 
 # ==========================================================================================
 # ==== Resolution
-fluid_particle_spacing = 0.008
+fluid_particle_spacing = 0.005
 
 # ==========================================================================================
 # ==== Experiment Setup
 gravity = 9.81
-nu = 0.01
+nu = 0.001
 fluid_density = 1000.0
 sound_speed = 50
 
 sphere1_radius = 0.05
 
-sphere1_center = (0.5, 0.5, 0.2)
-sphere2_center = (1.5, 0.5, 0.2)
+sphere1_center = (0.5, 0.5, 0.075)
+sphere2_center = (1.5, 0.5, 0.075)
 sphere1 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere1_center,
-                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0))
+                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0))
 sphere2 = SphereShape(fluid_particle_spacing, sphere1_radius, sphere2_center,
-                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -2.0))
+                      fluid_density, sphere_type=VoxelSphere(), velocity=(0.0, 0.0, -1.0))
 
 # `compact_support` needs to be `2.0 * particle_spacing` to be correct
 fluid_smoothing_length = 1.0 * fluid_particle_spacing
 
 trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "falling_water_spheres_2d.jl"),
-              fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.2),
+              fluid_particle_spacing=fluid_particle_spacing, tspan=(0.0, 0.1),
               initial_fluid_size=(0.0, 0.0, 0.0),
               tank_size=(2.0, 1.0, 0.1), sound_speed=sound_speed,
               faces=(true, true, true, true, true, false),
               acceleration=(0.0, 0.0, -gravity), sphere1=sphere1, sphere2=sphere2,
               fluid_smoothing_length=fluid_smoothing_length,
               fluid_smoothing_kernel=SchoenbergCubicSplineKernel{3}(),
               nu=nu, alpha=10 * nu / (fluid_smoothing_length * sound_speed),
-              surface_tension_coefficient=1.5, adhesion_coefficient=0.5)
+              surface_tension_coefficient=10, adhesion_coefficient=0.1)