Add Brazilian disk tensile test benchmark

Cargo.toml

@@ -118,6 +118,10 @@ path = "examples/polymer_melt/main.rs"
 name = "lj_binary"
 path = "examples/lj_binary/main.rs"
 
+[[example]]
+name = "bench_brazilian_disk"
+path = "examples/bench_brazilian_disk/main.rs"
+
 [[example]]
 name = "example_template"
 path = "examples/example_template/main.rs"

examples/bench_brazilian_disk/README.md

@@ -0,0 +1,71 @@
+# Brazilian Disk Tensile Test Benchmark
+
+## Physics
+
+The Brazilian disk test (indirect tensile test) compresses a cylindrical disk
+between two flat platens.  The loading induces a nearly uniform tensile stress
+along the vertical diametral plane, causing failure by a vertical crack through
+the center.
+
+The analytical tensile strength from the peak compressive load *P* is:
+
+$$\sigma_t = \frac{2P}{\pi D t}$$
+
+where *D* is the disk diameter and *t* is the thickness.
+
+## Model
+
+| Parameter | Value |
+|-----------|-------|
+| Disk radius | 10 mm |
+| Particle radius | 0.5 mm |
+| Particle density | 2500 kg/m³ |
+| Young's modulus | 50 MPa |
+| Bond normal stiffness | 10⁸ N/m |
+| Bond break strain | 0.5% |
+| Platen speed | 5 mm/s (each) |
+
+The simulation is **quasi-2D** (one particle layer in the y-direction) with
+~250 particles arranged on a hexagonal close-packed lattice.
+
+### Two-stage approach
+
+1. **Packing** — FIRE energy minimization relaxes the hex lattice so all
+   particles reach equilibrium before bonds are created.
+2. **Loading** — Flat platens move inward at constant velocity, compressing
+   the bonded disk until tensile failure.
+
+## Running
+
+```bash
+cargo run --release --example bench_brazilian_disk --no-default-features \
+    -- examples/bench_brazilian_disk/config.toml
+```
+
+Runtime: ~1–3 minutes in release mode.
+
+## Validation
+
+```bash
+python3 examples/bench_brazilian_disk/validate.py
+```
+
+Checks:
+1. No NaN/Inf in output data
+2. Load builds up (elastic phase exists)
+3. Brittle failure: post-peak load drops below 50% of peak
+4. Bonds break during the test
+5. Tensile strength falls within physically reasonable range (10²–10⁸ Pa)
+
+## Plots
+
+```bash
+python3 examples/bench_brazilian_disk/plot.py
+```
+
+Generates:
+- `load_displacement.png` — Load vs platen displacement with peak annotation
+- `load_bond_breakage.png` — Load and cumulative bond breakage vs time
+
+![Load-Displacement Curve](load_displacement.png)
+![Bond Breakage](load_bond_breakage.png)

examples/bench_brazilian_disk/config.toml

@@ -0,0 +1,94 @@
+# Brazilian Disk Tensile Test Benchmark
+#
+# A bonded hex-lattice disk compressed between two platens until tensile failure.
+# Bond stiffness k_n = 1e5 N/m is chosen to be stable with the auto-computed
+# Rayleigh timestep (~2.9e-6 s). Verlet stability limit: dt < 2/sqrt(k/m) ~ 7e-6 s.
+#
+# Disk: radius = 0.01 m, particle radius = 0.0005 m, hex lattice (~313 particles).
+# Walls start just outside the disk and move inward at 5 mm/s each.
+
+[domain]
+x_low  = -0.020
+x_high =  0.020
+y_low  = -0.002
+y_high =  0.002
+z_low  = -0.020
+z_high =  0.020
+periodic_x = false
+periodic_y = false
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1
+bin_size = 0.002
+
+[output]
+dir = "examples/bench_brazilian_disk"
+
+# ── Material ──────────────────────────────────────────────────────────────
+
+[[dem.materials]]
+name = "rock"
+youngs_mod = 5.0e7       # Young's modulus [Pa]
+poisson_ratio = 0.25
+restitution = 0.5
+friction = 0.5
+
+# ── Walls (platens) ──────────────────────────────────────────────────────
+# Platens start 50 um outside the disk surface (no initial overlap).
+
+# Bottom platen — normal +z, moves up at 5 mm/s.
+[[wall]]
+wall_type = "plane"
+point_x = 0.0
+point_y = 0.0
+point_z = -0.01005
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "rock"
+velocity = [0.0, 0.0, 0.005]
+
+# Top platen — normal -z, moves down at 5 mm/s.
+[[wall]]
+wall_type = "plane"
+point_x = 0.0
+point_y = 0.0
+point_z = 0.01005
+normal_x = 0.0
+normal_y = 0.0
+normal_z = -1.0
+material = "rock"
+velocity = [0.0, 0.0, -0.005]
+
+# ── Bonds ─────────────────────────────────────────────────────────────────
+# Soft bonds compatible with auto-computed Rayleigh timestep.
+# Verlet stability: dt < 2/sqrt(k/m) ~ 7.2e-6 s (auto dt ~ 2.9e-6 s: stable).
+
+[bonds]
+auto_bond = true
+bond_tolerance = 1.05           # 5% tolerance for bonding distance
+normal_stiffness = 1.0e5        # Bond normal stiffness [N/m]
+normal_damping = 5.0            # Bond normal damping [N*s/m]
+tangential_stiffness = 5.0e4    # Bond tangential stiffness [N/m]
+tangential_damping = 2.5        # Bond tangential damping [N*s/m]
+break_normal_stretch = 0.005    # Break at 0.5% normal strain
+break_shear = 0.000005          # Break at 5 um shear displacement [m]
+
+# ── Thermo output ─────────────────────────────────────────────────────────
+
+[thermo]
+columns = ["step", "ke", "bond_strain", "bonds_broken", "walltime"]
+
+# ── Run ───────────────────────────────────────────────────────────────────
+# Auto-computed dt ~ 2.9e-6 s. At 5 mm/s wall speed:
+# - Time to contact: 50 um / 5 mm/s = 0.01 s = ~3500 steps
+# - 500k steps = ~1.44 s total time = ~7.2 mm wall travel each
+# - Compression well exceeds bond break threshold
+
+[run]
+steps = 500000
+thermo = 5000
+
+[dump]
+interval = 50000

examples/bench_brazilian_disk/main.rs

@@ -0,0 +1,187 @@
+//! Brazilian disk (indirect tensile) test benchmark.
+//!
+//! A circular disk of bonded particles is compressed between two flat platens.
+//! The disk fails by a vertical tensile crack and the measured tensile
+//! strength is validated against the analytical Brazilian test formula:
+//! `sigma_t = 2P / (pi D t)`.
+//!
+//! Bond stiffness is chosen to be compatible with the auto-computed Rayleigh
+//! timestep, avoiding the need for manual dt override. Platens start just
+//! outside the disk surface to avoid initial overlaps.
+//!
+//! ```bash
+//! cargo run --release --example bench_brazilian_disk --no-default-features \
+//!     -- examples/bench_brazilian_disk/config.toml
+//! ```
+
+use std::f64::consts::PI;
+use std::fs::{self, File, OpenOptions};
+use std::io::Write as IoWrite;
+
+use mddem::dem_atom::DemAtom;
+use mddem::dem_bond::BondMetrics;
+use mddem::prelude::*;
+
+fn main() {
+    let mut app = App::new();
+    app.add_plugins(CorePlugins)
+        .add_plugins(GranularDefaultPlugins)
+        .add_plugins(DemBondPlugin)
+        .add_plugins(WallPlugin)
+        .add_plugins(FixesPlugin);
+
+    app.add_setup_system(
+        setup_disk.run_if(first_stage_only()),
+        ScheduleSetupSet::Setup,
+    );
+    app.add_update_system(record_load_displacement, ScheduleSet::PostForce);
+
+    app.start();
+}
+
+// ---------------------------------------------------------------------------
+// Setup: hexagonal close-packed disk
+// ---------------------------------------------------------------------------
+
+/// Create a circular disk of particles on a hexagonal lattice, centered at
+/// the origin in the xz-plane.  The y-direction is thin (quasi-2D).
+fn setup_disk(
+    mut atom: ResMut<Atom>,
+    registry: Res<AtomDataRegistry>,
+    material_table: Res<MaterialTable>,
+) {
+    let radius: f64 = 0.0005;  // particle radius [m]
+    let density: f64 = 2500.0; // density [kg/m^3]
+    let disk_radius: f64 = 0.010; // disk outer radius [m]
+
+    let mat_idx = material_table
+        .find_material("rock")
+        .expect("material 'rock' not found");
+
+    let mut dem = registry.expect_mut::<DemAtom>("setup_disk");
+
+    let spacing = 2.0 * radius;
+    let row_height = spacing * (3.0_f64).sqrt() / 2.0;
+
+    let n_rows = (2.0 * disk_radius / row_height).ceil() as i32 + 1;
+    let n_cols = (2.0 * disk_radius / spacing).ceil() as i32 + 1;
+
+    let mut count = 0u32;
+
+    for row in -n_rows..=n_rows {
+        let z = row as f64 * row_height;
+        let x_offset = if row.unsigned_abs() % 2 == 1 {
+            radius
+        } else {
+            0.0
+        };
+
+        for col in -n_cols..=n_cols {
+            let x = col as f64 * spacing + x_offset;
+
+            let dist_from_center = (x * x + z * z).sqrt();
+            if dist_from_center + radius > disk_radius {
+                continue;
+            }
+
+            let mass = density * (4.0 / 3.0) * PI * radius.powi(3);
+            let tag = atom.get_max_tag() + 1;
+
+            atom.natoms += 1;
+            atom.nlocal += 1;
+            atom.tag.push(tag);
+            atom.origin_index.push(0);
+            atom.cutoff_radius.push(radius);
+            atom.is_ghost.push(false);
+            atom.pos.push([x, 0.0, z]);
+            atom.vel.push([0.0; 3]);
+            atom.force.push([0.0; 3]);
+            atom.mass.push(mass);
+            atom.inv_mass.push(1.0 / mass);
+            atom.atom_type.push(mat_idx);
+
+            dem.radius.push(radius);
+            dem.density.push(density);
+            dem.inv_inertia.push(1.0 / (0.4 * mass * radius * radius));
+            dem.quaternion.push([1.0, 0.0, 0.0, 0.0]);
+            dem.omega.push([0.0; 3]);
+            dem.ang_mom.push([0.0; 3]);
+            dem.torque.push([0.0; 3]);
+
+            count += 1;
+        }
+    }
+
+    println!(
+        "BrazilianDisk: placed {} particles in disk of radius {:.4} m",
+        count, disk_radius
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Data recording
+// ---------------------------------------------------------------------------
+
+/// Record wall forces and platen displacement to a CSV file each thermo step.
+fn record_load_displacement(
+    walls: Res<Walls>,
+    atoms: Res<Atom>,
+    run_state: Res<RunState>,
+    run_config: Res<RunConfig>,
+    input: Res<Input>,
+    comm: Res<CommResource>,
+    bond_metrics: Res<BondMetrics>,
+) {
+    if comm.rank() != 0 {
+        return;
+    }
+    let stage = run_config.current_stage(0);
+    let thermo_interval = stage.thermo;
+    if thermo_interval == 0 || run_state.total_cycle % thermo_interval != 0 {
+        return;
+    }
+
+    let base_dir = match input.output_dir.as_deref() {
+        Some(dir) => format!("{}/data", dir),
+        None => "data".to_string(),
+    };
+    fs::create_dir_all(&base_dir).expect("failed to create data directory");
+    let path = format!("{}/load_displacement.csv", base_dir);
+
+    let step = run_state.total_cycle;
+    let time = step as f64 * atoms.dt;
+
+    let f_bottom = walls.planes[0].force_accumulator;
+    let f_top = walls.planes[1].force_accumulator;
+    let load = (f_bottom.abs() + f_top.abs()) / 2.0;
+
+    let z_bottom = walls.planes[0].point_z;
+    let z_top = walls.planes[1].point_z;
+    let gap = z_top - z_bottom;
+
+    let bonds_broken = bond_metrics.total_bonds_broken;
+    let bond_count = bond_metrics.bond_count;
+
+    let mut file = if step == 0 {
+        let mut f = File::create(&path).expect("failed to create load_displacement.csv");
+        writeln!(
+            f,
+            "step,time,load,gap,z_bottom,z_top,f_bottom,f_top,bonds_broken,bond_count"
+        )
+        .expect("write header");
+        f
+    } else {
+        OpenOptions::new()
+            .create(true)
+            .append(true)
+            .open(&path)
+            .expect("failed to open load_displacement.csv")
+    };
+
+    writeln!(
+        &mut file,
+        "{},{:.8e},{:.8e},{:.8e},{:.8e},{:.8e},{:.8e},{:.8e},{},{}",
+        step, time, load, gap, z_bottom, z_top, f_bottom, f_top, bonds_broken, bond_count
+    )
+    .expect("write data line");
+}