The grand reopening Part 3: Add basic DEM (#433)

* implement basic dem

* fix

* fix bnd

* get basic simulation working

* fix

* format

* fix spellcheck

* add docs

* update readme and news

* format

* add to examples run list

* add basic test

* ignore author name

* remove

* most of the review stuff

* format

* forgot to remove one thing

* one more

* rework interaction

* format

* update

* typo

* fix doctest

* fix test

* fix doc

* fix tests again

* fix doc and test

* fix test

* fix incorrect test

* review

* format

* typo

* review

* review

* format

* realign line

* fix test

* update

.typos.toml

@@ -2,4 +2,5 @@
 ba = "ba"
 Shepard = "Shepard"
 shepard = "shepard"
+Strack = "Strack"
 Lok = "Lok"

NEWS.md

@@ -4,8 +4,7 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](h
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability. 
 We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions. 
 
-
-## Version 0.1.x
+## Version 0.2.x
 
 ### Highlights
 
@@ -16,15 +15,28 @@ We aim at 3 to 4 month between major release versions and about 2 weeks between
 ### Deprecated
 
 
-## Pre Initial Release (v0.1.0)
-This section summarizes the initial features that TrixiParticles.jl was released with.
+## Version 0.1.2
 
 ### Highlights
-#### EDAC
+
+#### Discrete Element Method
+A basic implementation of the discrete element method was added.
+
+
+## Version 0.1.1
+
+### Added
+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
+
+# 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
+### WCSPH
 An implementation of WCSPH (Weakly Compressible Smoothed Particle Hydrodynamics), which is the classical SPH approach.
 
 Features:
@@ -36,5 +48,5 @@ Features:
 - Density diffusion based on the models by Molteni & Colagrossi (2009), Ferrari et al. (2009) and Antuono et al. (2010).
 
 
-#### TLSPH
+### TLSPH
 An implementation of TLSPH (Total Lagrangian Smoothed Particle Hydrodynamics) for solid bodies enabling FSI (Fluid Structure Interactions).

README.md

@@ -17,8 +17,9 @@ Its features include:
 ## Features
 - Incompressible Navier-Stokes
   - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC)
+  - Models: Surface Tension
 - Solid-body mechanics
-  - Methods:  Total Lagrangian SPH (TLSPH)
+  - Methods:  Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM)
 - Fluid-Structure Interaction
 - Output formats:
   - VTK

docs/make.jl

@@ -116,6 +116,8 @@ makedocs(sitename="TrixiParticles.jl",
                      "Util" => joinpath("general", "util.md"),
                  ],
                  "Systems" => [
+                     "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/boundary.md

@@ -4,6 +4,10 @@
     BoundarySPHSystem
 ```
 
+```@docs
+    BoundaryDEMSystem
+```
+
 ```@docs
     BoundaryMovement
 ```

docs/src/systems/dem.md

@@ -0,0 +1,33 @@
+# [Discrete Element Method](@id dem)
+The Discrete Element Method (DEM) is a computational technique widely used in physics, engineering,
+and applied mathematics for simulating the mechanical behavior of granular materials, such as powders,
+sand, soil, or rock, as well as other discontinua. Unlike continuum mechanics that treats materials as 
+continuous, DEM considers individual particles or elements and their interactions. This approach provides
+detailed insights into the micro-mechanical behavior of materials, making it particularly valuable
+in fields such as geomechanics, material science, and mechanical engineering.
+
+## Fundamental Principles
+The core idea behind DEM is the discretization of a material system into a finite set of distinct,
+interacting mass elements (particles). These elements (particles) can vary in shape, size, and properties, and 
+they interact with each other and possibly with their boundaries through contact forces and potential fields.
+The motion and behavior of each mass element are governed by Newton's laws of motion, accounting for the forces
+and moments acting upon them.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "solid", "discrete_element_method", "system.jl")]
+```
+
+## References
+- N. Bićanić. "Discrete element methods". 
+  In: Encyclopedia of Computational Mechanics (2007).
+  [doi: 10.1002/0470091355.ecm006.pub2](https://doi.org/10.1002/0470091355.ecm006.pub2)
+
+- P. Cundall and O. Strack. "A discrete numerical model for granular assemblies". 
+  In:  Géotechnique 29.1 (1979), pages 47--65.
+  [doi: 10.1680/geot.1979.29.1.47](https://doi.org/10.1680/geot.1979.29.1.47)
+
+- A. Renzo and F. Maio. "Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes"
+  In: Chemical Engineering Science 59.3 (2004), pages 525--541.
+  [doi: 10.1016/j.ces.2003.09.037](https://doi.org/10.1016/j.ces.2003.09.037)
+

examples/dem/rectangular_tank_2d.jl

@@ -0,0 +1,50 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+gravity = -9.81
+
+# ==========================================================================================
+# ==== Falling rocks
+
+particle_spacing = 0.1
+
+rock_width = 2.0
+rock_height = 2.0
+rock_density = 3000.0
+
+tank_width = 2.0
+tank_height = 4.0
+
+tank = RectangularTank(particle_spacing, (rock_width, rock_height),
+                       (tank_width, tank_height), rock_density,
+                       n_layers=2)
+
+# ==========================================================================================
+# ==== Systems
+
+# Move the rocks up to let them fall
+tank.fluid.coordinates[2, :] .+= 0.5
+rock_system = DEMSystem(tank.fluid, 2 * 10e5, 10e9, 0.3, acceleration=(0.0, gravity))
+boundary_system = BoundaryDEMSystem(tank.boundary, 10e7)
+
+# ==========================================================================================
+# ==== Simulation
+
+semi = Semidiscretization(rock_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch)
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=5000)
+saving_callback = SolutionSavingCallback(dt=0.02)
+
+callbacks = CallbackSet(info_callback, saving_callback)
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-3, # Limit stepsize to prevent crashing
+            dt=1e-7,  # Initial step size
+            save_everystep=false, callback=callbacks);

src/TrixiParticles.jl

@@ -42,7 +42,7 @@ include("visualization/recipes_plots.jl")
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
-       BoundarySPHSystem
+       BoundarySPHSystem, DEMSystem, BoundaryDEMSystem
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback
 export ContinuityDensity, SummationDensity

src/general/semidiscretization.jl

@@ -150,6 +150,15 @@ end
     return compact_support(smoothing_kernel, smoothing_length)
 end
 
+@inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem)
+    return 0.0
+end
+
+@inline function compact_support(system::BoundaryDEMSystem, neighbor::DEMSystem)
+    # Use the compact support of the DEMSystem
+    return compact_support(neighbor, system)
+end
+
 @inline function compact_support(system::TotalLagrangianSPHSystem,
                                  neighbor::TotalLagrangianSPHSystem)
     (; smoothing_kernel, smoothing_length) = system
@@ -578,6 +587,20 @@ end
 end
 
 # NHS updates
+# To prevent hard to spot errors there is not default version
+
+function nhs_coords(system::DEMSystem, neighbor::DEMSystem, u)
+    return current_coordinates(u, neighbor)
+end
+
+function nhs_coords(system::BoundaryDEMSystem,
+                    neighbor::Union{BoundaryDEMSystem, DEMSystem}, u)
+    return nothing
+end
+function nhs_coords(system::DEMSystem, neighbor::BoundaryDEMSystem, u)
+    return nothing
+end
+
 function nhs_coords(system::FluidSystem,
                     neighbor::FluidSystem, u)
     return current_coordinates(u, neighbor)

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -52,6 +52,7 @@ struct BoundaryModelDummyParticles{DC, ELTYPE <: Real, VECTOR, SE, K, V, COR, C}
                                          density_calculator, smoothing_kernel,
                                          smoothing_length; viscosity=nothing,
                                          state_equation=nothing, correction=nothing)
+        ELTYPE = eltype(initial_density)
         pressure = initial_boundary_pressure(initial_density, density_calculator,
                                              state_equation)
         NDIMS = ndims(smoothing_kernel)

src/schemes/boundary/rhs.jl

@@ -1,7 +1,8 @@
 # Interaction of boundary  with other systems
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system, neighborhood_search,
-                   particle_system::BoundarySPHSystem, neighbor_system)
+                   particle_system::BoundarySystem,
+                   neighbor_system)
     # TODO Solids and moving boundaries should be considered in the continuity equation
     return dv
 end

src/schemes/boundary/system.jl

@@ -40,6 +40,53 @@ struct BoundarySPHSystem{BM, NDIMS, IC, CO, M, IM, CA} <: BoundarySystem{NDIMS}
     end
 end
 
+"""
+    BoundaryDEMSystem(initial_condition, normal_stiffness)
+
+System for boundaries modeled by boundary particles.
+The interaction between fluid and boundary particles is specified by the boundary model.
+
+!!! warning "Experimental Implementation"
+    This is an experimental feature and may change in a future releases.
+
+"""
+struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D} <: BoundarySystem{NDIMS}
+    coordinates      :: ARRAY2D # [dimension, particle]
+    radius           :: ARRAY1D # [particle]
+    normal_stiffness :: ELTYPE
+
+    function BoundaryDEMSystem(initial_condition, normal_stiffness)
+        coordinates = initial_condition.coordinates
+        radius = 0.5 * initial_condition.particle_spacing *
+                 ones(length(initial_condition.mass))
+        NDIMS = size(coordinates, 1)
+
+        return new{NDIMS, eltype(coordinates), typeof(radius), typeof(coordinates)}(coordinates,
+                                                                                    radius,
+                                                                                    normal_stiffness)
+    end
+end
+
+function Base.show(io::IO, system::BoundaryDEMSystem)
+    @nospecialize system # reduce precompilation time
+
+    print(io, "BoundaryDEMSystem{", ndims(system), "}(")
+    print(io, system.boundary_model)
+    print(io, ") with ", nparticles(system), " particles")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", system::BoundaryDEMSystem)
+    @nospecialize system # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, system)
+    else
+        summary_header(io, "BoundaryDEMSystem{$(ndims(system))}")
+        summary_line(io, "#particles", nparticles(system))
+        summary_footer(io)
+    end
+end
+
 """
     BoundaryMovement(movement_function, is_moving; moving_particles=nothing)
 
@@ -119,7 +166,17 @@ function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem)
     end
 end
 
-timer_name(::BoundarySPHSystem) = "boundary"
+timer_name(::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary"
+
+@inline function Base.eltype(system::Union{BoundarySPHSystem, BoundaryDEMSystem})
+    eltype(system.coordinates)
+end
+
+# This does not account for moving boundaries, but it's only used to initialize the
+# neighborhood search, anyway.
+@inline function initial_coordinates(system::Union{BoundarySPHSystem, BoundaryDEMSystem})
+    system.coordinates
+end
 
 function (movement::BoundaryMovement)(system, t)
     (; coordinates, cache) = system
@@ -151,27 +208,28 @@ function (movement::Nothing)(system, t)
     return system
 end
 
-@inline function nparticles(system::BoundarySPHSystem)
-    length(system.boundary_model.hydrodynamic_mass)
+@inline function nparticles(system::Union{BoundaryDEMSystem, BoundarySPHSystem})
+    size(system.coordinates, 2)
 end
 
 # No particle positions are advanced for boundary systems,
 # except when using `BoundaryModelDummyParticles` with `ContinuityDensity`.
-@inline function n_moving_particles(system::BoundarySPHSystem)
+@inline function n_moving_particles(system::Union{BoundarySPHSystem, BoundaryDEMSystem})
     return 0
 end
 
 @inline function n_moving_particles(system::BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}})
     return nparticles(system)
 end
 
-@inline u_nvariables(system::BoundarySPHSystem) = 0
+@inline u_nvariables(system::Union{BoundarySPHSystem, BoundaryDEMSystem}) = 0
 
 # For BoundaryModelDummyParticles with ContinuityDensity, this needs to be 1.
 # For all other models and density calculators, it's irrelevant.
 @inline v_nvariables(system::BoundarySPHSystem) = 1
+@inline v_nvariables(system::BoundaryDEMSystem) = 0
 
-@inline function current_coordinates(u, system::BoundarySPHSystem)
+@inline function current_coordinates(u, system::Union{BoundarySPHSystem, BoundaryDEMSystem})
     return system.coordinates
 end
 
@@ -248,11 +306,13 @@ function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function write_u0!(u0, system::BoundarySPHSystem)
+function write_u0!(u0, system::Union{BoundarySPHSystem, BoundaryDEMSystem})
     return u0
 end
 
-function write_v0!(v0, system::BoundarySPHSystem)
+function write_v0!(v0,
+                   system::Union{BoundarySPHSystem,
+                                 BoundaryDEMSystem})
     return v0
 end
 

src/schemes/schemes.jl

@@ -3,6 +3,7 @@
 include("fluid/fluid.jl")
 include("boundary/boundary.jl")
 include("solid/total_lagrangian_sph/total_lagrangian_sph.jl")
+include("solid/discrete_element_method/discrete_element_method.jl")
 # Monaghan-Kajtar repulsive boundary particles require the `BoundarySPHSystem`
 # and the `TotalLagrangianSPHSystem`.
 include("boundary/monaghan_kajtar/monaghan_kajtar.jl")
@@ -12,3 +13,4 @@ include("fluid/weakly_compressible_sph/rhs.jl")
 include("fluid/entropically_damped_sph/rhs.jl")
 include("boundary/rhs.jl")
 include("solid/total_lagrangian_sph/rhs.jl")
+include("solid/discrete_element_method/rhs.jl")

src/schemes/solid/discrete_element_method/discrete_element_method.jl

@@ -0,0 +1 @@
+include("system.jl")