2. Feature Extraction from LS-DYNA d3plot Files for OOD Scoring

Why d3plot Is the Right Source for OOD Features

LS-DYNA writes its primary transient results to the binary d3plot family of files. Each state (time step written to disk) contains the full model state: nodal coordinates and displacements, velocities, accelerations, element stresses and strains, and — critically for OOD scoring — global and part-level energies plus contact interface forces. Because these quantities are produced by the solver itself, they are consistent across runs regardless of how the model was set up, which makes them a stable feature substrate for an out-of-distribution (OOD) detector.

The practical task in this case is to turn a directory of d3plot output into a compact, fixed-length feature vector per simulation, then score that vector against a nominal reference distribution. The scoring math (Mahalanobis distance) is deliberately kept simple here; the engineering effort is in extraction correctness — reading the right states, the right entities, and reducing time series without losing the signal that distinguishes a nominal crash from an off-nominal one.

Feature Set and Physical Rationale

A crash simulation's nominality is best captured by a small set of physically meaningful scalars sampled over the event timeline. The following features are extracted per run:

Each feature is reduced to a small number of scalars so that a run becomes a vector x ∈ ℝ^d with d ≈ 10–12. Keeping d small matters: with a nominal reference set of, say, 30 runs, a high-dimensional covariance estimate becomes ill-conditioned and the Mahalanobis score becomes unstable.

Extraction Pipeline

The extraction is done with lasso-python (the LS-DYNA Python reader), which exposes states, nodal results, and energy blocks without requiring a full LS-PrePost session. The pipeline below is self-contained: it reads a run directory, walks every state, and emits one feature row.

import numpy as np
import lasso
from pathlib import Path

# Node IDs chosen from the model: B-pillar base, front rail tip, door latch
TRACKED_NODES = [10432, 20871, 31550]

def extract_features(run_dir: Path) -> dict:
    d3plot = lasso.dyna.D3plot(str(run_dir / "d3plot"))
    states = d3plot.arrays

    ke, ie, hg, fc = [], [], [], []
    disp_max = {nid: 0.0 for nid in TRACKED_NODES}
    disp_final = {nid: 0.0 for nid in TRACKED_NODES}
    times = []

    for s in states:
        times.append(float(s["time"]))
        # Global energies (scalar per state)
        ke.append(float(s["kinetic_energy"]))
        ie.append(float(s["internal_energy"]))
        hg.append(float(s["hourglass_energy"]))
        # Contact interface forces: sum of resultant magnitudes
        if "interface_force" in s:
            fc.append(float(np.linalg.norm(s["interface_force"], axis=1).max()))
        else:
            fc.append(0.0)
        # Nodal displacement magnitude for tracked nodes
        node_ids = s["node_ids"]
        disp = s["node_displacement"]
        for nid in TRACKED_NODES:
            idx = np.where(node_ids == nid)[0]
            if idx.size:
                mag = float(np.linalg.norm(disp[idx[0]]))
                disp_max[nid] = max(disp_max[nid], mag)
                disp_final[nid] = mag

    ke, ie, hg, fc = map(np.asarray, (ke, ie, hg, fc))
    t_peak = times[int(np.argmax(ke))]

    feats = {
        "ke_peak": ke.max(),
        "ke_t_peak": t_peak,
        "ke_final": ke[-1],
        "ie_final": ie[-1],
        "ie_slope_at_peak": float(np.gradient(ie, times)[int(np.argmax(ke))]),
        "fc_max": fc.max(),
        "fc_t_max": times[int(np.argmax(fc))],
        "hg_final": hg[-1],
        "hg_ratio": hg[-1] / max(ie[-1], 1e-9),
    }
    for nid in TRACKED_NODES:
        feats[f"disp_max_{nid}"] = disp_max[nid]
        feats[f"disp_final_{nid}"] = disp_final[nid]
    return feats

Two extraction details are easy to get wrong and will silently corrupt the OOD score. First, d3plot states are not guaranteed to be uniformly spaced in time; always use the per-state time value rather than an assumed dt when computing slopes. Second, contact interface forces are stored per interface, so the reduction must decide whether to take the maximum across interfaces (as above) or the sum — mixing the two conventions between the reference set and the query run produces a systematic bias that looks like an OOD signal but is not.

Scoring Against the Nominal Reference

With features extracted, the OOD score is the Mahalanobis distance to the nominal reference set. The reference set is built once from a batch of known-nominal runs, and the query run is scored against it:

def fit_reference(feature_rows: list[dict]):
    keys = sorted(feature_rows[0].keys())
    X = np.array([[r[k] for k in keys] for r in feature_rows])
    mu = X.mean(axis=0)
    # Shrinkage stabilizes covariance when n_runs is small
    S = np.cov(X, rowvar=False) + 1e-6 * np.eye(X.shape[1])
    S_inv = np.linalg.inv(S)
    return keys, mu, S_inv

def ood_score(feats: dict, keys, mu, S_inv) -> float:
    x = np.array([feats[k] for k in keys])
    d = x - mu
    return float(np.sqrt(d @ S_inv @ d))

The score is a single non-negative number: small values indicate the run lies near the nominal manifold, large values indicate it does not. The threshold that separates the two is not set here — that calibration is a separate task. What this case must deliver is a trustworthy score, and the figure below shows how the extraction and scoring stages connect, including the two failure points that most often produce a misleading score.

Verifying the Extraction Before Trusting the Score

A score is only as good as the features behind it, so the extraction must be checked against observable solver output before it is used for OOD decisions. Three checks are sufficient in practice:

1. Energy balance. For a closed system, KE + IE + HG should be approximately conserved after the initial impact. If the extracted KE and IE do not sum to a near-constant total, the state indexing or the energy block selection is wrong.
2. Peak timing. The time of peak kinetic energy should fall within the impact window (typically the first few milliseconds for a frontal crash). A ke_t_peak at the final state indicates the extraction is reading states out of order or the run terminated before impact.
3. Reference-set spread. Compute the OOD score of each nominal reference run against the reference set itself. These self-scores should be small and tightly clustered; a wide spread means the nominal set is not homogeneous (mixed load cases, different mesh densities) and the covariance is not representative.

What This Case Delivers

The output of this case is a reproducible extraction function and a scoring function that together convert any LS-DYNA run directory into a single OOD score. The feature vector is deliberately compact and physically interpretable, so that when a run scores high, the individual feature deviations (for example, an elevated hg_ratio or a shifted ke_t_peak) can be inspected to understand why the run was flagged. That interpretability is what makes the score usable downstream: a fallback decision can be justified by pointing at the specific feature that moved, rather than at an opaque distance value.