Add Hertz contact rebound benchmark

Cargo.toml

@@ -30,6 +30,10 @@ path = "examples/granular_basic/main.rs"
 name = "granular_gas_benchmark"
 path = "examples/granular_gas_benchmark/main.rs"
 
+[[example]]
+name = "bench_hertz_rebound"
+path = "examples/bench_hertz_rebound/main.rs"
+
 [[example]]
 name = "toml_single"
 path = "examples/toml_single/main.rs"

examples/bench_hertz_rebound/.gitignore

@@ -0,0 +1,7 @@
+# Generated output — not tracked
+data/
+sweep/
+plots/
+output*/
+dump/
+GranularTemp.txt

examples/bench_hertz_rebound/README.md

@@ -0,0 +1,96 @@
+# Hertz Contact Rebound Benchmark
+
+Validates Hertzian contact mechanics by dropping a single sphere onto a rigid flat wall and measuring the coefficient of restitution (COR), contact duration, and peak overlap.
+
+## Physics
+
+A sphere of radius R, mass m, impacts a rigid flat wall at velocity v₀. The Hertz contact model predicts:
+
+- **Contact duration** (elastic):
+  ```
+  t_c = 2.87 × (m²/(R·E*²·v₀))^(1/5)
+  ```
+- **Peak overlap** (elastic):
+  ```
+  δ_max = (15·m·v₀² / (16·√R·E*))^(2/5)
+  ```
+- **COR**: The viscoelastic damping model should reproduce the input COR parameter.
+
+where E* = E/(2(1−ν²)) is the reduced modulus for sphere-on-flat contact.
+
+## Material Properties
+
+| Property | Value | Unit |
+|----------|-------|------|
+| Young's modulus E | 70 GPa | Pa |
+| Poisson's ratio ν | 0.22 | — |
+| Density ρ | 2500 | kg/m³ |
+| Radius R | 5 | mm |
+
+## Parameter Sweep
+
+- **Impact velocities**: 0.1, 0.5, 1.0, 2.0 m/s
+- **COR values**: 0.5, 0.7, 0.9, 0.95
+
+## Validation Criteria
+
+| Check | Tolerance | Notes |
+|-------|-----------|-------|
+| COR matches input (COR ≥ 0.7) | ≤ 3% relative error | |
+| COR matches input (COR < 0.7) | ≤ 12% relative error | Known Hertz nonlinearity effect* |
+| Contact duration vs Hertz | ≤ 10% relative error | |
+| Peak overlap vs Hertz | ≤ 10% relative error | |
+| All 16 cases complete | 16/16 | |
+
+\* The β damping coefficient is derived from linear (Hooke) contact theory. When applied with nonlinear Hertz stiffness, the achieved COR deviates from the input value, especially at low COR. This is a well-known limitation shared by LAMMPS and other DEM codes using the same model.
+
+## How to Run
+
+### Single case (default config)
+
+```bash
+cargo run --release --example bench_hertz_rebound --no-default-features -- examples/bench_hertz_rebound/config.toml
+```
+
+### Full parameter sweep
+
+```bash
+python3 examples/bench_hertz_rebound/run_sweep.py
+```
+
+### Validate results
+
+```bash
+python3 examples/bench_hertz_rebound/validate.py
+```
+
+### Generate plots
+
+```bash
+python3 examples/bench_hertz_rebound/plot.py
+```
+
+## Expected Plots
+
+### COR Validation
+![COR validation](plots/cor_validation.png)
+
+### Contact Duration
+![Contact duration](plots/contact_duration.png)
+
+### Peak Overlap
+![Peak overlap](plots/peak_overlap.png)
+
+## Assumptions
+
+- **3D simulation** with a single spherical particle
+- **No friction** (friction = 0) for clean normal-only rebound
+- **No gravity** effect on contact (gravity is off; particle given direct velocity)
+- **Monodisperse** — single particle size
+- **Hertz–Mindlin** contact model with viscoelastic damping (MDDEM default)
+- Wall is treated as **infinitely massive and rigid** (standard DEM wall)
+
+## References
+
+1. K.L. Johnson, *Contact Mechanics*, Cambridge University Press, 1985.
+2. L. Vu-Quoc and X. Zhang, "An accurate and efficient tangential force-displacement model for elastic frictional contact in particle-flow simulations", *Mechanics of Materials*, 31(4):235–269, 1999.

examples/bench_hertz_rebound/config.toml

@@ -0,0 +1,60 @@
+# Hertz Contact Rebound Benchmark
+# Single glass sphere dropped onto a rigid flat wall.
+# Validates: COR, contact duration, peak overlap against Hertz theory.
+#
+# Run: cargo run --release --example bench_hertz_rebound --no-default-features -- examples/bench_hertz_rebound/config.toml
+
+[comm]
+processors_x = 1
+processors_y = 1
+processors_z = 1
+
+[domain]
+# Simulation box [m] — small box, single particle
+x_low = -0.01
+x_high = 0.01
+y_low = -0.01
+y_high = 0.01
+z_low = 0.0
+z_high = 0.1
+periodic_x = false
+periodic_y = false
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1
+bin_size = 0.015          # Bin size [m] — larger than particle diameter
+every = 1
+
+[[dem.materials]]
+name = "glass"
+youngs_mod = 70.0e9       # Young's modulus [Pa] — soda-lime glass
+poisson_ratio = 0.22      # Poisson's ratio [dimensionless]
+restitution = 0.9         # Coefficient of restitution [0–1]
+friction = 0.0            # No friction for clean normal-only rebound test
+
+[[particles.insert]]
+material = "glass"
+count = 1                 # Single particle
+radius = 0.005            # Particle radius [m] = 5 mm
+density = 2500.0          # Particle density [kg/m³] — glass
+velocity_z = -1.0         # Impact velocity [m/s] — downward
+# Place particle just above the wall so it impacts quickly (z = radius + small gap)
+region = { type = "block", min = [-0.001, -0.001, 0.007], max = [0.001, 0.001, 0.008] }
+
+# Floor wall (z = 0, normal +z)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "glass"
+
+[output]
+dir = "examples/bench_hertz_rebound"
+
+[run]
+steps = 20000             # Enough steps for impact and rebound (no gravity, particle close to wall)
+thermo = 1000

examples/bench_hertz_rebound/main.rs

@@ -0,0 +1,164 @@
+//! Hertz contact rebound benchmark — validates coefficient of restitution,
+//! contact duration, and peak overlap against Hertz contact theory.
+//!
+//! Drops a single sphere onto a rigid flat wall and records the impact
+//! velocity, rebound velocity, contact duration, and peak overlap.
+//!
+//! ```bash
+//! cargo run --release --example bench_hertz_rebound --no-default-features -- examples/bench_hertz_rebound/config.toml
+//! ```
+
+use mddem::prelude::*;
+use mddem::dem_atom::DemAtom;
+use std::fs;
+use std::io::Write as IoWrite;
+
+/// Tracks contact state for the rebound measurement.
+struct ReboundTracker {
+    /// True if the particle has been in contact with the wall.
+    was_in_contact: bool,
+    /// True once the particle has separated after contact.
+    finished: bool,
+    /// z-velocity just before first contact (impact velocity, negative = downward).
+    v_impact: f64,
+    /// z-velocity just after separation (rebound velocity, positive = upward).
+    v_rebound: f64,
+    /// Timestep when contact first occurs.
+    step_contact_start: usize,
+    /// Timestep when contact ends (separation).
+    step_contact_end: usize,
+    /// Maximum overlap during contact.
+    max_overlap: f64,
+    /// z-velocity at the previous step (to capture pre-contact velocity).
+    prev_vz: f64,
+    /// Output directory for results.
+    output_dir: String,
+}
+
+impl ReboundTracker {
+    fn new() -> Self {
+        Self {
+            was_in_contact: false,
+            finished: false,
+            v_impact: 0.0,
+            v_rebound: 0.0,
+            step_contact_start: 0,
+            step_contact_end: 0,
+            max_overlap: 0.0,
+            prev_vz: 0.0,
+            output_dir: String::new(),
+        }
+    }
+}
+
+fn main() {
+    let mut app = App::new();
+    app.add_plugins(CorePlugins)
+        .add_plugins(GranularDefaultPlugins)
+        .add_plugins(WallPlugin);
+
+    app.add_resource(ReboundTracker::new());
+
+    app.add_update_system(track_rebound, ScheduleSet::PostFinalIntegration);
+
+    app.start();
+}
+
+/// System that monitors the single particle's contact with the floor wall
+/// and records impact/rebound velocity, contact duration, and peak overlap.
+fn track_rebound(
+    atoms: Res<Atom>,
+    registry: Res<AtomDataRegistry>,
+    run_state: Res<RunState>,
+    input: Res<Input>,
+    mut tracker: ResMut<ReboundTracker>,
+) {
+    if tracker.finished || atoms.nlocal == 0 {
+        return;
+    }
+
+    let dem = registry.expect::<DemAtom>("track_rebound");
+    let step = run_state.total_cycle;
+
+    // Single particle at index 0
+    let z = atoms.pos[0][2];
+    let vz = atoms.vel[0][2];
+    let r = dem.radius[0];
+
+    // Check overlap with the floor wall (z = 0, normal +z)
+    // The floor is the first wall defined in config
+    // Overlap = radius - z (positive when overlapping)
+    let overlap = r - z;
+    let in_contact = overlap > 0.0;
+
+    if !tracker.was_in_contact && !in_contact {
+        // Pre-contact: record velocity for next step
+        tracker.prev_vz = vz;
+    } else if !tracker.was_in_contact && in_contact {
+        // First contact! Record impact velocity from previous step
+        tracker.was_in_contact = true;
+        tracker.v_impact = tracker.prev_vz;
+        tracker.step_contact_start = step;
+        tracker.max_overlap = overlap;
+    } else if tracker.was_in_contact && in_contact {
+        // During contact: track max overlap
+        if overlap > tracker.max_overlap {
+            tracker.max_overlap = overlap;
+        }
+    } else if tracker.was_in_contact && !in_contact {
+        // Separation: contact has ended
+        tracker.finished = true;
+        tracker.v_rebound = vz;
+        tracker.step_contact_end = step;
+
+        let dt = atoms.dt;
+        let contact_steps = tracker.step_contact_end - tracker.step_contact_start;
+        let contact_time = contact_steps as f64 * dt;
+        let cor_measured = (tracker.v_rebound / tracker.v_impact).abs();
+
+        // Determine output directory
+        let out_dir = if let Some(ref dir) = input.output_dir {
+            dir.clone()
+        } else {
+            "examples/bench_hertz_rebound/data".to_string()
+        };
+        tracker.output_dir = out_dir.clone();
+
+        // Ensure data directory exists
+        let data_dir = format!("{}/data", out_dir);
+        fs::create_dir_all(&data_dir).ok();
+
+        // Write results to file
+        let results_file = format!("{}/data/rebound_results.csv", out_dir);
+        let mut f = fs::File::create(&results_file)
+            .unwrap_or_else(|e| panic!("Cannot create {}: {}", results_file, e));
+        writeln!(f, "v_impact,v_rebound,cor_measured,contact_time,max_overlap,dt,radius,density")
+            .unwrap();
+        writeln!(
+            f,
+            "{:.10e},{:.10e},{:.10e},{:.10e},{:.10e},{:.10e},{:.10e},{:.10e}",
+            tracker.v_impact.abs(),
+            tracker.v_rebound.abs(),
+            cor_measured,
+            contact_time,
+            tracker.max_overlap,
+            dt,
+            r,
+            dem.density[0],
+        )
+        .unwrap();
+
+        println!("=== Hertz Rebound Results ===");
+        println!("  Impact velocity:  {:.6e} m/s", tracker.v_impact.abs());
+        println!("  Rebound velocity: {:.6e} m/s", tracker.v_rebound.abs());
+        println!("  COR (measured):   {:.6}", cor_measured);
+        println!("  Contact duration: {:.6e} s ({} steps)", contact_time, contact_steps);
+        println!("  Peak overlap:     {:.6e} m", tracker.max_overlap);
+        println!("  Timestep dt:      {:.6e} s", dt);
+        println!("  Results saved to: {}", results_file);
+    }
+
+    if !tracker.finished && !in_contact {
+        tracker.prev_vz = vz;
+    }
+}