Add lunar regolith JKR adhesion angle of repose benchmark

Cargo.toml

@@ -126,6 +126,10 @@ path = "examples/example_template/main.rs"
 name = "cone_hopper"
 path = "examples/cone_hopper/main.rs"
 
+[[example]]
+name = "bench_lunar_regolith"
+path = "examples/bench_lunar_regolith/main.rs"
+
 [profile.release]
 opt-level = 3
 lto = "fat"

examples/bench_lunar_regolith/README.md

@@ -0,0 +1,98 @@
+# Lunar Regolith Cohesive Angle of Repose Benchmark
+
+Validates JKR adhesion in reduced gravity by simulating a funnel-pour of cohesive particles and measuring the resulting angle of repose.
+
+## Physics
+
+Cohesive particles form steeper piles than non-cohesive ones. The key dimensionless parameter is the **granular Bond number**:
+
+$$\text{Bo} = \frac{F_{\text{adhesion}}}{mg} = \frac{\frac{3}{2}\pi\gamma R^*}{mg}$$
+
+where $\gamma$ is JKR surface energy, $R^* = R/2$ for equal spheres, and $g$ is gravitational acceleration.
+
+| Regime | Bond Number | Behavior |
+|--------|-------------|----------|
+| Gravity-dominated | Bo << 1 | Angle ≈ internal friction angle (~25-30°) |
+| Transitional | Bo ~ 1 | Angle increases significantly |
+| Cohesion-dominated | Bo >> 1 | Very steep piles, near-vertical cliffs |
+
+Lunar gravity (1.62 m/s²) gives ~6× higher Bo than Earth gravity (9.81 m/s²) for the same particles, explaining why lunar regolith can maintain steeper slopes.
+
+## Setup
+
+- **Geometry**: Quasi-2D (periodic in y, 3 particle diameters thick)
+- **Particles**: 350 spheres, R = 1 mm, ρ = 1500 kg/m³ (scaled up from real regolith for tractability)
+- **Contact model**: Hertz-Mindlin with JKR adhesion
+- **Material**: E = 5 MPa (soft for fast simulation), μ = 0.5, μ_r = 0.1
+- **Insertion**: Rate-based pouring from narrow slot above center (10 particles every 500 steps)
+- **Floor**: Flat wall at z = 0
+- **Boundaries**: shrink-wrap in z, periodic in y, fixed in x (wide enough pile doesn't reach walls)
+
+## Angle Measurement
+
+The pile angle is measured by:
+1. Binning particles by |x| (exploiting symmetry about x=0)
+2. Finding the surface height in each bin (max z + radius)
+3. Fitting a line to the surface profile
+4. Angle = atan(|slope|)
+
+## Running
+
+### Single case (default: lunar gravity, medium adhesion)
+
+```bash
+cargo run --release --no-default-features --example bench_lunar_regolith -- examples/bench_lunar_regolith/config.toml
+```
+
+### Full parametric study (8 cases: {Earth, Moon} × {0, 5, 20, 50 mJ/m²})
+
+```bash
+python examples/bench_lunar_regolith/run_benchmark.py
+```
+
+### Validation
+
+```bash
+python examples/bench_lunar_regolith/validate.py
+```
+
+### Plots
+
+```bash
+python examples/bench_lunar_regolith/plot.py
+```
+
+## Validation Checks
+
+1. **Non-cohesive angle**: 15-45° (consistent with μ = 0.5)
+2. **Angle vs adhesion**: Monotonically increases with surface energy
+3. **Lunar vs Earth**: Steeper piles on Moon for same adhesion (higher Bo)
+4. **Significant cohesion effect**: High-adhesion angles > no-adhesion + 3°
+
+## Expected Results
+
+| Gravity | γ [mJ/m²] | Bo | Expected Angle |
+|---------|-----------|-----|---------------|
+| Earth | 0 | 0 | ~25-30° |
+| Earth | 5 | 0.19 | ~27-32° |
+| Earth | 20 | 0.76 | ~32-40° |
+| Earth | 50 | 1.91 | ~40-55° |
+| Moon | 0 | 0 | ~25-30° |
+| Moon | 5 | 1.16 | ~35-45° |
+| Moon | 20 | 4.63 | ~50-65° |
+| Moon | 50 | 11.6 | ~65-80° |
+
+## Output
+
+- `results.csv` — Tabulated results from parametric study
+- `angle_vs_surface_energy.png` — Angle vs γ for Earth and Moon
+- `angle_vs_bond_number.png` — Universal scaling with Bond number
+
+![Angle vs Surface Energy](angle_vs_surface_energy.png)
+![Angle vs Bond Number](angle_vs_bond_number.png)
+
+## References
+
+- Castellanos, A. (2005). "The relationship between attractive interparticle forces and bulk behaviour in dry and uncharged fine powders." *Advances in Physics*, 54(4), 263-376.
+- Johnson, K.L., Kendall, K., Roberts, A.D. (1971). "Surface energy and the contact of elastic solids." *Proc. R. Soc. Lond. A*, 324, 301-313.
+- Carrier, W.D. (2003). "Particle size distribution of lunar soil." *Journal of Geotechnical and Geoenvironmental Engineering*, 129(10), 956-959.

examples/bench_lunar_regolith/config.toml

@@ -0,0 +1,77 @@
+# Lunar Regolith Cohesive Angle of Repose — Default Config
+#
+# Funnel-pour test: particles are dropped from a narrow slot above the center
+# of a flat floor. They accumulate into a pile whose angle depends on
+# friction, JKR adhesion, and gravity.
+#
+# This is the "demo" config. Use run_benchmark.py for the full parametric study.
+
+[comm]
+processors_x = 1
+processors_y = 1
+processors_z = 1
+
+[domain]
+x_low = -0.06           # [m] Wide enough that pile doesn't reach walls
+x_high = 0.06
+y_low = 0.0             # [m] Quasi-2D: periodic in y
+y_high = 0.006          # 3 particle diameters thick
+z_low = 0.0             # [m] Floor level
+z_high = 0.10           # [m] Tall enough for insertion region
+periodic_x = false
+periodic_y = true       # Quasi-2D slice
+periodic_z = false
+boundary_z = "shrink-wrap"  # Track highest particle
+
+[neighbor]
+skin_fraction = 1.2     # Extra skin for JKR extended range
+bin_size = 0.004        # [m] ~2x particle diameter
+every = 5
+rebuild_on_pbc_wrap = true
+
+[gravity]
+gx = 0.0
+gy = 0.0
+gz = -1.62              # [m/s^2] Lunar surface gravity
+
+[dem]
+contact_model = "hertz"
+
+[[dem.materials]]
+name = "regolith"
+youngs_mod = 5e6         # [Pa] Soft for computational tractability
+poisson_ratio = 0.3
+restitution = 0.3        # Low restitution for fast energy dissipation
+friction = 0.5           # Typical for angular lunar regolith
+rolling_friction = 0.1   # Resists rolling for realistic pile formation
+surface_energy = 0.02    # [J/m^2] JKR adhesion — moderate cohesion
+
+# Rate-based insertion: drop particles from narrow slot above center
+# Slot: x [-0.005, 0.005], z [0.04, 0.06]
+# Insert 10 particles every 500 steps, up to 350 total
+[[particles.insert]]
+material = "regolith"
+radius = 0.001           # [m] 1 mm
+density = 1500.0         # [kg/m^3]
+rate = 10                # Particles per insertion event
+rate_interval = 500      # Insert every 500 timesteps
+rate_limit = 350         # Total particles to insert
+velocity_z = -0.5        # [m/s] Initial downward velocity
+region = { type = "block", min = [-0.005, 0.0, 0.04], max = [0.005, 0.006, 0.06] }
+
+# Floor wall
+[[wall]]
+point_z = 0.0
+normal_z = 1.0
+material = "regolith"
+
+[output]
+dir = "examples/bench_lunar_regolith/output"
+
+# Single stage: pour and settle
+[[run]]
+name = "pour"
+steps = 200000           # 1.0s at dt=5e-6 — enough to pour and settle
+thermo = 10000
+dt = 5e-6                # [s] Stable timestep for E=5e6
+dump_interval = 200000   # Only dump at end

examples/bench_lunar_regolith/main.rs

@@ -0,0 +1,24 @@
+//! Lunar Regolith Cohesive Angle of Repose Benchmark
+//!
+//! Validates JKR adhesion in reduced gravity by simulating a funnel-pour
+//! of cohesive particles. Particles are inserted from a narrow slot above
+//! the center of a flat floor and naturally form a conical/triangular pile.
+//! The pile angle depends on friction, cohesion, and gravity.
+//!
+//! Run with:
+//! ```bash
+//! cargo run --release --no-default-features --example bench_lunar_regolith \
+//!   -- examples/bench_lunar_regolith/config.toml
+//! ```
+
+use mddem::prelude::*;
+
+fn main() {
+    let mut app = App::new();
+    app.add_plugins(CorePlugins)
+        .add_plugins(GranularDefaultPlugins)
+        .add_plugins(GravityPlugin)
+        .add_plugins(WallPlugin);
+
+    app.start();
+}

examples/bench_lunar_regolith/plot.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+"""
+Generate publication-quality plots for the lunar regolith benchmark.
+
+Reads results.csv and produces:
+  1. angle_vs_surface_energy.png — Angle of repose vs surface energy (Earth & Lunar)
+  2. angle_vs_bond_number.png — Angle vs Bond number (universal scaling)
+
+Usage:
+    python examples/bench_lunar_regolith/plot.py
+"""
+
+import os
+import sys
+import numpy as np
+
+try:
+    import matplotlib
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+except ImportError:
+    print("ERROR: matplotlib is required. Install with: pip install matplotlib")
+    sys.exit(1)
+
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+RESULTS_FILE = os.path.join(SCRIPT_DIR, "results.csv")
+DATA_RESULTS = os.path.join(SCRIPT_DIR, "data", "results.csv")
+
+
+def load_results(path):
+    """Load results CSV into arrays."""
+    labels, gravities, gammas, bond_numbers, angles = [], [], [], [], []
+    with open(path) as f:
+        f.readline()  # skip header
+        for line in f:
+            parts = line.strip().split(",")
+            if len(parts) < 5:
+                continue
+            try:
+                angle = float(parts[4])
+            except ValueError:
+                continue
+            if np.isnan(angle):
+                continue
+            labels.append(parts[0])
+            gravities.append(float(parts[1]))
+            gammas.append(float(parts[2]))
+            bond_numbers.append(float(parts[3]))
+            angles.append(angle)
+    return (np.array(labels), np.array(gravities), np.array(gammas),
+            np.array(bond_numbers), np.array(angles))
+
+
+def castellanos_trend(bo, theta0=27.0):
+    """
+    Empirical cohesive angle of repose trend (inspired by Castellanos 2005).
+
+    theta = theta0 + k * Bo / (1 + Bo)
+    where theta0 is the non-cohesive angle and k sets the max increase.
+    Captures: Bo<<1 -> theta0, Bo~1 -> intermediate, Bo>>1 -> plateau.
+    """
+    k = 50.0  # Maximum additional angle [deg]
+    return theta0 + k * bo / (1.0 + bo)
+
+
+def main():
+    # Find results
+    if os.path.isfile(RESULTS_FILE):
+        results_path = RESULTS_FILE
+    elif os.path.isfile(DATA_RESULTS):
+        results_path = DATA_RESULTS
+    else:
+        print(f"ERROR: No results file found at {RESULTS_FILE}")
+        print("Run: python run_benchmark.py first.")
+        sys.exit(1)
+
+    labels, gravities, gammas, bond_numbers, angles = load_results(results_path)
+
+    if len(labels) == 0:
+        print("ERROR: No valid results to plot.")
+        sys.exit(1)
+
+    # --- Plot 1: Angle vs Surface Energy ---
+    fig, ax = plt.subplots(figsize=(8, 5.5))
+
+    earth_mask = gravities < -5.0
+    lunar_mask = gravities > -5.0
+
+    if np.any(earth_mask):
+        idx = np.argsort(gammas[earth_mask])
+        ax.plot(gammas[earth_mask][idx] * 1000, angles[earth_mask][idx],
+                "s-", color="#2196F3", markersize=8, linewidth=2,
+                label="Earth (g = 9.81 m/s$^2$)", markerfacecolor="white",
+                markeredgewidth=2)
+
+    if np.any(lunar_mask):
+        idx = np.argsort(gammas[lunar_mask])
+        ax.plot(gammas[lunar_mask][idx] * 1000, angles[lunar_mask][idx],
+                "o-", color="#FF5722", markersize=8, linewidth=2,
+                label="Moon (g = 1.62 m/s$^2$)", markerfacecolor="white",
+                markeredgewidth=2)
+
+    ax.set_xlabel("Surface Energy, $\\gamma$ [mJ/m$^2$]", fontsize=13)
+    ax.set_ylabel("Angle of Repose [deg]", fontsize=13)
+    ax.set_title("Cohesive Angle of Repose: Earth vs Moon", fontsize=14)
+    ax.legend(fontsize=11, framealpha=0.9)
+    ax.set_ylim(0, 90)
+    ax.set_xlim(-2, 55)
+    ax.grid(True, alpha=0.3)
+    ax.tick_params(labelsize=11)
+
+    # Annotate Bond numbers
+    for i in range(len(labels)):
+        if bond_numbers[i] > 0.01:
+            offset_y = 3 if gravities[i] < -5 else -5
+            ax.annotate(f"Bo={bond_numbers[i]:.2f}",
+                        (gammas[i] * 1000, angles[i]),
+                        textcoords="offset points", xytext=(5, offset_y),
+                        fontsize=8, color="gray")
+
+    fig.tight_layout()
+    out1 = os.path.join(SCRIPT_DIR, "angle_vs_surface_energy.png")
+    fig.savefig(out1, dpi=150)
+    print(f"Saved: {out1}")
+    plt.close(fig)
+
+    # --- Plot 2: Angle vs Bond Number (universal scaling) ---
+    fig, ax = plt.subplots(figsize=(8, 5.5))
+
+    if np.any(earth_mask):
+        ax.scatter(bond_numbers[earth_mask], angles[earth_mask],
+                   s=100, marker="s", facecolors="white", edgecolors="#2196F3",
+                   linewidths=2, zorder=5, label="DEM — Earth")
+
+    if np.any(lunar_mask):
+        ax.scatter(bond_numbers[lunar_mask], angles[lunar_mask],
+                   s=100, marker="o", facecolors="white", edgecolors="#FF5722",
+                   linewidths=2, zorder=5, label="DEM — Moon")
+
+    # Theoretical / empirical curve
+    bo_theory = np.linspace(0, 15, 200)
+    theta_theory = castellanos_trend(bo_theory)
+    ax.plot(bo_theory, theta_theory, "k--", linewidth=1.5, alpha=0.7,
+            label="Empirical trend (Castellanos 2005)")
+
+    ax.set_xlabel("Bond Number, Bo = $F_{adhesion}$ / ($mg$)", fontsize=13)
+    ax.set_ylabel("Angle of Repose [deg]", fontsize=13)
+    ax.set_title("Angle of Repose vs Granular Bond Number", fontsize=14)
+    ax.legend(fontsize=11, framealpha=0.9)
+    ax.set_ylim(0, 90)
+    ax.set_xlim(-0.5, max(15, np.max(bond_numbers) * 1.2) if len(bond_numbers) > 0 else 15)
+    ax.grid(True, alpha=0.3)
+    ax.tick_params(labelsize=11)
+
+    # Regime annotations
+    ax.axvspan(-0.5, 1.0, alpha=0.05, color="blue")
+    ax.axvspan(1.0, 5.0, alpha=0.05, color="orange")
+    ax.axvspan(5.0, 20.0, alpha=0.05, color="red")
+    ax.text(0.3, 85, "Gravity\ndominated", fontsize=9, ha="center", color="blue", alpha=0.6)
+    ax.text(2.5, 85, "Transitional", fontsize=9, ha="center", color="orange", alpha=0.6)
+    ax.text(10.0, 85, "Cohesion\ndominated", fontsize=9, ha="center", color="red", alpha=0.6)
+
+    fig.tight_layout()
+    out2 = os.path.join(SCRIPT_DIR, "angle_vs_bond_number.png")
+    fig.savefig(out2, dpi=150)
+    print(f"Saved: {out2}")
+    plt.close(fig)
+
+    print("\nAll plots generated successfully.")
+
+
+if __name__ == "__main__":
+    main()