diff --git a/examples/benchmarks/sd_static.cpp b/examples/benchmarks/sd_static.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c526bf1d43786d86f08d1695fc351d6c66b7117d
--- /dev/null
+++ b/examples/benchmarks/sd_static.cpp
@@ -0,0 +1,106 @@
+#include <iostream>
+#include <memory>
+
+#include "spring_dashpot.hpp"
+#include <chrono>
+
+void set_feature_properties(std::shared_ptr<PairsAccessor> &ac){
+ ac->setTypeStiffness(0,0, 1e6);
+ ac->syncTypeStiffness();
+
+ ac->setTypeDampingNorm(0,0, 300);
+ ac->syncTypeDampingNorm();
+
+ ac->setTypeFriction(0,0, 1.2);
+ ac->syncTypeFriction();
+
+ ac->setTypeDampingTan(0,0, 300);
+ ac->syncTypeDampingTan();
+}
+
+int main(int argc, char **argv) {
+ if(argc!=5){
+ std::cerr << "4 args are required: Domain size (i.e. number of particles) in x,y,z and #timesteps." << std::endl;
+ exit(-1);
+ }
+
+ double domain_size[3] = {std::stod(argv[1]), std::stod(argv[2]), std::stod(argv[3])};
+ uint64_t num_timesteps = std::stoull(argv[4]);
+
+ auto pairs_sim = std::make_shared<PairsSimulation>();
+ pairs_sim->initialize();
+
+ auto ac = std::make_shared<PairsAccessor>(pairs_sim.get());
+ set_feature_properties(ac); // Arbitray properties (all forces cancel out due to PBC in all dims)
+
+ auto pairs_runtime = pairs_sim->getPairsRuntime();
+
+ pairs_runtime->initDomain(&argc, &argv, 0, 0, 0, domain_size[0], domain_size[1], domain_size[2], false);
+
+ double particle_spacing = 1.0;
+
+ // Particle overlap is required for force calculation
+ double peneration_depth = 0.01;
+ double diameter = particle_spacing + peneration_depth;
+
+ double initial_velocity = 0.0; // Stationary
+ double density = 1000; // Arbitrary
+
+ pairs::dem_sc_grid(pairs_runtime, domain_size[0], domain_size[1], domain_size[2],
+ particle_spacing,
+ diameter, diameter, diameter,
+ initial_velocity, density, 1);
+
+
+ // Cell width is here smaller than sphere diameter only for convenience to have everything aligned on a grid, but this
+ // doesn't affect the interactions computed. All spheres are on cell centers and are in contact with 6 neighbors.
+ double cell_width = particle_spacing;
+ pairs_sim->setup_cells(cell_width, cell_width, cell_width, cell_width);
+
+ // Inertia update is required for euler updates to be valid (but particles remain stationary)
+ pairs_sim->update_mass_and_inertia();
+ double dt = 0.001; // Arbitrary
+
+ int rank = pairs_sim->rank();
+ if(rank==0) std::cout << "NUM_PROC: " << pairs_runtime->getDomainPartitioner()->getWorldSize() << std::endl;
+ if(rank==0) std::cout << "NUM_NEIGH_AABBS: " << pairs_runtime->getDomainPartitioner()->getNumberOfNeighborAABBs() << std::endl;
+ if(rank==0) std::cout << "NUM_TIMESTEPS: " << num_timesteps << std::endl;
+ int print_interval = (num_timesteps >= 5) ? (num_timesteps / 5) : 1;
+
+ // pairs::vtk_write_subdom(pairs_runtime, "output/subdom_init", pairs_runtime->getDomainPartitioner()->getWorldSize());
+
+
+ // ------------------------------------------------------------------------------
+ MPI_Barrier(MPI_COMM_WORLD);
+ auto start = std::chrono::high_resolution_clock::now();
+
+ for (int t=0; t<num_timesteps; ++t){
+ if ((t%print_interval==0) && rank==0) std::cout << "Timestep: " << t << std::endl;
+ pairs_sim->communicate(t);
+ pairs_sim->update_cells(t);
+ pairs_sim->spring_dashpot();
+ pairs_sim->euler(dt);
+ }
+
+ auto end = std::chrono::high_resolution_clock::now();
+ // ------------------------------------------------------------------------------
+
+
+ uint64_t nlocal = pairs_sim->nlocal();
+ uint64_t global_nparticles;
+ MPI_Reduce(&nlocal, &global_nparticles, 1, MPI_UINT64_T, MPI_SUM, 0, MPI_COMM_WORLD);
+ if(rank==0) {
+ auto duration = std::chrono::duration<double>(end - start);
+ double total_runtime = duration.count(); // seconds
+ std::cout << "TOTAL_RUNTIME: " << total_runtime << std::endl;
+ std::cout << "GLOBAL_NPARTICLES: " << global_nparticles << std::endl;
+
+ double pups = global_nparticles * num_timesteps / total_runtime; // particle updates per second
+ std::cout << "PUPS: " << pups << std::endl;
+ }
+
+ // if (rank==1) pairs::vtk_write_data(pairs_runtime, "output/ghost_spheres", pairs_sim->nlocal(), pairs_sim->size(), 0);
+ // if (rank==2) pairs::vtk_write_data(pairs_runtime, "output/ghost_spheres", pairs_sim->nlocal(), pairs_sim->size(), 0);
+
+ pairs_sim->end();
+}
\ No newline at end of file
diff --git a/examples/benchmarks/spring_dashpot.py b/examples/benchmarks/spring_dashpot.py
new file mode 100644
index 0000000000000000000000000000000000000000..bfc7aed159450cd8c3f91b47c94dbe2541aaf501
--- /dev/null
+++ b/examples/benchmarks/spring_dashpot.py
@@ -0,0 +1,96 @@
+import math
+import pairs
+import sys
+import os
+
+def update_mass_and_inertia(i):
+ rotation_matrix[i] = diagonal_matrix(1.0)
+ rotation[i] = default_quaternion()
+
+ if is_sphere(i):
+ inv_inertia[i] = inversed(diagonal_matrix(0.4 * mass[i] * radius[i] * radius[i]))
+
+ else:
+ mass[i] = infinity
+ inv_inertia[i] = 0.0
+
+def spring_dashpot(i, j):
+ delta_ij = -penetration_depth(i, j)
+ skip_when(delta_ij < 0.0)
+
+ velocity_wf_i = linear_velocity[i] + cross(angular_velocity[i], contact_point(i, j) - position[i])
+ velocity_wf_j = linear_velocity[j] + cross(angular_velocity[j], contact_point(i, j) - position[j])
+
+ rel_vel = -(velocity_wf_i - velocity_wf_j)
+ rel_vel_n = dot(rel_vel, contact_normal(i, j))
+ rel_vel_t = rel_vel - rel_vel_n * contact_normal(i, j)
+
+ fNabs = stiffness[i,j] * delta_ij + damping_norm[i,j] * rel_vel_n
+ fN = fNabs * contact_normal(i, j)
+
+ fTabs = min(damping_tan[i,j] * length(rel_vel_t), friction[i, j] * fNabs)
+ fT = fTabs * normalized(rel_vel_t)
+
+ partial_force = fN + fT
+ apply(force, partial_force)
+ apply(torque, cross(contact_point(i, j) - position, partial_force))
+
+def euler(i):
+ skip_when(is_fixed(i) or is_infinite(i))
+ inv_mass = 1.0 / mass[i]
+ position[i] += 0.5 * inv_mass * force[i] * dt * dt + linear_velocity[i] * dt
+ linear_velocity[i] += inv_mass * force[i] * dt
+ wdot = rotation_matrix[i] * (inv_inertia[i] * torque[i]) * transposed(rotation_matrix[i])
+ phi = angular_velocity[i] * dt + 0.5 * wdot * dt * dt
+ rotation[i] = quaternion(phi, length(phi)) * rotation[i]
+ rotation_matrix[i] = quaternion_to_rotation_matrix(rotation[i])
+ angular_velocity[i] += wdot * dt
+
+
+file_name = os.path.basename(__file__)
+file_name_without_extension = os.path.splitext(file_name)[0]
+
+psim = pairs.simulation(
+ file_name_without_extension,
+ [pairs.sphere()],
+ double_prec=True,
+ particle_capacity=10000000,
+ debug=True)
+
+target = sys.argv[1] if len(sys.argv[1]) > 1 else "none"
+if target == 'gpu':
+ psim.target(pairs.target_gpu())
+elif target == 'cpu':
+ psim.target(pairs.target_cpu())
+else:
+ print(f"Invalid target, use {sys.argv[0]} <cpu/gpu>")
+
+psim.add_position('position')
+psim.add_property('mass', pairs.real())
+psim.add_property('linear_velocity', pairs.vector())
+psim.add_property('angular_velocity', pairs.vector())
+psim.add_property('force', pairs.vector(), volatile=True)
+psim.add_property('torque', pairs.vector(), volatile=True)
+psim.add_property('radius', pairs.real())
+psim.add_property('inv_inertia', pairs.matrix())
+psim.add_property('rotation_matrix', pairs.matrix())
+psim.add_property('rotation', pairs.quaternion())
+
+ntypes = 1
+psim.add_feature('type', ntypes)
+psim.add_feature_property('type', 'stiffness', pairs.real())
+psim.add_feature_property('type', 'damping_norm', pairs.real())
+psim.add_feature_property('type', 'damping_tan', pairs.real())
+psim.add_feature_property('type', 'friction', pairs.real())
+
+# psim.set_domain_partitioner(pairs.regular_domain_partitioner())
+psim.set_domain_partitioner(pairs.block_forest())
+psim.pbc([True, True, True])
+psim.build_cell_lists(use_halo_cells=False)
+
+psim.compute(update_mass_and_inertia, symbols={'infinity': math.inf })
+psim.compute(spring_dashpot, profile=False)
+psim.compute(euler, parameters={'dt': pairs.real()})
+
+psim.generate()
+