Merge branch 'main' into update-callback

.github/workflows/Documenter.yml

@@ -28,7 +28,7 @@ jobs:
       - name: Check out project
         uses: actions/checkout@v4
       - name: Set up Julia
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: '1'
           show-versioninfo: true

.github/workflows/FormatCheck.yml

@@ -15,7 +15,7 @@ jobs:
       - name: Check out project
         uses: actions/checkout@v4
       - name: Set up Julia
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.19.0
+        uses: crate-ci/typos@v1.21.0

.github/workflows/ci.yml

@@ -61,7 +61,7 @@ jobs:
       - name: Check out project
         uses: actions/checkout@v4
       - name: Set up Julia
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
       - uses: julia-actions/cache@v1

.typos.toml

@@ -2,3 +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).

Project.toml

@@ -1,9 +1,10 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"]
-version = "0.1.0"
+version = "0.1.2-dev"
 
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

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/development.md

@@ -20,4 +20,49 @@ Alternatively, run
 ```bash
 python3 -m http.server -d docs/build
 ```
-and open `localhost:8000` in your webbrowser.
+and open `localhost:8000` in your webbrowser.
+
+
+## Release management
+
+To create a new release for TrixiParticles.jl, perform the following steps:
+1) Make sure that all PRs and changes that you want to go into the release are merged to
+   `main` and that the latest commit on `main` has passed all CI tests.
+2) Determine the currently released version of TrixiParticles.jl, e.g., on the
+   [release page](https://github.com/trixi-framework/TrixiParticles.jl/releases). For this manual,
+   we will assume that the latest release was `v0.2.3`.
+3) Decide on the next version number. We follow [semantic versioning](https://semver.org/),
+   thus each version is of the form `vX.Y.Z` where `X` is the major version, `Y` the minor
+   version, and `Z` the patch version. In this manual, we assume that the major version is
+   always `0`, thus the decision process on the new version is as follows:
+   * If the new release contains *breaking changes* (i.e., user code might not work as
+     before without modifications), increase the *minor* version by one and set the
+     *patch* version to zero. In our example, the new version should thus be `v0.3.0`.
+   * If the new release only contains minor modifications and/or bug fixes, the *minor*
+     version is kept as-is and the *patch* version is increased by one. In our example, the
+     new version should thus be `v0.2.4`.
+4) Edit the `version` string in the
+   [`Project.toml`](https://github.com/trixi-framework/TrixiParticles.jl/blob/main/Project.toml)
+   and set it to the new version. Push/merge this change to `main`.
+5) Go to GitHub and add a comment to the commit that you would like to become the new
+   release (typically this will be the commit where you just updated the version). You can
+   comment on a commit by going to the
+   [commit overview](https://github.com/trixi-framework/TrixiParticles.jl/commits/main/) and clicking
+   on the title of the commit. The comment should contain the following text:
+   ```
+   @JuliaRegistrator register
+   ```
+6) Wait for the magic to happen! Specifically, JuliaRegistrator will create a new PR to the
+   Julia registry with the new release information. After a grace period of ~15 minutes,
+   this PR will be merged automatically. A short while after,
+   [TagBot](https://github.com/trixi-framework/TrixiParticles.jl/blob/main/.github/workflows/TagBot.yml)
+   will create a new release of TrixiParticles.jl in our GitHub repository.
+7) Once the new release has been created, the new version can be obtained through the Julia
+   package manager as usual.
+8) To make sure people do not mistake the latest state of `main` as the latest release, we
+   set the version in the `Project.toml` to a *development* version. The development version
+   should be the latest released version, with the patch version incremented by one, and the
+   `-dev` suffix added. For example, if you just released `v0.3.0`, the new development
+   version should be `v0.3.1-dev`. If you just released `v0.2.4`, the new development
+   version should be `v0.2.5-dev`.
+

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

@@ -2,6 +2,7 @@ module TrixiParticles
 
 using Reexport: @reexport
 
+using Adapt: Adapt
 using CSV: CSV
 using Dates
 using DataFrames: DataFrame
@@ -34,15 +35,16 @@ include("setups/setups.jl")
 include("schemes/schemes.jl")
 
 # Note that `semidiscretization.jl` depends on the system types and has to be
-# included separately.
+# included separately. `gpu.jl` in turn depends on the semidiscretization type.
 include("general/semidiscretization.jl")
+include("general/gpu.jl")
 include("visualization/write2vtk.jl")
 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, UpdateCallback
 export ContinuityDensity, SummationDensity

src/general/gpu.jl

@@ -0,0 +1,34 @@
+# Adapt.jl provides a function `adapt(to, x)`, which adapts a value `x` to `to`.
+# In practice, this means that we can use `adapt(CuArray, system)` to adapt a system to
+# the `CuArray` type.
+# What this does is that it converts all `Array`s inside this system to `CuArray`s,
+# therefore copying them to the GPU.
+# In order to run a simulation on a GPU, we want to call `adapt(T, semi)` to adapt the
+# `Semidiscretization` `semi` to the GPU array type `T` (e.g. `CuArray`).
+#
+# `Adapt.@adapt_structure` automatically generates the `adapt` function for our custom types.
+Adapt.@adapt_structure Semidiscretization
+Adapt.@adapt_structure WeaklyCompressibleSPHSystem
+Adapt.@adapt_structure DensityDiffusionAntuono
+Adapt.@adapt_structure BoundarySPHSystem
+Adapt.@adapt_structure BoundaryModelDummyParticles
+Adapt.@adapt_structure BoundaryModelMonaghanKajtar
+
+# The initial conditions are only used for initialization, which happens before `adapt`ing
+# the semidiscretization, so we don't need to store `InitialCondition`s on the GPU.
+# To save precious GPU memory, we replace initial conditions by `nothing`.
+function Adapt.adapt_structure(to, ic::InitialCondition)
+    return nothing
+end
+
+# `adapt(CuArray, ::SVector)::SVector`, but `adapt(Array, ::SVector)::Vector`.
+# We don't want to change the type of the `SVector` here.
+function Adapt.adapt_structure(to::typeof(Array), svector::SVector)
+    return svector
+end
+
+# `adapt(CuArray, ::UnitRange)::UnitRange`, but `adapt(Array, ::UnitRange)::Vector`.
+# We don't want to change the type of the `UnitRange` here.
+function Adapt.adapt_structure(to::typeof(Array), range::UnitRange)
+    return range
+end