implement basic dem

examples/dem/rectangular_tank_2d.jl

@@ -0,0 +1,65 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+gravity = -9.81
+
+# ==========================================================================================
+# ==== Fluid
+
+particle_spacing = 0.1
+
+# Ratio of fluid particle spacing to boundary particle spacing
+beta = 1
+boundary_layers = 2
+
+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=boundary_layers, spacing_ratio=beta)
+
+# ==========================================================================================
+# ==== Boundary models
+
+boundary_model = BoundaryModelDummyParticles(tank.boundary.radius)
+
+# ==========================================================================================
+# ==== Systems
+
+# let them fall
+tank.fluid.coordinates[2,:] .+=0.5
+solid_system = DEMSystem(tank.fluid, 2*10^5, acceleration=(0.0, gravity))
+boundary_system = BoundaryDEMSystem(tank.boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Simulation
+
+semi = Semidiscretization(solid_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch)
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=50)
+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.
+# Enable threading of the RK method for better performance on multiple threads.
+# Limiting of the maximum stepsize is necessary to prevent crashing.
+# When particles are approaching a wall in a uniform way, they can be advanced
+# with large time steps. Close to the wall, the stepsize has to be reduced drastically.
+# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces
+# become extremely large when fluid particles are very close to boundary particles,
+# and the time integration method interprets this as an instability.
+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
+            save_everystep=false, callback=callbacks);

examples/fluid/rectangular_tank_2d.jl

@@ -24,7 +24,8 @@ state_equation = StateEquationCole(sound_speed, 7, water_density, 100000.0,
                                    background_pressure=100000.0,
                                    clip_negative_pressure=false)
 
-smoothing_length = 1.2 * particle_spacing
+smoothing_factor = 1.2
+smoothing_length = smoothing_factor * particle_spacing
 smoothing_kernel = SchoenbergCubicSplineKernel{2}()
 
 viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
@@ -39,7 +40,8 @@ tank = RectangularTank(particle_spacing, (water_width, water_height),
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
                                              AdamiPressureExtrapolation(),
-                                             smoothing_kernel, smoothing_length)
+                                             smoothing_kernel, particle_spacing,
+                                             smoothing_factor)
 
 # K = 9.81 * water_height
 # boundary_model = BoundaryModelMonaghanKajtar(K, beta, particle_spacing / beta,

src/TrixiParticles.jl

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

src/general/initial_condition.jl

@@ -63,6 +63,7 @@ struct InitialCondition{ELTYPE}
     mass             :: Array{ELTYPE, 1}
     density          :: Array{ELTYPE, 1}
     pressure         :: Array{ELTYPE, 1}
+    radius           :: Array{ELTYPE, 1}
 
     function InitialCondition(coordinates, velocities, masses, densities; pressure=0.0,
                               particle_spacing=-1.0)
@@ -82,8 +83,10 @@ struct InitialCondition{ELTYPE}
                                 "same length as `masses`"))
         end
 
+        radius = 0.5 * particle_spacing * ones(length(masses))
+
         return new{eltype(coordinates)}(particle_spacing, coordinates, velocities, masses,
-                                        densities, pressure)
+                                        densities, pressure, radius)
     end
 end
 

src/general/semidiscretization.jl

@@ -106,6 +106,15 @@ end
     return compact_support(smoothing_kernel, smoothing_length)
 end
 
+@inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem)
+    return 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
@@ -379,7 +388,10 @@ end
     return dv
 end
 
-@inline add_acceleration!(dv, particle, system::BoundarySPHSystem) = dv
+@inline function add_acceleration!(dv, particle,
+                                   system::Union{BoundarySPHSystem, BoundaryDEMSystem})
+    dv
+end
 
 @inline function add_damping_force!(dv, damping_coefficient, v, particle,
                                     system::FluidSystem)
@@ -420,6 +432,12 @@ function system_interaction!(dv_ode, v_ode, u_ode, semi)
 end
 
 # NHS updates
+
+# Default to update
+function nhs_coords(system, neighbor, u)
+    return current_coordinates(u, neighbor)
+end
+
 function nhs_coords(system::FluidSystem,
                     neighbor::FluidSystem, u)
     return current_coordinates(u, neighbor)

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -68,28 +68,52 @@ We provide three options to compute the boundary density and pressure, determine
 struct BoundaryModelDummyParticles{ELTYPE <: Real, SE, DC, K, V, C}
     pressure           :: Vector{ELTYPE}
     hydrodynamic_mass  :: Vector{ELTYPE}
+    radius             :: Vector{ELTYPE}
     state_equation     :: SE
     density_calculator :: DC
     smoothing_kernel   :: K
     smoothing_length   :: ELTYPE
+    smoothing_factor   :: ELTYPE
     viscosity          :: V
     cache              :: C
 
     function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,
                                          density_calculator, smoothing_kernel,
-                                         smoothing_length; viscosity=NoViscosity(),
+                                         particle_spacing,
+                                         smoothing_factor; viscosity=NoViscosity(),
                                          state_equation=nothing)
+        ELTYPE = eltype(initial_density)
         pressure = similar(initial_density)
 
         n_particles = length(initial_density)
 
+        radius = 0.5 * ones(ELTYPE, n_particles) * particle_spacing
+        smoothing_length = smoothing_factor * particle_spacing
+
         cache = (; create_cache(viscosity, n_particles, ndims(smoothing_kernel))...,
                  create_cache(initial_density, density_calculator)...)
 
-        new{eltype(initial_density), typeof(state_equation),
+        new{ELTYPE, typeof(state_equation),
             typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity),
-            typeof(cache)}(pressure, hydrodynamic_mass, state_equation, density_calculator,
-                           smoothing_kernel, smoothing_length, viscosity, cache)
+            typeof(cache)}(pressure, hydrodynamic_mass, radius, state_equation,
+                           density_calculator,
+                           smoothing_kernel, smoothing_length, smoothing_factor, viscosity,
+                           cache)
+    end
+
+    function BoundaryModelDummyParticles(radius)
+        ELTYPE = eltype(radius)
+        pressure = similar(radius)
+        n_particles = length(radius)
+
+        mass = ones(ELTYPE, n_particles)
+
+        smoothing_length = 2 * maximum(radius)
+
+        new{ELTYPE, Nothing, Nothing, Nothing, Nothing, Nothing}(pressure, mass, radius,
+                                                                 nothing, nothing, nothing,
+                                                                 smoothing_length, 1.0,
+                                                                 nothing, nothing)
     end
 end
 

src/schemes/boundary/rhs.jl

@@ -1,7 +1,7 @@
 # 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,
+                   particle_system::Union{BoundarySPHSystem, BoundaryDEMSystem},
                    neighbor_system)
     # TODO Solids and moving boundaries should be considered in the continuity equation
     return dv

src/schemes/boundary/system.jl

@@ -27,6 +27,46 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, M, C} <: System{NDIMS}
     end
 end
 
+"""
+    BoundaryDEMSystem(inititial_condition)
+
+System for boundaries modeled by boundary particles.
+The interaction between fluid and boundary particles is specified by the boundary model.
+
+"""
+struct BoundaryDEMSystem{BM, NDIMS, ELTYPE <: Real} <: System{NDIMS}
+    coordinates    :: Array{ELTYPE, 2}
+    boundary_model :: BM
+
+    function BoundaryDEMSystem(inititial_condition, model)
+        coordinates = inititial_condition.coordinates
+        NDIMS = size(coordinates, 1)
+
+        return new{typeof(model), NDIMS, eltype(coordinates)}(coordinates, model)
+    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_line(io, "boundary model", system.boundary_model)
+        summary_footer(io)
+    end
+end
+
 """
     BoundaryMovement(movement_function, is_moving)
 
@@ -90,13 +130,17 @@ function Base.show(io::IO, ::MIME"text/plain", system::BoundarySPHSystem)
     end
 end
 
-timer_name(::BoundarySPHSystem) = "boundary"
+timer_name(::Union{BoundarySPHSystem, BoundaryDEMSystem}) = "boundary"
 
-@inline Base.eltype(system::BoundarySPHSystem) = eltype(system.coordinates)
+@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 initial_coordinates(system::BoundarySPHSystem) = system.coordinates
+@inline function initial_coordinates(system::Union{BoundarySPHSystem, BoundaryDEMSystem})
+    system.coordinates
+end
 
 function (movement::BoundaryMovement)(system, t)
     (; coordinates, cache) = system
@@ -128,13 +172,13 @@ 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
 
@@ -152,13 +196,14 @@ end
     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
 
@@ -237,11 +282,13 @@ function update_final!(system::BoundarySPHSystem, system_index, v, u, v_ode, u_o
     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{<:BoundaryModelMonaghanKajtar})
+function write_v0!(v0,
+                   system::Union{BoundarySPHSystem{<:BoundaryModelMonaghanKajtar},
+                                 BoundaryDEMSystem})
     return v0
 end
 

src/schemes/schemes.jl

@@ -3,9 +3,11 @@
 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")
 
 # Include rhs for all schemes
 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")