validate.py

#!/usr/bin/env python3
"""Verification & validation suite for the particlesci solvers.

Each benchmark compares simulation output against an analytical solution or
published experimental data and prints PASS/FAIL with the measured error.

  1. hydrostatic   — incompressibility, volume conservation, settling,
                     hydrostatic (linear) constraint-pressure profile
  2. dambreak      — dam-break front position vs Martin & Moyce (1952)
  3. repose        — angle of repose vs internal friction angle arctan(mu)
  4. infiltration  — flow through a static coarse bed vs the Ergun equation
                     (KNOWN OPEN ITEM: boundary coupling not yet calibrated)

Run:  python3 validate.py [--quick] [--only NAME]
Exit code = number of failed benchmarks.
"""

import argparse
import math
import sys
import time

import numpy as np

from particlesci.fluid import PBFluid
from particlesci.granular import Granular

G = 9.81
RESULTS = []


def check(bench, quantity, measured, expected, tol_rel, unit="", note=""):
    err = abs(measured - expected) / max(abs(expected), 1e-12)
    ok = err <= tol_rel
    RESULTS.append(dict(bench=bench, quantity=quantity, measured=measured,
                        expected=expected, err=err, tol=tol_rel, ok=ok,
                        unit=unit, note=note))
    return ok


def check_abs(bench, quantity, measured, tol_abs, unit="", note=""):
    """For quantities whose ideal value is zero (relative error undefined)."""
    ok = abs(measured) <= tol_abs
    RESULTS.append(dict(bench=bench, quantity=quantity, measured=measured,
                        expected=0.0, err=abs(measured), tol=tol_abs, ok=ok,
                        unit=unit, note=note))
    return ok


# --------------------------------------------------------------------------- #
# 1. Hydrostatic tank
# --------------------------------------------------------------------------- #
def bench_hydrostatic(quick=False):
    """Water block settles in a tank: it must conserve volume, hold its rest
    density (incompressibility), come to rest, and show a linear (hydrostatic)
    constraint-pressure profile with depth."""
    f = PBFluid(0.0, 0.0, 0.6, 0.6, h=0.03 if quick else 0.024)
    f.add_block(0.05, 0.0, 0.55, 0.30)
    poured_volume = f.volume()
    f.settle(max_t=3.0 if quick else 6.0, tol=0.008)

    # incompressibility (interior particles: full kernel support)
    surf = float(np.percentile(f.pos[:, 1], 97))
    inner = f.rho[(f.pos[:, 1] < surf - f.h)
                  & (f.pos[:, 0] > 0.05 + f.h) & (f.pos[:, 0] < 0.55 - f.h)]
    rho_err = abs(inner.mean() - f.rho0) / f.rho0
    check_abs("hydrostatic", "interior density error", rho_err, 0.02,
              note="|rho-rho0|/rho0; PBF target <2%")

    # volume conservation: settled column height == poured area / width
    width = 0.6 - 2 * f.pr
    h_meas = float(np.percentile(f.pos[:, 1], 99)) + f.pr
    h_pred = poured_volume / width
    check("hydrostatic", "column height", h_meas, h_pred, 0.05, "m")

    # settling (a liquid at rest is at rest)
    check_abs("hydrostatic", "settled mean speed", f.mean_speed(), 0.02, "m/s")

    # hydrostatic profile: the constraint multiplier lambda is the solver's
    # pressure response; time-averaged over substeps it must be linear in
    # depth. Average over 0.5 s to beat per-substep noise.
    lam_acc = np.zeros(f.n)
    nacc = 0
    for _ in range(30):
        f.step(1 / 60.0)
        lam_acc += f.lam
        nacc += 1
    lam_avg = lam_acc / nacc
    surf = float(np.percentile(f.pos[:, 1], 97))
    sel = f.pos[:, 1] < surf - f.h
    depth = surf - f.pos[sel, 1]
    lam = lam_avg[sel]
    if len(lam) > 20:
        # bin-average to particle-noise-free profile, then fit
        bins = np.linspace(depth.min(), depth.max(), 10)
        idx = np.digitize(depth, bins)
        d_b, l_b = [], []
        for k in range(1, len(bins) + 1):
            m = idx == k
            if m.sum() >= 3:
                d_b.append(depth[m].mean())
                l_b.append(lam[m].mean())
        d_b, l_b = np.array(d_b), np.array(l_b)
        A = np.vstack([d_b, np.ones_like(d_b)]).T
        coef, res, *_ = np.linalg.lstsq(A, l_b, rcond=None)
        ss_tot = np.sum((l_b - l_b.mean()) ** 2)
        r2 = 1.0 - (res[0] / ss_tot if len(res) and ss_tot > 0 else 1.0)
    else:
        r2 = 0.0
    check("hydrostatic", "pressure-profile linearity R^2", r2, 1.0, 0.10,
          note="time+bin-averaged lambda vs depth")


# --------------------------------------------------------------------------- #
# 2. Dam break (Martin & Moyce 1952)
# --------------------------------------------------------------------------- #
def bench_dambreak(quick=False):
    """Collapse of a water column: dimensionless front position Z = x/a at
    T = t*sqrt(2g/a) compared with Martin & Moyce's experiments (Z(T=2) ~ 3.0,
    the standard SPH validation). The front speed under-predicts and converges
    from below with resolution (Z = 2.26 at h=30mm, 2.39 at h=18mm)."""
    a, h0 = 0.15, 0.30
    f = PBFluid(0.0, 0.0, 1.4, 0.6, h=0.024 if quick else 0.018, xsph=0.0)
    f.add_block(0.0, 0.0, a, h0)
    tscale = math.sqrt(2 * G / a)
    target_T = 2.0
    t, dt = 0.0, 1 / 240.0
    series = []
    while t * tscale < target_T + 0.2:
        f.step(dt)
        t += dt
        series.append((t * tscale, float(f.pos[:, 0].max()) / a))
    Ts = np.array([s[0] for s in series])
    Zs = np.array([s[1] for s in series])
    Z2 = float(np.interp(target_T, Ts, Zs))
    check("dambreak", "front Z at T=2", Z2, 3.0, 0.25,
          note="Martin & Moyce 1952; converges from below with resolution")


# --------------------------------------------------------------------------- #
# 3. Angle of repose
# --------------------------------------------------------------------------- #
def measure_repose(mu, mu_roll, n=320, quick=False):
    g = Granular(0.0, 0.0, 1.6, 1.2, mu=mu, mu_roll=mu_roll, seed=1)
    cx = 0.8
    for i in range(n):        # cheap benchmark: full grain count even in --quick
        r = g.rng.uniform(0.008, 0.012)
        # drop from low height: angular (pentagon) grains poured gently
        g.add_grain(cx + g.rng.uniform(-0.02, 0.02), 0.5, r, rho=1600)
        if i % 4 == 3:
            g.step(1 / 60.0)
    g.settle(max_t=10.0)

    pos, rad = g.positions(), g.radii()
    bins = np.arange(0.0, 1.6, 0.04)
    surface = []
    for b0, b1 in zip(bins[:-1], bins[1:]):
        m = (pos[:, 0] >= b0) & (pos[:, 0] < b1)
        if m.any():
            surface.append(((b0 + b1) / 2, float((pos[m, 1] + rad[m]).max())))
    surface = np.array(surface)
    peak = surface[:, 1].max()
    angles = []
    for side in (surface[surface[:, 0] <= cx], surface[surface[:, 0] >= cx]):
        m = (side[:, 1] > 0.25 * peak) & (side[:, 1] < 0.85 * peak)
        if m.sum() >= 3:
            slope = np.polyfit(side[m, 0], side[m, 1], 1)[0]
            angles.append(math.degrees(math.atan(abs(slope))))
    return float(np.mean(angles)) if angles else 0.0


def bench_repose(quick=False):
    """A poured pile's angle of repose should approximate the internal
    friction angle arctan(mu). Requires ANGULAR grains: circles roll and give
    ~5-11 deg regardless of rolling-resistance surrogates (measured); the
    solver's irregular pentagons interlock and reach ~26 deg."""
    mu = 0.6
    angle = measure_repose(mu, mu_roll=30.0, quick=quick)
    target = math.degrees(math.atan(mu))       # 31.0 deg
    check("repose", "angle of repose", angle, target, 0.25, "deg",
          note=f"mu={mu} -> phi={target:.1f} deg; pentagon grains")


# --------------------------------------------------------------------------- #
# 4. Infiltration through a static coarse bed vs Ergun equation
# --------------------------------------------------------------------------- #
def sample_disc(cx, cy, r, dx):
    """Boundary-particle sampling of a solid disc at spacing dx: a grain must
    be represented by its geometry, not its centre point (Akinci-style)."""
    pts = []
    n_rings = max(1, int(r / dx))
    for k in range(n_rings + 1):
        rr = r * k / n_rings
        n = max(1, int(2 * math.pi * rr / dx))
        for j in range(n):
            a = 2 * math.pi * j / n
            pts.append((cx + rr * math.cos(a), cy + rr * math.sin(a)))
    return pts


def bench_infiltration(quick=False):
    """Gravity-driven flow of water through a STATIC bed of coarse grains,
    compared against the Ergun equation (inertial regime at d=100 mm).

    Geometry matters: grains are disc-sampled with boundary particles, and
    the lattice is STAGGERED — a square lattice has straight vertical pore
    channels (zero tortuosity) and measures ~2.7x too fast. The Darcy flux is
    measured volumetrically (infiltrated volume rate / width)."""
    r_g, spacing = 0.05, 0.125
    porosity = 1.0 - math.pi * r_g ** 2 / spacing ** 2
    d = 2 * r_g
    mu_w = 1.0e-3
    A = 150 * mu_w * (1 - porosity) ** 2 / (porosity ** 3 * d ** 2)
    B = 1.75 * (1 - porosity) * 1000.0 / (porosity ** 3 * d)
    v_pred = (-A + math.sqrt(A * A + 4 * B * 1000.0 * G)) / (2 * B)

    h_fluid = 0.018 if quick else 0.015
    dx_fluid = 0.6 * h_fluid
    pts, top_c = [], 0.10
    for row, cy in enumerate(np.arange(0.10, 0.45, spacing * 0.87)):
        off = (spacing / 2) if row % 2 else 0.0
        top_c = max(top_c, cy)
        for cx in np.arange(0.10 + off, 1.15, spacing):
            pts.extend(sample_disc(cx, cy, r_g, dx_fluid))
    bed = np.array(pts)
    bed_top = top_c + r_g

    f = PBFluid(0.0, 0.0, 1.2, 0.9, h=h_fluid)
    f.add_block(0.3, bed_top + 0.02, 0.9, bed_top + 0.20)
    width = 0.6
    part_vol = f.mass / f.rho0
    t, dt = 0.0, 1 / 120.0
    flux = []
    while t < (0.9 if quick else 1.2):
        f.step(dt, bpos=bed, bmass_mult=1.0)
        t += dt
        flux.append((t, float((f.pos[:, 1] < bed_top).sum()) * part_vol))
    fl = np.array(flux)
    v05, v95 = np.percentile(fl[:, 1], [15, 75])
    adv = fl[(fl[:, 1] > v05) & (fl[:, 1] < v95)]
    q = float(np.polyfit(adv[:, 0], adv[:, 1], 1)[0]) / width if len(adv) > 5 else 0.0
    check("infiltration", "Darcy flux (superficial velocity)", q, v_pred,
          1.0, "m/s",
          note=f"Ergun (eps={porosity:.2f}, d={d*1000:.0f}mm); factor-2 goal")


# --------------------------------------------------------------------------- #
BENCHES = {
    "hydrostatic": bench_hydrostatic,
    "dambreak": bench_dambreak,
    "repose": bench_repose,
    "infiltration": bench_infiltration,
}


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--quick", action="store_true", help="coarser/faster run")
    ap.add_argument("--only", choices=BENCHES, help="run one benchmark")
    args = ap.parse_args()

    names = [args.only] if args.only else list(BENCHES)
    for name in names:
        t0 = time.time()
        print(f"running {name} ...", flush=True)
        BENCHES[name](args.quick)
        print(f"  done in {time.time() - t0:.1f}s")

    print()
    print(f"{'benchmark':<14} {'quantity':<32} {'measured':>10} "
          f"{'expected':>10} {'error':>8}  result")
    print("-" * 88)
    fails = 0
    for r in RESULTS:
        status = "PASS" if r["ok"] else "FAIL"
        fails += 0 if r["ok"] else 1
        print(f"{r['bench']:<14} {r['quantity']:<32} "
              f"{r['measured']:>10.4g} {r['expected']:>10.4g} "
              f"{r['err']*100:>7.1f}%  {status}"
              + (f"   [{r['note']}]" if r["note"] else ""))
    print("-" * 88)
    print(f"{len(RESULTS) - fails}/{len(RESULTS)} checks passed")
    return fails


if __name__ == "__main__":
    sys.exit(main())