PeriDEM

PeriDEM - High-fidelity modeling of granular media consisting of deformable complex-shaped particles

Table of contents

• Introduction
• Documentation
• Examples
• Brief implementation details
• Installation
  • Dependencies
  • Building the code
  • Install & use as a CMake package
  • Parallelism (MPI)
• Running simulations
  • Deck layout (JSON)
  • Setting up a problem in Python
  • Two-particle contact
  • Compressive test
  • Attrition
  • Impact and fracture
  • Single-particle Peridynamics
• Visualizing results
• Contributing
• Citations
• Developers

Introduction

Implementation of the high-fidelity model of granular media that combines the advantages of peridynamics and the discrete element method (DEM). The model has the following advantages over existing mechanical models for granular media:

• handle intra-particle deformation and breakage/damage
• handle the arbitrary shape of the particle. Inter-particle contact is not specific to any shape of the particle
• tunable inter-particle contact parameters
• easy to add different mechanical constitutive laws within peridynamics for individual particle deformation

For more details about the model and results, we refer to the paper:

Prashant K. Jha, Prathamesh S. Desai, Debdeep Bhattacharya, Robert P Lipton (2020). Peridynamics-based discrete element method (PeriDEM) model of granular systems involving breakage of arbitrarily shaped particles. Journal of the Mechanics and Physics of Solids, 151, p.104376. Doi https://doi.org/10.1016/j.jmps.2021.104376. Download pdf here.

PeriDEM is published as a software article in the Journal of Open Source Software:

Prashant K. Jha (2025). PeriDEM – High-fidelity modeling of granular media consisting of deformable complex-shaped particles Journal of Open Source Software, vol. 10, 116, p.7525, DOI 10.21105/joss.07525. Download pdf here.

We have created these channels:

• PeriDEM on Gitter
  • Gitter is absolutely open and easy to join.
• PeriDEM on slack
  • Email us if interested in joining the workspace.

Documentation

Doxygen generated documentation details functions and objects in the library.

Examples

We next highlight some key examples. Further details are available in examples/README.md.

Two-particle tests

Circular without damping	Circular with damping

Different materials	Different radius	Different radius different material

Two-particle with wall test

Concave particles

Compressive tests

Paper setup (Jha et al. 2021): 502 circular and hexagonal particles in a rectangle container; the top wall moves downward at fixed speed. Reaction on the moving wall rises with penetration; damage then concentrates along force chains and the pack yields. Runnable decks: small pack n12 and paper-scale two-stage n500. Details: Jha et al. 2021.

Pack geometry (N≈502)	Wall reaction and damage frames

Compressive test simulation

Attrition tests

Mix of circular, triangular, hexagonal, and drum-shaped grains in a rotating container (size and toughness vary). Portable JSON decks:

Rotating cylinder (setup)	Thin container, offset rotation (setup)

Impact and fracture

Silling KW 3D notched plate (setup, ./run_3d.sh)

Single particle deformation

Model.Particle_Sim_Type = Single_Particle. JSON demos: examples/Peridynamics.

Circle (setup)	Rectangle / CreateMesh (setup)

Brief implementation details

The simulation driver is class PeriDEMModel in PeriDEM/. Libraries live under src/. PeriDEMModel::run() initializes the simulation, optionally restarts, then hands the time loop to time_int::Integrator.

PeriDEMModel::run()

void PeriDEMModel::run(std::shared_ptr<inp::Input> &deck) {
    init();
    if (d_modelDeck_p->d_isRestartActive)
      restart(deck);
    integrate();  // time_int::Integrator().integrate(*this)
    close();
}

init() creates particles, sets up contact and quadrature data, builds neighbor lists and peridynamic bonds, and initializes loading.

Time integration

PeriDEMModel::integrate() calls time_int::Integrator. The integrator applies initial conditions, displacement BCs, and forces, then advances with central difference or velocity Verlet using data::ModelData kinematics accessors. After each step it writes output and calls checkStop().

PeriDEMModel::computeForces()

void PeriDEMModel::computeForces() {
    // reset nodal force
    pd::computeForces(*this);
    if (multi-particle)
      d_contact_p->computeForces(*this);
    computeExternalForces();
}

Contact::computeForces loops over the contact neighbour list. The node-node relation is contact::PairForce; damping is contact::Damping. A different pair law is a PairForce subclass set with Contact::setPairForce. Do not copy contact.cpp.

Further reading

See periDEMModel.cpp, src/time_int/integrator.h, and src/contact.

Installation

The pixi.toml file defines the dependencies and build instructions for the library. It should be used to create a reproducible environment and build the code using Pixi and CMake.

To install the Pixi package manager, follow the instructions at the official installation page.

Dependencies

Core dependencies are:

• cmake (>= 3.10.2)
• vtk (>= 7.1.1)
• metis (>= 5.1.0)
• MPI

Following dependencies are included in the PeriDEM library in external folder (see external/README.md for more details):

• fast-cpp-csv-parser (version included - master)
• nanoflann (>= 1.3.2, version included - v1.5.5)
• taskflow (>= 3.7.0)
• doxygen-awesome-css (>= v2.3.3)
• nlohmann_json (>= 3.12.0)

Pixi: Easiest method to build the library on ubuntu and mac

# run ubuntu using docker (we are using the same image we use to test the library)
docker run -it prashjha/peridem-base-noble

# we install pixi and add it to the path
curl -fsSL https://pixi.sh/install.sh | sh
export PATH="/root/.pixi/bin:$PATH"

# assuming we are now in root of docker image
cd user/
git clone [email protected]:prashjha/PeriDEM.git
cd PeriDEM/
pixi run test

Building the code

If all the dependencies are installed on the global path (e.g., /usr/local/), commands for building the PeriDEM code is as simple as

cmake   -DEnable_Documentation=OFF \
        -DEnable_Tests=ON \
        -DEnable_High_Load_Tests=OFF \
        -DDisable_Docker_MPI_Tests=ON \
        -DVTK_DIR="${VTK_DIR}" \
        -DMETIS_DIR="${METIS_DIR}" \
        -DCMAKE_BUILD_TYPE=Release \
        <PeriDEM source directory>
        
make -j 4

:exclamation: cmake and make commands should be run inside the build directory. You can create the build directory either inside or outside the repository.

Python

-DEnable_Python=ON builds an in-process import peridem (nanobind). It is a library interface and does not call bin/PeriDEM as a subprocess: you build the problem, run it, and read the node fields as numpy arrays in the same process.

cmake -S . -B build -DEnable_Python=ON
cmake --build build --target peridem_core PeriDEM -j
export PYTHONPATH="$PWD/build/python:$PYTHONPATH"

# or, to install it
pip install -e . --no-build-isolation

See Setting up a problem in Python below and python/README.md for the full API.

Install & use as a CMake package

• Build and install (starting from a fresh clone, e.g., git clone ... && cd PeriDEM; create a build directory wherever you like, build inside the source is assumed below):

  # from the source root
  mkdir -p build
  cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
  cmake --build build -- -j$(sysctl -n hw.ncpu)
  
  # install library in /tmp/peridem-install
  cmake --install build --prefix /tmp/peridem-install
  

  This installs bin/PeriDEM (source: PeriDEM/), shared libs in lib/, headers in include/, and the CMake package files under lib/cmake/PeriDEM.
• Consume in another CMake project:

  cmake_minimum_required(VERSION 3.18)
  project(peridem_consumer LANGUAGES CXX)
  find_package(PeriDEM REQUIRED)
  add_executable(hello main.cpp)
  target_link_libraries(hello PRIVATE PeriDEM::PeriDEMModel)
  

  Place this in, e.g., /tmp/peridem-consumer/CMakeLists.txt. A minimal main.cpp in the same folder:

  #include <iostream>
  #include <PeriDEMConfig.h>
  int main() {
    std::cout << "PeriDEM version: "
              << PERIDEM_VERSION_MAJOR << "."
              << PERIDEM_VERSION_MINOR << "."
              << PERIDEM_VERSION_PATCH << "\n";
  }
  

  Configure and build the consumer (run these inside /tmp/peridem-consumer):

  cmake -S . -B build -DCMAKE_PREFIX_PATH=/tmp/peridem-install
  cmake --build build -- -j$(sysctl -n hw.ncpu)
  ./build/hello
  

• External dependencies required on the target system: MPI, Threads, VTK (CommonCore/DataModel/IOXML), BLAS/LAPACK (Accelerate on macOS), Metis and Gmsh (found via bundled FindMetis.cmake and FindGmsh.cmake), plus their transitive libraries. Ensure these are installed and discoverable (e.g., via CMAKE_PREFIX_PATH or system paths) when configuring consumers. Bundled headers (nlohmann_json, nanoflann, csv, taskflow) are installed with the package.

Parallelism (MPI)

Decks are independent of MPI mode. Set Model.MPI_Strategy in the JSON input:

Value	Meaning
auto (default)	Multi_Particle → Particle-MPI; Single_Particle → DOF-MPI
none	No domain split (mpirun -n 1)
particle	Particle-MPI: distribute whole particles across ranks
dof	DOF-MPI: distribute nodes/DOFs across ranks

• Particle-MPI: each rank owns whole particles; near-contact and wall neighbors are exchanged as ghosts.
• DOF-MPI: each rank owns a subset of nodes (any body, including packs with walls). Before contact, Multi_Particle syncs owned nodal u/v to every rank.
• Threads (-nThreads) combine with MPI.

Identity checks (serial vs particle@2 vs dof@2) live under test/test_data/peridem/twop_circ_inbuilt/, jha2021_comp_n50/, and mpi_identity_twop_wall/.

Not covered here: GPU offload; larger-scale weak scaling.

Ask for help

Earlier releases depended on large libraries such as HPX, PCL, and Boost. Those are gone. Current configure needs VTK, MPI, Metis, Gmsh (for built-in meshing / some tests), and BLAS/LAPACK (Accelerate on macOS), plus a C++20 toolchain. Use pixi.toml (see Pixi above) or the CMake steps in Building the code.

Feel free to reach out or open an issue. For more open discussion of issues and ideas, contact via PeriDEM on Gitter or PeriDEM on slack (for slack, email us to join). If you like some help, want to contribute, extend the code, or discuss new ideas, please do reach out to us.

Running simulations

Input is JSON only (bin/PeriDEM -i input.json). Mesh files (.msh) and particle-location CSVs are referenced from the deck. Example:

<path of PeriDEM>/bin/PeriDEM -i input.json -nThreads 4
# or with MPI
mpirun -n 4 <path of PeriDEM>/bin/PeriDEM -i input.json -nThreads 2

Most example folders provide ./run.sh (or run_stage1.sh / run_stage2.sh) that locate bin/PeriDEM under build/. Index: examples/README.md.

The same problems can be set up and run from Python. Every example folder has a problem.py that builds its deck through the peridem interface. See Setting up a problem in Python.

Deck layout (JSON)

A multi-particle deck has these top-level blocks:

Block	Role
Model	Dimension, time, Particle_Sim_Type (Multi_Particle / Single_Particle), MPI_Strategy
Particle / Mesh / Material	Geometry sets, meshes (File or CreateMesh), PD material
Displacement_BC / Force_BC	Regions, directions, time/space functions
Contact / Neighbor	Inter-particle (and wall) contact; neighbor list
Particle_Generation	Pack / container / wall placement when not listing every body by hand
Output	Path, tags (Displacement, Velocity, Damage_Z, …), optional PVD_Collection
Restart	Optional settled IC for two-stage runs

Copy a short deck from examples/PeriDEM/compressive/n12/ or examples/Peridynamics/circle/ and change geometry, BCs, and time. Full block details: Doxygen and the checked-in example JSON files.

Setting up a problem in Python

With -DEnable_Python=ON, a problem can be built and run entirely in Python. No JSON file and no .msh need exist: geometry, mesh, material, contact, boundary conditions and particle placement are all set through the API, and the mesh is generated in-process by Gmsh.

import peridem
from peridem import Deck, Geometry
from peridem.deck import MeshSpec, contact_stiffness

peridem.init(n_threads=4)

R, h = 0.001, 0.0002                       # grain radius, mesh size
horizon = 3.0 * h

d = Deck(dim=2, t_final=0.012, n_steps=36000)
d.set_output("runs/out", tags=["Displacement", "Velocity", "Damage_Z"],
             interval=3600)
d.set_gravity(0.0, -10.0)

grain = d.add_particle_type(Geometry("circle", [R, 0, 0, 0]),
                            MeshSpec(size=h))          # Gmsh, in memory
steel = d.add_material(horizon=horizon, density=1200.0, K=2.16e7, nu=0.25,
                       Gc=50.0, influence_fn_type=1)
d.add_contact_pair(0, 0, Kn=contact_stiffness(2.16e7, 2.16e7, horizon),
                   K=2.16e7, eps=0.9)

d.add_displacement_bc(particles=[0], direction=[1, 2], zero_displacement=True)
d.add_initial_velocity([0.0, -0.089, 0.0], particles=[1])  # free fall, 0.4 mm

d.place(R, R, geometry=grain, material=steel)           # fixed grain
d.place(R, 4 * R, geometry=grain, material=steel)       # falling grain

sim = d.run()

u = sim.displacement                        # numpy (N, 3) view into the model
print(sim.n_particles, "grains,", sim.n_nodes, "nodes at t =", sim.time)
for p in sim.particles:
    print(p.id, p.geometry_name, p.center_of_mass, float(p.damage.max()))

2 grains, 246 nodes at t = 0.012
0 circle [0.001, 0.001, 0.0] 0.0
1 circle [0.001, 0.003797, 0.0] 0.001737

d.write("input.json") writes the same deck out, so a problem prototyped in Python can be handed to bin/PeriDEM unchanged. Conversely peridem.Simulation.from_file("input.json") runs an existing deck in-process, and python -m peridem -i input.json -nThreads 4 does it from the shell (mpirun works on both).

Beyond deck building, the interface exposes the model itself: per-particle node slices (sim.particle(i).velocity), writable field views, the peridynamic bond graph (sim.break_bonds_in_slots(...), used to seed a pre-notch), and the individual steps of the time loop so it can be driven from Python. Full API: python/README.md.

Examples in Python

Every example folder has a problem.py that sets that problem up through this interface. run.py next to it runs that file:

cd examples/PeriDEM/ellipse_triangle
./run.py                                  # in-process, Python-authored deck
./problem.py --num-steps 2000 --snapshot crack.png
./problem.py --write-deck /tmp/input.json # then: bin/PeriDEM -i /tmp/input.json

Where a specific mesh is needed, which is the case for the attrition packs and for runs compared against archived C++ output, problem.py reads the committed .msh. Otherwise it generates the mesh in the calling process and writes no file.

Checking that both paths agree

The Python side is a binding, not a second implementation. peridem.decks.* calls the same inp::*Deck::getExampleJson factories as the C++ drivers, and peridem.to_E, peridem.contact_stiffness, Geometry and the bond-cutting calls are the C++ material::, util::, geom:: and geometry:: functions.

python/tests/test_example_parity.py builds each example’s deck in Python and compares it key by key against the deck the matching C++ driver writes, then runs that deck through bin/PeriDEM and through the in-process interface and compares Displacement, Velocity, Force, Damage_Z and Damage node by node. Neither side reads a deck the other produced. The cases run at reduced step counts; pass --full to include the two largest.

PYTHONPATH=build/python python3 python/tests/test_example_parity.py
ctest --test-dir build -R Test_Python -V

Naming deck values

Deck values are given by name on both sides. Deck(2, 1.0, 10) and decks.model(2, 1.0, 10) are rejected:

d = Deck(dim=2, t_final=1.0, n_steps=10)

Each deck declares its fields once, as a table of key, type, default and description, from which its reader, writer, printer and schema are generated. A key the table does not declare is rejected where the deck is built, with the closest declared name given. A quantity that can be given in more than one way, such as a horizon as a length or as a multiple of the mesh size, is selected by naming the key that applies. decks.schema(name="model") returns the table for one deck.

Two-particle contact

JSON via bin/PeriDEM: start from a compressive or attrition short deck. C++ driver (shares the twop inbuilt test): examples/PeriDEM/twop_circ_contact.

cd examples/PeriDEM/twop_circ_contact
./run.py                          # Python; both discs meshed in-process

Compressive test

Path	Role
examples/PeriDEM/compressive/n12	Small 4×3 pack; short / MPI identity decks
examples/PeriDEM/compressive/n500	Paper N≈502 two-stage settle → compress

cd examples/PeriDEM/compressive/n12
./run.sh
DECK=input_quick_dof.json NP=4 ./run.sh

cd examples/PeriDEM/compressive/n500
NP=4 ./run_stage1.sh              # or use checked-in settled restart
NP=1 ./run_stage2.sh

In Python the pack, the U-channel cup and the moving plate are derived from R and the pack size rather than read from a deck:

cd examples/PeriDEM/compressive/n12
./problem.py --ncols 4 --nrows 3 --in-process-mesh --num-steps 2000

Attrition

Path	Role
attrition/sim1_rotating_cylinder	Thick rotating drum
attrition/sim2_thin_container	Thin drum, offset rotation

Each folder has ./run.sh (and mesh/CSV setup scripts). Keep outputs under runs/ (gitignored).

cd examples/PeriDEM/attrition/sim1_rotating_cylinder
./problem.py --preset short --snapshot damage.png

problem.py builds the drum as a complex geometry (outer circle, inner circle removed, protrusion added) and places it at its signed-volume centroid, which the geometry object computes. The grain packs come from the committed CSV and .msh files so the runs stay comparable with the archived results.

Impact and fracture

Path	Role
silling_kw	Silling KW 3D (./run_3d.sh) and 2D (./run_2d.sh)
ellipse_triangle	Hollow ellipse dropped on a tip

The Kalthoff–Winkler case needs more than a deck: after the model is built, the peridynamic bonds spanning each notch slot have to be cut. That is a Python call, and the time loop can then be driven step by step to record when each node first becomes damaged, which is how the crack speed is measured.

cd examples/PeriDEM/silling_kw
./problem.py --quick                      # smaller, softer plate
./problem.py --arrival-times              # 2D, reports crack speed
./problem.py --dim3                       # Silling's 200x100x9 mm plate

Single-particle Peridynamics

Path	Role
Peridynamics/circle	File mesh; fixed / pull BC
Peridynamics/rectangle	In-process CreateMesh; fixed / pull BC

cd examples/Peridynamics/circle
./run.sh                          # short deck
DECK=input.json NP=2 ./run.sh     # full; auto → DOF-MPI on multi-rank

./run.py                          # same problem, set up in Python
./problem.py --mesh-size 3e-4     # ... meshing in-process instead

cd ../rectangle
./problem.py --snapshot pull.png  # uniform grid built in memory

Visualizing results

Simulation files output_*.vtu (and output.pvd when PVD_Collection is on) can be loaded in either ParaView or VisIt.

By default, in all tests and examples, we only output the particle mesh, i.e., a pair of nodal coordinate and nodal volume, and not the finite element mesh (it can be enabled by setting Perform_FE_Out to true within the Output block in the JSON deck). After loading the file in ParaView, the first thing to do is to change the plot type from Surface to Point Gaussian. Next, a couple of things to do are:

• Adjust the radius of circle/sphere at the nodes by going to the Properties tab on the left side and change the value of Gaussian Radius
• You may also want to choose the field to display. For starter, you could select the Damage_Z variable, a ratio of maximum bond strain in the neighborhood of a node and critical bond strain. When the Damage_Z value is below one at a given node, the deformation in the vicinity of that node is elastic, whereas when the value is above 1, it indicates there is at least one node in the neighborhood which has bond strain above critical strain (meaning the bond between these two nodes is broken)
• You may also need to rescale the plot by clicking on the Zoom to Data button in ParaView
• Lastly, when the Damage_Z is very high at few nodes, you may want to rescale the data to the range, say [0,2] or [0,10], so that it is easier to identify regions with elastic deformation and region with fracture.

Contributing

We welcome contributions to the code. Limitations under Parallelism (MPI) are noted there. Please fork this repository, make changes, and make a pull request to the source branch.

Citations

If this library was useful in your work, we recommend citing the following article:

Jha, P.K., Desai, P.S., Bhattacharya, D. and Lipton, R., 2021. Peridynamics-based discrete element method (PeriDEM) model of granular systems involving breakage of arbitrarily shaped particles. Journal of the Mechanics and Physics of Solids, 151, p.104376.

You can also cite the PeriDEM using zenodo doi:

Prashant K., J. (2024). Peridynamics-based discrete element method (PeriDEM) model of granular systems. Zenodo. https://doi.org/10.5281/zenodo.13888588

Developers

• Prashant K. Jha ([email protected], [email protected])