check_health Namespace Reference

Functions
int	frame_key (Path p)
int	wall_id_from_deck (Path here)
	load_frame (Path out, int fr, int wall_id)
int	main ()
tuple[int, float]	wall_id_and_T (Path here, str|None deck_name)
int	frame_index (Path p)
list[tuple[np.ndarray, int]]	spatial_clusters (np.ndarray X, np.ndarray|None D=None, float link=CLUSTER_LINK, float dam_cut=0.75)
	load_frames (Path out)

Variables
float	R_IN = 0.02
float	R_OUT = 0.021
float	PEN_TOL = 8.0e-5
float	MAX_V = 40.0
float	EARLY_MEAN_DAMAGE_MAX = 0.08
float	MIN_FINAL_DAMAGE = 0.1
int	MIN_DAMAGED_NODES = 10
int	MIN_DAMAGED_PARTICLES = 2
float	MIN_FRAC_DAMAGE = 0.001
	AXIS = np.array([0.0, 0.0])
float	EARLY_DAMAGE_CAP = 0.98
float	EARLY_MEAN_DAMAGE_CAP = 0.25
float	MAX_FINAL_MEAN_DAMAGE = 0.50
int	MIN_DAMAGE_RISE_FRAMES = 1
float	MAX_COM_RADIUS = 0.012
float	MAX_COM_DRIFT = 0.012
float	TETHER_DIST = 1.5e-4
float	TETHER_COM_LIFT = 5.0e-4
float	TETHER_DAMAGE = 0.9
float	CLUSTER_LINK = 1.9e-4
float	MIN_CLUSTER_FRAC = 0.15
float	MIN_CLUSTER_SEP = 6.0e-4
float	MESH_SIZE = 1.0e-4

Detailed Description

Health gates for attrition sim1 (increased Kn_Factor, R_out=0.021).

Containment is measured about the drum axis (0,0) — the rotation center —
not the wall-node COM (protrusion biases COM and falsely flags escape).
IC also checks that the wall mesh was not translated into the packing.

Health gates for attrition sim2 (thin container, R_out=0.0203).

Containment about drum axis (0,0). Wall rotates about an offset, but the
cavity geometry is still centered at the origin.

Health gates for ellipse × tip-up triangle.

Desired: tip-seeded crack that grows into two spatially separated pieces
(not a glued tip tether, not fragment spray).

Function Documentation

frame_index()

int check_health.frame_index

(

Path

p

)

Definition at line 36 of file check_health.py.

def frame_index(p: Path) -> int:
    m = re.search(r"output_(\d+)\.vtu$", p.name)
    return int(m.group(1)) if m else -1
 
 

Referenced by load_frames().

Here is the caller graph for this function:

frame_key()

int check_health.frame_key

(

Path

p

)

Definition at line 30 of file check_health.py.

def frame_key(p: Path) -> int:
    m = re.search(r"output_0_(\d+)_r", p.name)
    return int(m.group(1)) if m else -1
 
 

Referenced by main(), and wall_id_and_T().

Here is the caller graph for this function:

load_frame()

check_health.load_frame	(	Path	out,
		int	fr,
		int	wall_id
	)

Definition at line 46 of file check_health.py.

def load_frame(out: Path, fr: int, wall_id: int):
    xs, ids, dm, vs = [], [], [], []
    for r in sorted(out.glob(f"output_0_{fr}_r*.vtu")):
        m = meshio.read(r)
        xs.append(m.points[:, :2])
        ids.append(m.point_data["Particle_ID"].ravel())
        if "Damage" in m.point_data:
            dm.append(m.point_data["Damage"].ravel())
        if "Velocity" in m.point_data:
            vs.append(m.point_data["Velocity"])
    X = np.concatenate(xs)
    I = np.concatenate(ids).astype(int)
    wall = I == wall_id
    grain = ~wall
    R = np.hypot(X[grain, 0] - AXIS[0], X[grain, 1] - AXIS[1])
    Rw = np.hypot(X[wall, 0] - AXIS[0], X[wall, 1] - AXIS[1])
    D = np.concatenate(dm)[grain] if dm else np.zeros(grain.sum())
    V = np.concatenate(vs)[grain] if vs else np.zeros((grain.sum(), 3))
    Ig = I[grain]
    vmax = float(np.linalg.norm(V, axis=1).max()) if len(V) else 0.0
    n_part_dmg = 0
    for pid in set(Ig.tolist()):
        if D[Ig == pid].max() > 0.1:
            n_part_dmg += 1
    n_com_out = 0
    for pid in set(Ig.tolist()):
        com = X[grain][Ig == pid].mean(axis=0)
        if np.hypot(com[0] - AXIS[0], com[1] - AXIS[1]) > R_IN:
            n_com_out += 1
    return R, D, vmax, Rw, n_part_dmg, n_com_out
 
 

Referenced by main(), and wall_id_and_T().

Here is the caller graph for this function:

load_frames()

check_health.load_frames

(

Path

out

)

Definition at line 83 of file check_health.py.

def load_frames(out: Path):
    files = sorted(out.glob("output_*.vtu"), key=frame_index)
    files = [p for p in files if frame_index(p) >= 0]
    rows = []
    for p in files:
        m = meshio.read(p)
        pid = m.point_data["Particle_ID"].ravel().astype(int)
        X = m.points[:, :2]
        D_all = np.asarray(m.point_data["Damage"]).ravel()
        V_all = np.asarray(m.point_data["Velocity"])
        ell = pid == 1
        tri = pid == 0
        tip = X[tri][np.argmax(X[tri, 1])]
        D = D_all[ell]
        V = V_all[ell]
        Xe = X[ell]
        com = Xe.mean(axis=0)
        R = np.hypot(Xe[:, 0] - com[0], Xe[:, 1] - com[1])
        rmax = float(np.percentile(R, 99)) if len(R) else 0.0
        d_tip = np.linalg.norm(Xe - tip, axis=1)
        j = int(np.argmin(d_tip))
        clusters = spatial_clusters(Xe, D, CLUSTER_LINK)
        n_ref = max(1, int((D < 0.75).sum()) if len(D) else len(Xe))
        n_big = sum(1 for _, sz in clusters if sz >= MIN_CLUSTER_FRAC * n_ref)
        sep = 0.0
        if len(clusters) >= 2 and clusters[1][1] >= MIN_CLUSTER_FRAC * n_ref:
            sep = float(np.linalg.norm(clusters[0][0] - clusters[1][0]))
        rows.append(
            {
                "t_idx": frame_index(p),
                "dmax": float(D.max()) if len(D) else 0.0,
                "dmean": float(D.mean()) if len(D) else 0.0,
                "nd05": int((D > 0.05).sum()),
                "nd50": int((D > 0.5).sum()),
                "vmax": float(np.linalg.norm(V, axis=1).max()) if len(V) else 0.0,
                "rmax": rmax,
                "com": com,
                "tip": tip,
                "min_d_tip": float(d_tip[j]),
                "stuck_dam": float(D[j]),
                "n": int(ell.sum()),
                "n_big_clusters": n_big,
                "cluster_sep": sep,
                "cluster_sizes": [sz for _, sz in clusters[:4]],
            }
        )
    return rows
 
 

References frame_index(), load_frames(), main(), and spatial_clusters().

Referenced by load_frames().

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int check_health.main

(

)

Definition at line 78 of file check_health.py.

def main() -> int:
    here = Path(__file__).resolve().parent
    out = Path(sys.argv[1] if len(sys.argv) > 1 else here / "runs" / "out").resolve()
    wall_id = wall_id_from_deck(here)
 
    ranks0 = sorted(out.glob("output_0_*_r0.vtu"), key=frame_key)
    frames = [frame_key(p) for p in ranks0 if frame_key(p) >= 0]
    if len(frames) < 4:
        print(f"FAIL: need ≥4 frames in {out}, got {len(frames)}")
        return 1
 
    print(
        f"wall_id={wall_id}  R_in={R_IN}  R_out={R_OUT}  pen_tol={PEN_TOL}  "
        f"(radius about drum axis (0,0))"
    )
    print(
        "fr   maxR_ax  n_pen  n_thru  n_com_out  DamMax  DamMean  n_part_dmg  |v|max  wallR"
    )
    series = []
    fails = []
    for fr in frames:
        R, D, vmax, Rw, n_part, n_com = load_frame(out, fr, wall_id)
        n_pen = int((R > R_IN + PEN_TOL).sum())
        n_thru = int((R > R_OUT).sum())
        series.append(
            (
                fr,
                float(R.max()),
                n_pen,
                n_thru,
                n_com,
                float(D.max()),
                float(D.mean()),
                n_part,
                vmax,
                float(Rw.min()),
                float(Rw.max()),
            )
        )
        s = series[-1]
        print(
            f"{fr:4d}  {s[1]:.5f}  {n_pen:5d}  {n_thru:5d}  {n_com:5d}  "
            f"{s[5]:.4f}  {s[6]:.4f}  {n_part:5d}  {vmax:.3f}  "
            f"[{s[9]:.4f},{s[10]:.4f}]"
        )
        if fr == frames[0]:
            if s[9] < 0.013 or s[10] > 0.0225:
                fails.append(
                    f"IC wall mesh shifted: wall R=[{s[9]:.5f},{s[10]:.5f}] "
                    f"(expect ~[0.014,0.021]) — wall site must equal geom.center()"
                )
            if n_pen > 0 or n_thru > 0:
                fails.append(f"IC grain/wall overlap: n_pen={n_pen} n_thru={n_thru}")
        if n_thru > 0:
            fails.append(f"fr={fr}: {n_thru} nodes through wall (R>{R_OUT})")
        elif n_pen > 0:
            fails.append(f"fr={fr}: {n_pen} nodes past R_in+tol (maxR={R.max():.6f})")
        if n_com > 0:
            fails.append(f"fr={fr}: {n_com} particle COMs outside R_in")
 
    early = series[1] if len(series) > 1 else series[0]
    late = series[-1]
    peak_v = max(s[8] for s in series)
 
    if early[8] > MAX_V:
        fails.append(f"early |v|={early[8]:.1f} > {MAX_V}")
    if peak_v > MAX_V:
        fails.append(f"peak |v|={peak_v:.1f} > {MAX_V}")
    if early[6] > EARLY_MEAN_DAMAGE_MAX:
        fails.append(f"early DamMean={early[6]:.4f} > {EARLY_MEAN_DAMAGE_MAX}")
    if late[5] < MIN_FINAL_DAMAGE:
        fails.append(f"final DamMax={late[5]:.4f} < {MIN_FINAL_DAMAGE}")
    n_dmg = int((load_frame(out, late[0], wall_id)[1] > 0.05).sum())
    if n_dmg < MIN_DAMAGED_NODES:
        fails.append(f"damaged nodes={n_dmg} < {MIN_DAMAGED_NODES}")
    if late[7] < MIN_DAMAGED_PARTICLES:
        fails.append(f"damaged particles={late[7]} < {MIN_DAMAGED_PARTICLES}")
    if late[6] < MIN_FRAC_DAMAGE * 0.5 and n_dmg < MIN_DAMAGED_NODES:
        fails.append("insufficient attrition damage")
    if early[5] > 0.95 and early[6] > 0.1:
        fails.append("early frame already mass-ruptured")
 
    if fails:
        seen = set()
        print("FAIL:")
        for f in fails:
            if f not in seen:
                print(" ", f)
                seen.add(f)
        return 1
    print(
        f"PASS: zero escape maxR={late[1]:.5f}, DamMax {early[5]:.3f}→{late[5]:.3f}, "
        f"n_part_dmg={late[7]}, n_dmg_nodes={n_dmg}, peak|v|={peak_v:.2f}"
    )
    return 0
 
 

References frame_key(), load_frame(), main(), and wall_id_from_deck().

Referenced by load_frames(), main(), and wall_id_and_T().

Here is the call graph for this function:

Here is the caller graph for this function:

spatial_clusters()

list[tuple[np.ndarray, int]] check_health.spatial_clusters	(	np.ndarray	X,
		np.ndarray | None	D = None,
		float	link = CLUSTER_LINK,
		float	dam_cut = 0.75
	)

Union-find clusters. Highly damaged nodes are removed so a crack is a cut.

Definition at line 41 of file check_health.py.

) -> list[tuple[np.ndarray, int]]:
    """Union-find clusters. Highly damaged nodes are removed so a crack is a cut."""
    if len(X) == 0:
        return []
    if D is not None:
        alive = D < dam_cut
        if alive.sum() < 4:
            alive = np.ones(len(X), dtype=bool)
        X = X[alive]
    n = len(X)
    parent = np.arange(n)
 
    def find(a: int) -> int:
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a
 
    def union(a: int, b: int) -> None:
        ra, rb = find(a), find(b)
        if ra != rb:
            parent[rb] = ra
 
    for i in range(n):
        for j in range(i + 1, n):
            if np.hypot(X[i, 0] - X[j, 0], X[i, 1] - X[j, 1]) < link:
                union(i, j)
 
    groups: dict[int, list[int]] = {}
    for i in range(n):
        r = find(i)
        groups.setdefault(r, []).append(i)
    out = []
    for idxs in groups.values():
        xi = X[np.asarray(idxs)]
        out.append((xi.mean(axis=0), len(idxs)))
    out.sort(key=lambda t: t[1], reverse=True)
    return out
 
 

Referenced by load_frames().

Here is the caller graph for this function:

wall_id_and_T()

tuple[int, float] check_health.wall_id_and_T	(	Path	here,
		str | None	deck_name
	)

Pick wall id / Final_Time from the deck that was actually run.

Definition at line 35 of file check_health.py.

def wall_id_and_T(here: Path, deck_name: str | None) -> tuple[int, float]:
    """Pick wall id / Final_Time from the deck that was actually run."""
    names = []
    if deck_name:
        names.append(Path(deck_name).name)
    names.extend(["input_short.json", "input_medium.json", "input.json"])
    seen = set()
    for name in names:
        if name in seen:
            continue
        seen.add(name)
        p = here / name
        if not p.is_file():
            continue
        j = json.loads(p.read_text())
        wid = int(j["Displacement_BC"]["Set_1"]["Particle_List"][0])
        T = float(j["Model"]["Final_Time"])
        return wid, T
    return 58, 0.1
 
 

References frame_key(), load_frame(), main(), and wall_id_and_T().

Referenced by wall_id_and_T().

Here is the call graph for this function:

Here is the caller graph for this function:

wall_id_from_deck()

int check_health.wall_id_from_deck

(

Path

here

)

Definition at line 35 of file check_health.py.

def wall_id_from_deck(here: Path) -> int:
    import json
 
    for name in ("input_short.json", "input.json"):
        p = here / name
        if p.is_file():
            j = json.loads(p.read_text())
            return int(j["Displacement_BC"]["Set_1"]["Particle_List"][0])
    return 40
 
 

Referenced by main().

Here is the caller graph for this function:

Variable Documentation

AXIS

check_health.AXIS = np.array([0.0, 0.0])

Definition at line 27 of file check_health.py.

CLUSTER_LINK

float check_health.CLUSTER_LINK = 1.9e-4

Definition at line 30 of file check_health.py.

EARLY_DAMAGE_CAP

float check_health.EARLY_DAMAGE_CAP = 0.98

Definition at line 18 of file check_health.py.

EARLY_MEAN_DAMAGE_CAP

float check_health.EARLY_MEAN_DAMAGE_CAP = 0.25

Definition at line 19 of file check_health.py.

EARLY_MEAN_DAMAGE_MAX

float check_health.EARLY_MEAN_DAMAGE_MAX = 0.08

Definition at line 22 of file check_health.py.

MAX_COM_DRIFT

float check_health.MAX_COM_DRIFT = 0.012

Definition at line 24 of file check_health.py.

MAX_COM_RADIUS

float check_health.MAX_COM_RADIUS = 0.012

Definition at line 23 of file check_health.py.

MAX_FINAL_MEAN_DAMAGE

float check_health.MAX_FINAL_MEAN_DAMAGE = 0.50

Definition at line 21 of file check_health.py.

MAX_V

float check_health.MAX_V = 40.0

Definition at line 21 of file check_health.py.

MESH_SIZE

float check_health.MESH_SIZE = 1.0e-4

Definition at line 33 of file check_health.py.

MIN_CLUSTER_FRAC

float check_health.MIN_CLUSTER_FRAC = 0.15

Definition at line 31 of file check_health.py.

MIN_CLUSTER_SEP

float check_health.MIN_CLUSTER_SEP = 6.0e-4

Definition at line 32 of file check_health.py.

MIN_DAMAGE_RISE_FRAMES

int check_health.MIN_DAMAGE_RISE_FRAMES = 1

Definition at line 22 of file check_health.py.

MIN_DAMAGED_NODES

int check_health.MIN_DAMAGED_NODES = 10

Definition at line 24 of file check_health.py.

MIN_DAMAGED_PARTICLES

int check_health.MIN_DAMAGED_PARTICLES = 2

Definition at line 25 of file check_health.py.

MIN_FINAL_DAMAGE

float check_health.MIN_FINAL_DAMAGE = 0.1

Definition at line 23 of file check_health.py.

MIN_FRAC_DAMAGE

float check_health.MIN_FRAC_DAMAGE = 0.001

Definition at line 26 of file check_health.py.

PEN_TOL

float check_health.PEN_TOL = 8.0e-5

Definition at line 20 of file check_health.py.

R_IN

float check_health.R_IN = 0.02

Definition at line 18 of file check_health.py.

R_OUT

float check_health.R_OUT = 0.021

Definition at line 19 of file check_health.py.

TETHER_COM_LIFT

float check_health.TETHER_COM_LIFT = 5.0e-4

Definition at line 27 of file check_health.py.

TETHER_DAMAGE

float check_health.TETHER_DAMAGE = 0.9

Definition at line 28 of file check_health.py.

TETHER_DIST

float check_health.TETHER_DIST = 1.5e-4

Definition at line 26 of file check_health.py.