2. LS-DYNA and PAM-CRASH output formats

Crash simulation solvers produce results in several distinct container formats, and a data pipeline that cannot reliably decode all of them will stall at the first handoff between teams. LS-DYNA and PAM-CRASH — the two dominant explicit solvers in vehicle crash engineering — share a common conceptual model: a binary state database for full-field results, a binary time-history database for scalar channels, and a family of ASCII files for model metadata, energy summaries, and interface forces. This chapter examines each format's internal structure and shows how to parse them in Python without losing precision or misreading time bases.

The Three Result Tiers and Their Roles

Before writing any parser, it helps to classify what each format is actually for. The d3plot family (LS-DYNA) and its PAM-CRASH equivalent, the DSY binary database, store nodal coordinates, displacements, velocities, accelerations, element stresses, and shell resultants at discrete state times. These are the files you need for contour plots, section cuts, and intrusion measurements. The binout file (LS-DYNA) and THS time-history file (PAM-CRASH) store channel-oriented scalar data — accelerometer traces, interface forces, energy components — sampled at every timestep or at a user-defined interval. Finally, ASCII outputs such as glstat, matsum, nodout, rcforc, and secforc provide human-readable tables that are easy to parse but expensive to store at scale.

A practical pipeline usually reads binout/THS first for time-history screening, then selectively opens d3plot/DSY states only for the time windows and parts that matter. This tiered access pattern keeps memory bounded even for models with millions of nodes.

Parsing d3plot with lasso-python

The d3plot format is a Fortran-style binary stream with a fixed header block followed by per-state records. The header encodes the model title, node count, element counts, and the number of states. Each state record contains nodal coordinates (only in the first state), displacements, velocities, and accelerations, followed by element results. Because the record layout depends on flags in the header, hand-rolling a parser is error-prone; the lasso-python library handles the format variations across LS-DYNA versions.

import lasso
from lasso.dyna import D3plot, ArrayType

d3 = D3plot("run01/d3plot")
print("states:", d3.header.n_states)
print("nodes :", d3.header.n_nodes)

# nodal displacement for all states: shape (n_states, n_nodes, 3)
disp = d3.arrays[ArrayType.node_displacement]
print(disp.shape)

# pick the state closest to 80 ms
times = d3.get_timesteps()
idx = int(abs(times - 0.080).argmin())
print("state", idx, "at", times[idx], "s")

# intrusion of node 4521 along global X
node_id = 4521
print("X intrusion:", disp[idx, node_id, 0])

Two details matter in production. First, get_timesteps() returns the state times in seconds, but LS-DYNA writes them in the model's unit system — if the deck uses millisecond-tonne-millimetre units, the returned values are in milliseconds. Always confirm the unit system from the keyword deck before comparing against a target time. Second, node indices in the array are zero-based positions, not the node IDs printed in the deck. The mapping is available through d3.arrays[ArrayType.node_ids]; using a raw node ID as an array index silently returns the wrong node.

Reading binout Channels

The binout file is a self-describing binary tree. The root contains named branches such as glstat, nodout, rcforc, and matsum; each branch holds one or more channels, and each channel is a 2D array whose first column is time. The lasso.dyna.Binout class exposes this tree directly.

from lasso.dyna import Binout
import numpy as np

bo = Binout("run01/binout")
print(bo.read())                      # list all branches

# global energy balance
t_gl, ke, ie = bo.read("glstat", "kinetic_energy",
                       "internal_energy")
print("final KE:", ke[-1], " IE:", ie[-1])

# interface force on crossmember contact
t_rc, fx = bo.read("rcforc", "x_force", "x_force")
peak = np.max(np.abs(fx))
print("peak |Fx| =", peak, "N at t =", t_rc[np.argmax(np.abs(fx))])

Note that bo.read() with no arguments returns the branch names, while bo.read("rcforc") returns the channel names inside that branch. Channel names are case-sensitive and vary between solver versions — x_force in one release may appear as x-force in another. A robust pipeline enumerates channels at runtime and matches with a normalized key rather than hard-coding strings.

ASCII Result Files

ASCII outputs are the easiest to inspect but the most fragile to parse. The glstat file begins with a title line, then a header row of column names, then one row per timestep. Column positions are fixed-width, and some columns contain * or blank fields when a quantity is undefined. A tolerant parser reads the header, strips whitespace, and coerces non-numeric tokens to NaN.

import pandas as pd
import numpy as np

def read_glstat(path):
    with open(path) as f:
        lines = f.readlines()
    # find the header row containing "time"
    hdr = next(i for i, l in enumerate(lines)
               if l.strip().lower().startswith("time"))
    cols = lines[hdr].split()
    df = pd.read_csv(path, skiprows=hdr + 1, sep=r"\s+",
                     names=cols, na_values=["*"],
                     engine="python")
    return df.apply(pd.to_numeric, errors="coerce")

g = read_glstat("run01/glstat")
print(g[["time", "kinetic energy",
         "internal energy"]].tail())

For nodout, the file is organized in blocks: a header line naming the node, then a column header, then rows of time and displacement/velocity/acceleration. Parsing requires detecting block boundaries by the keyword NODE and accumulating each block into a separate DataFrame. The same block structure appears in rcforc and secforc, where each block corresponds to one interface or cross-section.

PAM-CRASH DSY and THS

PAM-CRASH writes its binary state database as a set of DSY files and time histories as THS files. The DSY format uses a directory block that lists dataset offsets, followed by fixed-size records per state. Because no widely adopted pure-Python reader exists, most pipelines either call the vendor's pamcrash2lsdyna converter or parse the directory block manually with struct.unpack.

import struct, numpy as np

def read_ths_header(path):
    with open(path, "rb") as f:
        raw = f.read(64)
    # example layout: int32 n_channels, int32 n_steps,
    # float32 dt, then channel names
    n_ch, n_steps, dt = struct.unpack("<iif", raw[:12])
    return n_ch, n_steps, dt

n_ch, n_steps, dt = read_ths_header("run01/ths")
print(f"{n_ch} channels, {n_steps} steps, dt={dt}")

The exact byte layout depends on the PAM-CRASH version and the endianness of the writing platform. Always validate the header against a known channel count from the deck before trusting the offsets; a mismatch of four bytes will shift every subsequent record and produce plausible but wrong traces.

Choosing the Right Format for Each Pipeline Stage

The common failure mode across all three format families is silent misalignment: a wrong node index, a unit mismatch, or a byte-offset error produces numbers that look reasonable but are wrong. Guard against this by asserting known invariants after every parse — total mass from glstat should match the deck, the first state time should be zero, and the number of states should equal the number of time entries in the binout. These cheap checks catch the majority of parsing defects before they propagate into a report.