diff --git a/Cargo.toml b/Cargo.toml index 93467f5..e9f5ac7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -126,6 +126,10 @@ path = "examples/example_template/main.rs" name = "cone_hopper" path = "examples/cone_hopper/main.rs" +[[example]] +name = "bench_multisphere_segregation" +path = "examples/bench_multisphere_segregation/main.rs" + [profile.release] opt-level = 3 lto = "fat" diff --git a/crates/dem_clump/README.md b/crates/dem_clump/README.md index 5a91151..a7033f1 100644 --- a/crates/dem_clump/README.md +++ b/crates/dem_clump/README.md @@ -23,10 +23,10 @@ A **clump** is a rigid body composed of multiple overlapping spheres. Each spher ## TOML Configuration -Define clump types in `[[dem.clumps]]`: +Define clump types in `[[clump.definitions]]`: ```toml -[[dem.clumps]] +[[clump.definitions]] name = "dimer" spheres = [ { offset = [-0.0003, 0.0, 0.0], radius = 0.001 }, @@ -58,7 +58,7 @@ app.run(); ``` The plugin automatically: -- Loads clump definitions from `dem.clumps` config +- Loads clump definitions from `clump.definitions` config - Registers per-atom clump data - Aggregates forces from sub-spheres to parents - Updates sub-sphere positions/velocities after parent integration diff --git a/crates/dem_clump/src/lib.rs b/crates/dem_clump/src/lib.rs index 144f03f..c4c560a 100644 --- a/crates/dem_clump/src/lib.rs +++ b/crates/dem_clump/src/lib.rs @@ -8,7 +8,7 @@ //! # Configuration //! //! ```toml -//! [[dem.clumps]] +//! [[clump.definitions]] //! name = "dimer" //! spheres = [ //! { offset = [-0.0003, 0.0, 0.0], radius = 0.001 }, @@ -49,7 +49,7 @@ pub struct ClumpSphereConfig { pub radius: f64, } -/// A clump type definition from `[[dem.clumps]]`. +/// A clump type definition from `[[clump.definitions]]`. #[derive(Deserialize, Clone, Debug)] #[serde(deny_unknown_fields)] pub struct ClumpDef { @@ -59,17 +59,24 @@ pub struct ClumpDef { pub spheres: Vec, } -/// Extended DEM config with optional clump definitions. -/// Uses `#[serde(flatten)]` with a HashMap to absorb unknown `[dem]` keys -/// (materials, contact_model, etc.) without failing deserialization. +/// Top-level `[clump]` configuration section with clump type definitions. +/// +/// Loaded from the `[clump]` TOML section (separate from `[dem]`) so that +/// `DemConfig`'s `deny_unknown_fields` doesn't reject the `clumps` key. +/// +/// ```toml +/// [[clump.definitions]] +/// name = "dimer" +/// spheres = [ +/// { offset = [-0.0003, 0.0, 0.0], radius = 0.001 }, +/// { offset = [0.0003, 0.0, 0.0], radius = 0.001 }, +/// ] +/// ``` #[derive(Deserialize, Clone, Default)] -pub struct DemClumpConfig { +pub struct ClumpPluginConfig { /// Clump type definitions. #[serde(default)] - pub clumps: Option>, - /// Pass-through fields we don't care about. - #[serde(flatten)] - pub _rest: std::collections::HashMap, + pub definitions: Option>, } // ── Per-atom clump data ───────────────────────────────────────────────────── @@ -238,9 +245,9 @@ impl Plugin for ClumpPlugin { // Load clump definitions from config let mut registry = ClumpRegistry::new(); - // Try to load dem.clumps from config - let dem_clump_config = Config::load::(app, "dem"); - if let Some(clumps) = dem_clump_config.clumps { + // Load clump definitions from [clump] section + let clump_config = Config::load::(app, "clump"); + if let Some(clumps) = clump_config.definitions { for def in clumps { assert!( !def.spheres.is_empty(), diff --git a/examples/bench_multisphere_segregation/README.md b/examples/bench_multisphere_segregation/README.md new file mode 100644 index 0000000..62ffdba --- /dev/null +++ b/examples/bench_multisphere_segregation/README.md @@ -0,0 +1,77 @@ +# Multisphere Segregation Benchmark + +## Physics + +This benchmark validates the **clump/multisphere** implementation by simulating +size-driven segregation of single spheres vs dimer clumps under vertical +vibration — the **Brazil nut effect** for particle shape. + +A 50/50 mixture of: +- **75 single spheres** (type 0, radius 1.0 mm) +- **75 dimer clumps** (type 1, two overlapping spheres with 1.0 mm offset) + +is placed in a box with: +- Periodic lateral boundaries (x, y) +- A vertically oscillating floor (sinusoidal, z-direction) +- A static ceiling + +### Vibration Parameters + +| Parameter | Value | +|-----------|-------| +| Amplitude A | 2.0 mm | +| Frequency f | 15 Hz | +| Angular frequency ω | 94.25 rad/s | +| Dimensionless acceleration Γ = Aω²/g | **3.6** | + +### Expected Behavior + +Under vibration at Γ ≈ 3.6, the dimers have a larger effective size than single +spheres. The **Brazil nut effect** (also called granular convection / geometric +segregation) causes larger effective particles to rise to the top of the bed. + +The **segregation index** is defined as: + +$$S = \frac{z_{\text{dimer}} - z_{\text{sphere}}}{z_{\text{dimer}} + z_{\text{sphere}}}$$ + +where $z$ denotes the mass-weighted center-of-mass height. **S > 0** indicates +dimers are preferentially at the top. + +### Assumptions / Simplifications + +- 3D simulation with periodic lateral boundaries (quasi-infinite horizontal extent) +- Monodisperse spheres and identical dimers (no size distribution) +- Softened Young's modulus (5×10⁷ Pa) for computational efficiency +- Hertz–Mindlin contact with rolling friction +- No cohesion or adhesion + +## Running + +```bash +# Build and run (release mode, ~2-3 min) +cargo run --release --example bench_multisphere_segregation --no-default-features \ + -- examples/bench_multisphere_segregation/config.toml + +# Validate results (PASS/FAIL) +python3 examples/bench_multisphere_segregation/validate.py + +# Generate plots +python3 examples/bench_multisphere_segregation/plot.py +``` + +## Output + +- `data/segregation.csv` — step, time, z_sphere, z_dimer, segregation_index +- `data/Thermo.txt` — standard thermo output +- `segregation_index.png` — segregation index vs time +- `com_trajectories.png` — COM height trajectories + +## Validation Criteria + +| Check | Criterion | +|-------|-----------| +| No NaN/Inf | All data must be finite | +| Physical z | Both COM heights must be positive | +| S > 0 at steady state | Mean S in final 20% of run must be positive | +| Significant segregation | Mean final S > 0.005 | +| Positive trend | S increases (or is already high) in second half | diff --git a/examples/bench_multisphere_segregation/config.toml b/examples/bench_multisphere_segregation/config.toml new file mode 100644 index 0000000..b368ec2 --- /dev/null +++ b/examples/bench_multisphere_segregation/config.toml @@ -0,0 +1,104 @@ +# Multisphere Segregation Benchmark +# ================================== +# A mixture of single spheres and dimer clumps in a vertically shaken box. +# Under vibration at Gamma ≈ 3.6, dimers (larger effective size) segregate +# upward due to the Brazil nut effect. +# +# Run with: +# cargo run --release --example bench_multisphere_segregation --no-default-features \ +# -- examples/bench_multisphere_segregation/config.toml + +[comm] +processors_x = 1 +processors_y = 1 +processors_z = 1 + +[domain] +# Box: 20mm x 20mm x 60mm (tall to allow segregation) +x_low = 0.0 +x_high = 0.02 # [m] +y_low = 0.0 +y_high = 0.02 # [m] +z_low = 0.0 +z_high = 0.06 # [m] +periodic_x = true # Periodic lateral boundaries +periodic_y = true +periodic_z = false + +[neighbor] +skin_fraction = 1.1 # [dimensionless] +bin_size = 0.004 # [m] — slightly larger than 2x particle diameter +rebuild_on_pbc_wrap = true + +[gravity] +gz = -9.81 # Gravitational acceleration [m/s²] + +# --- Material properties: glass-like --- +[[dem.materials]] +name = "glass" +youngs_mod = 5.0e7 # [Pa] — softened for faster timestep +poisson_ratio = 0.3 # [dimensionless] +restitution = 0.5 # [dimensionless] — moderate dissipation +friction = 0.5 # Sliding friction [dimensionless] +rolling_friction = 0.1 # Rolling friction [dimensionless] + +# --- Dimer clump definition --- +# Two overlapping spheres offset along x, each with radius 1.0 mm. +# Center-to-center distance 1.0 mm → aspect ratio ≈ 1.5 +# Note: clump definitions live under [clump], not [dem], because DemConfig +# uses deny_unknown_fields. +[[clump.definitions]] +name = "dimer" +spheres = [ + { offset = [-0.0005, 0.0, 0.0], radius = 0.001 }, + { offset = [ 0.0005, 0.0, 0.0], radius = 0.001 }, +] + +# --- Particle insertion --- +# Type 0: single spheres (75 particles, r = 1.0 mm) +[[particles.insert]] +material = "glass" +count = 75 +radius = 0.001 # [m] +density = 2500.0 # [kg/m³] +region = { type = "block", min = [0.0, 0.0, 0.005], max = [0.02, 0.02, 0.04] } + +# Type 1: dimer clumps (75 dimers) +[[particles.insert]] +material = "glass" +clump = "dimer" +count = 75 +density = 2500.0 # [kg/m³] +region = { type = "block", min = [0.0, 0.0, 0.005], max = [0.02, 0.02, 0.04] } + +# --- Walls --- +# Floor wall — oscillates sinusoidally along z +# Oscillation: A = 0.002 m, f = 15 Hz → ω = 2π·15 ≈ 94.25 rad/s +# Γ = A·ω²/g = 0.002 × 94.25² / 9.81 ≈ 3.6 +[[wall]] +type = "plane" +point_z = 0.0 +normal_z = 1.0 +material = "glass" +name = "floor" +oscillate = { amplitude = 0.002, frequency = 15.0 } + +# Ceiling wall — static, high above particles +[[wall]] +type = "plane" +point_z = 0.06 +normal_z = -1.0 +material = "glass" + +[output] +dir = "examples/bench_multisphere_segregation/data" + +[vtp] +interval = 10000 # Write VTP files every 10000 steps for visualization + +# Single run: particles settle under gravity while floor oscillates. +# The initial transient includes both settling and the onset of segregation. +# Total runtime at ~1e-6 s timestep: ~0.5 s of physical time. +[run] +steps = 500000 +thermo = 2000 diff --git a/examples/bench_multisphere_segregation/main.rs b/examples/bench_multisphere_segregation/main.rs new file mode 100644 index 0000000..e661ec2 --- /dev/null +++ b/examples/bench_multisphere_segregation/main.rs @@ -0,0 +1,148 @@ +//! Multisphere segregation benchmark — validates that dimer clumps rise above +//! single spheres under vertical vibration (Brazil nut effect for shape). +//! +//! A 50/50 mixture of single spheres (type 0) and dimer clumps (type 1) is +//! placed in a box with a vertically oscillating floor, periodic lateral walls, +//! and a static ceiling. Under vibration at Γ = A·ω²/g ≈ 3.6, the dimers' +//! larger effective size causes them to segregate upward. +//! +//! A custom system writes `segregation.csv` every N steps with the center-of-mass +//! heights of spheres vs dimers and the segregation index S. +//! +//! ```bash +//! cargo run --release --example bench_multisphere_segregation --no-default-features \ +//! -- examples/bench_multisphere_segregation/config.toml +//! ``` + +use std::fs; +use std::io::Write; + +use mddem::prelude::*; + +fn main() { + let mut app = App::new(); + app.add_plugins(CorePlugins) + .add_plugins(GranularDefaultPlugins) + .add_plugins(GravityPlugin) + .add_plugins(WallPlugin) + .add_plugins(ClumpPlugin); + + // Measure segregation periodically + app.add_update_system(measure_segregation, ScheduleSet::PostFinalIntegration); + + app.start(); +} + +/// Measure and write segregation data: average z of spheres vs dimers. +/// +/// Single spheres have type 0, clump parents (dimers) have type 1. +/// Sub-spheres are excluded (they are not independent particles). +fn measure_segregation( + atoms: Res, + registry: Res, + run_state: Res, + comm: Res, + input: Res, +) { + let step = run_state.total_cycle; + + // Write every 500 steps + if step % 500 != 0 { + return; + } + + let clump = registry.get::(); + let clump = match clump { + Some(c) => c, + None => return, + }; + + let nlocal = atoms.nlocal as usize; + + // Accumulate z-positions weighted by mass for each species + let mut sphere_z_sum = 0.0_f64; + let mut sphere_mass_sum = 0.0_f64; + let mut dimer_z_sum = 0.0_f64; + let mut dimer_mass_sum = 0.0_f64; + + for i in 0..nlocal { + // Skip sub-spheres (they are part of clumps, not independent) + if i < clump.clump_id.len() && clump.clump_id[i] > 0.0 && clump.is_parent_flag[i] < 0.5 { + continue; + } + + let z = atoms.pos[i][2]; + let m = atoms.mass[i]; + + if atoms.atom_type[i] == 0 { + // Single sphere + sphere_z_sum += m * z; + sphere_mass_sum += m; + } else if atoms.atom_type[i] == 1 { + // Clump parent (dimer) + if i < clump.clump_id.len() && clump.is_parent_flag[i] > 0.5 { + dimer_z_sum += m * z; + dimer_mass_sum += m; + } + } + } + + // Global reduction + let sphere_z_total = comm.all_reduce_sum_f64(sphere_z_sum); + let sphere_mass_total = comm.all_reduce_sum_f64(sphere_mass_sum); + let dimer_z_total = comm.all_reduce_sum_f64(dimer_z_sum); + let dimer_mass_total = comm.all_reduce_sum_f64(dimer_mass_sum); + + if comm.rank() != 0 { + return; + } + + let z_sphere = if sphere_mass_total > 0.0 { + sphere_z_total / sphere_mass_total + } else { + 0.0 + }; + let z_dimer = if dimer_mass_total > 0.0 { + dimer_z_total / dimer_mass_total + } else { + 0.0 + }; + + // Segregation index: S = (z_dimer - z_sphere) / (z_dimer + z_sphere) + // S > 0 means dimers are preferentially at the top + let seg_index = if (z_dimer + z_sphere) > 0.0 { + (z_dimer - z_sphere) / (z_dimer + z_sphere) + } else { + 0.0 + }; + + // Compute time from step and timestep + let time = step as f64 * atoms.dt; + + // Output directory from Input resource + let output_dir = match input.output_dir.as_deref() { + Some(dir) => dir.to_string(), + None => "data".to_string(), + }; + let filepath = format!("{}/segregation.csv", output_dir); + + // On first call, write header + use std::sync::Once; + static INIT: Once = Once::new(); + INIT.call_once(|| { + fs::create_dir_all(&output_dir).ok(); + let mut f = fs::File::create(&filepath).expect("Cannot create segregation.csv"); + writeln!(f, "step,time,z_sphere,z_dimer,segregation_index").unwrap(); + }); + + let mut f = fs::OpenOptions::new() + .append(true) + .open(&filepath) + .expect("Cannot open segregation.csv"); + writeln!( + f, + "{},{:.8e},{:.8e},{:.8e},{:.8e}", + step, time, z_sphere, z_dimer, seg_index + ) + .unwrap(); +} diff --git a/examples/bench_multisphere_segregation/plot.py b/examples/bench_multisphere_segregation/plot.py new file mode 100644 index 0000000..9a997d2 --- /dev/null +++ b/examples/bench_multisphere_segregation/plot.py @@ -0,0 +1,78 @@ +#!/usr/bin/env python3 +""" +Plot multisphere segregation benchmark results. + +Generates publication-quality figures: + 1. Segregation index S vs time + 2. Center-of-mass height trajectories for spheres vs dimers + +Reads: data/segregation.csv +Saves: segregation_index.png, com_trajectories.png +""" + +import os +import sys +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +script_dir = os.path.dirname(os.path.abspath(__file__)) +data_file = os.path.join(script_dir, "data", "segregation.csv") + +if not os.path.isfile(data_file): + print(f"ERROR: {data_file} not found. Run simulation first.") + sys.exit(1) + +data = np.loadtxt(data_file, delimiter=",", skiprows=1) +if data.ndim == 1: + data = data.reshape(1, -1) + +steps = data[:, 0] +times = data[:, 1] * 1000 # Convert to ms +z_sphere = data[:, 2] * 1000 # Convert to mm +z_dimer = data[:, 3] * 1000 # Convert to mm +seg_index = data[:, 4] + +# ── Figure 1: Segregation Index vs Time ────────────────────────────────── + +fig, ax = plt.subplots(figsize=(8, 5)) +ax.plot(times, seg_index, "b-", linewidth=1.5, label="Segregation index $S$") +ax.axhline(y=0, color="k", linestyle="--", linewidth=0.8, alpha=0.5) + +# Add running average +if len(seg_index) > 20: + window = max(len(seg_index) // 20, 5) + running_avg = np.convolve(seg_index, np.ones(window) / window, mode="valid") + t_avg = times[window - 1:] + ax.plot(t_avg, running_avg, "r-", linewidth=2.5, alpha=0.8, label=f"Running average (window={window})") + +ax.set_xlabel("Time [ms]", fontsize=13) +ax.set_ylabel("Segregation index $S = (z_d - z_s) / (z_d + z_s)$", fontsize=13) +ax.set_title("Sphere/Dimer Segregation Under Vertical Vibration", fontsize=14) +ax.legend(fontsize=11, loc="lower right") +ax.grid(True, alpha=0.3) +ax.tick_params(labelsize=11) + +fig.tight_layout() +fig.savefig(os.path.join(script_dir, "segregation_index.png"), dpi=150) +print("Saved segregation_index.png") + +# ── Figure 2: COM Trajectories ────────────────────────────────────────── + +fig2, ax2 = plt.subplots(figsize=(8, 5)) +ax2.plot(times, z_sphere, "bo-", markersize=3, linewidth=1, alpha=0.7, label="Spheres (COM $z$)") +ax2.plot(times, z_dimer, "rs-", markersize=3, linewidth=1, alpha=0.7, label="Dimers (COM $z$)") + +ax2.set_xlabel("Time [ms]", fontsize=13) +ax2.set_ylabel("Center-of-mass height [mm]", fontsize=13) +ax2.set_title("COM Height: Spheres vs Dimers", fontsize=14) +ax2.legend(fontsize=11) +ax2.grid(True, alpha=0.3) +ax2.tick_params(labelsize=11) + +fig2.tight_layout() +fig2.savefig(os.path.join(script_dir, "com_trajectories.png"), dpi=150) +print("Saved com_trajectories.png") + +plt.close("all") diff --git a/examples/bench_multisphere_segregation/validate.py b/examples/bench_multisphere_segregation/validate.py new file mode 100644 index 0000000..19cbd7c --- /dev/null +++ b/examples/bench_multisphere_segregation/validate.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +""" +Validate multisphere segregation benchmark. + +Checks on data/segregation.csv: + 1. No NaN/Inf in data + 2. Both sphere and dimer z-positions are positive (physical) + 3. Segregation index S > 0 in the final 20% of the run (dimers above spheres) + 4. Final S is significantly positive (S > 0.005) + 5. Segregation trend: S increases over time (positive slope in second half) +""" + +import os +import sys +import numpy as np + +script_dir = os.path.dirname(os.path.abspath(__file__)) +data_file = os.path.join(script_dir, "data", "segregation.csv") + +if not os.path.isfile(data_file): + print(f"ERROR: {data_file} not found. Run simulation first.") + sys.exit(1) + +data = np.loadtxt(data_file, delimiter=",", skiprows=1) +if data.ndim == 1: + data = data.reshape(1, -1) + +steps = data[:, 0] +times = data[:, 1] +z_sphere = data[:, 2] +z_dimer = data[:, 3] +seg_index = data[:, 4] + +print("=" * 55) +print("Multisphere Segregation Benchmark Validation") +print("=" * 55) + +passed = 0 +total = 0 + +# 1. No NaN/Inf +total += 1 +if np.all(np.isfinite(data)): + print(" No NaN/Inf: PASS") + passed += 1 +else: + print(" No NaN/Inf: FAIL") + +# 2. Physical z-positions (positive) +total += 1 +if np.all(z_sphere > 0) and np.all(z_dimer > 0): + print(f" Positive z: PASS (z_s={z_sphere[-1]:.4e}, z_d={z_dimer[-1]:.4e})") + passed += 1 +else: + print(f" Positive z: FAIL (min z_s={np.min(z_sphere):.4e}, min z_d={np.min(z_dimer):.4e})") + +# 3. S > 0 in final 20% (dimers above spheres at steady state) +total += 1 +n = len(seg_index) +final_portion = seg_index[int(0.8 * n):] +mean_final_S = np.mean(final_portion) +if mean_final_S > 0: + print(f" Final S > 0: PASS (mean final S = {mean_final_S:.6f})") + passed += 1 +else: + print(f" Final S > 0: FAIL (mean final S = {mean_final_S:.6f})") + +# 4. S is significantly positive at end (S > 0.005) +total += 1 +if mean_final_S > 0.005: + print(f" S > 0.005: PASS (S = {mean_final_S:.6f})") + passed += 1 +else: + print(f" S > 0.005: FAIL (S = {mean_final_S:.6f}, expected > 0.005)") + +# 5. Segregation trend: positive slope in second half +total += 1 +if n > 4: + half = n // 2 + s2 = steps[half:] + si2 = seg_index[half:] + coeffs = np.polyfit(s2, si2, 1) + slope = coeffs[0] + # Accept if slope is non-negative or if S is already high + if slope >= 0 or mean_final_S > 0.01: + print(f" Segregation trend: PASS (slope={slope:.3e}, final S={mean_final_S:.6f})") + passed += 1 + else: + print(f" Segregation trend: FAIL (slope={slope:.3e}, expected >= 0)") +else: + print(" Segregation trend: SKIP (not enough data points)") + +print(f"\nResults: {passed}/{total} checks passed") +if passed == total: + print("ALL CHECKS PASSED") +else: + print(f"WARNING: {total - passed} check(s) failed") + +sys.exit(0 if passed == total else 1)