Dynamic SENT: inspect a complete public PhAST example¶

import json
import subprocess
import sys
from pathlib import Path

from IPython.display import Image, display
import matplotlib.pyplot as plt
import numpy as np
import yaml

here = Path.cwd().resolve()
repo_root = next(
    (
        candidate
        for candidate in (here, *here.parents)
        if (candidate / "pyproject.toml").is_file()
        and (candidate / "src" / "phast").is_dir()
    ),
    None,
)
if repo_root is None:
    raise RuntimeError(
        "Start Jupyter from a PhAST repository checkout. "
        "See docs/getting-started.md for installation."
    )

example_dir = repo_root / "examples" / "dynamic" / "B3_dynamic_sent"
required = [
    example_dir / "config.yaml",
    example_dir / "mesh.msh",
    example_dir / "initial_conditions.png",
    example_dir / "damage_final.png",
    example_dir / "damage_evolution.gif",
    example_dir / "energy.csv",
    example_dir / "crack_tip.csv",
    example_dir / "run_metadata.json",
]
missing = [str(path.relative_to(repo_root)) for path in required if not path.is_file()]
if missing:
    raise FileNotFoundError(f"Missing public B3 artifacts: {missing}")

print("Repository:", repo_root)
print("Example:", example_dir.relative_to(repo_root))
print("Required public artifacts:", len(required))

config = yaml.safe_load((example_dir / "config.yaml").read_text(encoding="utf-8"))

print("Problem:", config["problem"]["name"])
print("Geometry units:", config["geometry"]["units"])
print("Named regions:", sorted(config["geometry"]["named_groups"]))
print("Material:")
for key in ("E", "nu", "Gc", "l0", "rho", "eta_residual", "energy_split", "pf_model"):
    print(f"  {key}: {config['material'][key]}")
print("Loading:", config["loading"])
print("Solver:", config["solver"])

import phast

mesh_summary = phast.inspect_mesh(example_dir / "mesh.msh")
print("Points:", mesh_summary["n_points"])
print("Cell blocks:", mesh_summary["cells"])
print("Named groups:")
print(json.dumps(mesh_summary["named_groups"], indent=2))

display(Image(filename=str(example_dir / "initial_conditions.png")))

for index, condition in enumerate(config["boundary_conditions"], start=1):
    print(
        f"{index}. region={condition['nodes']!r}, "
        f"type={condition['type']!r}, "
        f"component={condition.get('component')!r}, "
        f"value={condition.get('value', 0.0)!r}"
    )

expected_regions = {"left", "right", "top", "bottom"}
used_regions = {condition["nodes"] for condition in config["boundary_conditions"]}
if used_regions != expected_regions:
    raise RuntimeError(
        f"Boundary-condition regions {used_regions} do not match {expected_regions}"
    )

solver = config["solver"]
material = config["material"]

assert solver["solver_type"] == "explicit"
assert material["pf_model"] == "AT2"
assert material["energy_split"] == "spectral"
assert material["rho"] > 0.0

print("Mechanics route:", solver["solver_type"])
print("Stable-step safety factor:", solver["dt_safety"])
print("Phase-field model:", material["pf_model"])
print("Energy split:", material["energy_split"])
print("Density:", material["rho"])

validate_command = [
    sys.executable,
    "-m",
    "phast",
    "run",
    str(example_dir / "config.yaml"),
    "--validate-only",
]
print(" ".join(validate_command))
subprocess.run(validate_command, check=True, cwd=repo_root)

RUN_FULL = False
output_dir = repo_root / "runs" / "B3_dynamic_sent"

run_command = [
    sys.executable,
    "-m",
    "phast",
    "run",
    str(example_dir / "config.yaml"),
    "--output_dir",
    str(output_dir),
]

if RUN_FULL:
    subprocess.run(run_command, check=True, cwd=repo_root)
else:
    print("Full solve disabled. To run deliberately, set RUN_FULL = True.")
    print("Command:", " ".join(run_command))

display(Image(filename=str(example_dir / "damage_final.png")))
display(Image(filename=str(example_dir / "damage_evolution.gif")))

energy = np.genfromtxt(
    example_dir / "energy.csv",
    delimiter=",",
    names=True,
    dtype=float,
)
crack_tip = np.genfromtxt(
    example_dir / "crack_tip.csv",
    delimiter=",",
    names=True,
    dtype=float,
)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))
time_us = energy["t_s"] * 1.0e6
for name in ("elastic", "kinetic", "fracture", "total"):
    axes[0].plot(time_us, energy[name], label=name)
axes[0].set_xlabel("Time [microseconds]")
axes[0].set_ylabel("Energy [consistent model units]")
axes[0].set_title("Retained energy history")
axes[0].grid(alpha=0.25)
axes[0].legend()

axes[1].plot(crack_tip["t_us"], crack_tip["crack_tip_x_mm"])
axes[1].set_xlabel("Time [microseconds]")
axes[1].set_ylabel("Crack-tip x position [mm]")
axes[1].set_title("Retained inferred crack-tip position")
axes[1].grid(alpha=0.25)

fig.tight_layout()
plt.show()

metadata = json.loads(
    (example_dir / "run_metadata.json").read_text(encoding="utf-8")
)

print("Retained source revision:", metadata.get("git_hash"))
print("Retained device:", metadata.get("device"))
print("Retained mesh:", metadata.get("mesh"))
print("Retained solver:", metadata.get("solver"))
print("Retained runtime [s]:", metadata.get("total_time_s"))
print("Retained crack time [microseconds]:", metadata.get("crack_time_us"))
print(
    "Residual stiffness: retained =",
    metadata["material"]["eta_residual"],
    ", current YAML =",
    config["material"]["eta_residual"],
)
if metadata["material"]["eta_residual"] != config["material"]["eta_residual"]:
    print(
        "PROVENANCE NOTE: retained evidence and current YAML use different "
        "residual-stiffness values; they must not be described as identical runs."
    )