2. Core Python stack for CAE automation

A production CAE automation platform sits between two very different worlds. On one side are solvers such as Abaqus, Nastran, LS-DYNA, or OpenFOAM, which consume structured input decks and emit large binary or ASCII result files. On the other side are engineers who need parameter sweeps, DOE matrices, regression checks, and post-processing summaries. The Python scientific stack — NumPy, SciPy, pandas, and h5py — is what makes that bridge practical. Each library owns a distinct layer of the data pipeline, and understanding where one ends and the next begins prevents the most common automation failures: silent unit mismatches, memory blow-ups, and fragile file parsing.

Layer Responsibilities in the Automation Pipeline

NumPy is the numerical substrate. It provides the ndarray, a contiguous, typed, strided memory block that maps directly onto how solvers think about nodal fields: a displacement array of shape (n_nodes, 3), a stress tensor field of shape (n_elements, 6), or a time-history of shape (n_steps, n_channels). Because the buffer is contiguous and typed, it can be memory-mapped, sliced without copying, and passed to compiled kernels with near-zero overhead. In CAE automation, NumPy is where you reshape raw solver output into consistent axes and where you perform vectorized operations such as computing von Mises stress from six-component tensors across a million elements.

SciPy builds on NumPy and supplies the algorithms that turn raw arrays into engineering decisions: scipy.spatial.cKDTree for nearest-node mapping between dissimilar meshes, scipy.interpolate for resampling time histories onto a common time base, scipy.optimize for calibrating material parameters against test data, and scipy.signal for filtering noisy accelerometer channels. SciPy is also the natural home for sparse linear algebra when you need to assemble and solve reduced systems during surrogate model construction.

pandas handles the tabular, labeled side of the workflow. Solver inputs are rarely pure arrays; they are tables of load cases, material cards, thickness zones, and boundary condition sets. pandas gives those tables an index, named columns, and join/groupby semantics, which is exactly what a DOE driver needs when it must cross a parameter matrix with a load-case list and emit one input deck per combination. It is also the right tool for aggregating scalar results — max displacement, first natural frequency, peak stress — into a run summary that can be filtered and compared across hundreds of jobs.

h5py provides persistent, hierarchical storage for the large result datasets that would be wasteful to keep in memory or re-parse from ASCII. HDF5 files are self-describing, support chunking and compression, and allow partial reads via hyperslabs. A single .h5 file can hold the mesh, the time vector, and per-step field data, with attributes recording units, solver version, and coordinate system conventions.

Reading Solver Output Without Loading Everything

The most consequential design decision in a CAE automation platform is how result data is read. A transient run with 500 steps and 2 million nodes produces displacement data on the order of tens of gigabytes. Loading that into a single NumPy array is possible on a large workstation but wasteful when the automation only needs peak values or a handful of probe locations.

h5py solves this with hyperslab selection. If the solver writes results into a structured HDF5 file, you can read a single node across all time steps, or a single step across all nodes, without touching the rest of the dataset:

import h5py
import numpy as np

with h5py.File("run_042/results.h5", "r") as f:
    # Inspect structure before reading
    f.visit(lambda name: print(name))

    # Read one node's displacement history: shape (n_steps, 3)
    probe = f["/results/disp"][:, 18432, :]

    # Read one step's full stress field: shape (n_elem, 6)
    stress_step = f["/results/stress"][120, :, :]

    # Read metadata stored as attributes
    units = f["/results/disp"].attrs["units"]      # e.g. "mm"
    cs    = f["/results/disp"].attrs["coord_sys"]  # e.g. "global_cartesian"

The visit call is not decoration. Solver-exported HDF5 files vary widely in group naming, and a production platform should discover the structure at runtime rather than hard-code paths that break when the solver version changes. Storing units and coordinate-system conventions as attributes is equally important: a displacement array in millimeters and a stress array in MPa are meaningless together unless the platform tracks those conventions explicitly.

Vectorized Field Computation

Once data is in NumPy arrays, engineering quantities should be computed with vectorized expressions rather than Python loops. The following example computes von Mises stress from a six-component tensor field stored as (n_elem, 6) in Voigt order [s11, s22, s33, s12, s23, s13]:

def von_mises(stress):
    s11, s22, s33, s12, s23, s13 = stress.T
    return np.sqrt(0.5 * ((s11 - s22)**2 + (s22 - s33)**2 +
                          (s33 - s11)**2) +
                   3.0 * (s12**2 + s23**2 + s13**2))

vm = von_mises(stress_step)          # shape (n_elem,)
peak = vm.max()
hot_elem = int(vm.argmax())

On a million-element field this executes in milliseconds. The equivalent nested loop in pure Python would take minutes and obscure the intent. The same pattern applies to principal stresses via np.linalg.eigvalsh on the symmetric tensor, to strain energy density, and to any custom failure criterion.

Mapping Between Dissimilar Meshes with SciPy

Automation frequently needs to transfer results from a fine solver mesh to a coarser reporting mesh, or compare a simulation field against experimental measurement points. scipy.spatial.cKDTree handles nearest-neighbor mapping efficiently:

from scipy.spatial import cKDTree

tree = cKDTree(mesh_nodes)                 # (n_nodes, 3) solver mesh
dist, idx = tree.query(probe_points)       # probe_points: (n_probe, 3)
mapped = disp_step[idx]                    # (n_probe, 3)

# Flag probes that fell outside a tolerance
out_of_tol = dist > 2.5   # mm, matching mesh units

The distance check is the verification step. If a probe point maps to a node more than a few element lengths away, the mapped value is not trustworthy and should be reported as a gap rather than silently accepted. This is the kind of check that separates a production platform from a script.

Aggregating Results with pandas

After each run, scalar metrics are collected into a DataFrame. A typical run summary table has one row per job and columns for the DOE parameters plus the extracted responses:

With this table, pandas answers the questions that drive design decisions: df.groupby("thickness_mm")["max_vm_MPa"].max() shows the worst stress per thickness; df[df.max_vm_MPa > 300] isolates overloaded configurations; df.sort_values("f1_Hz") ranks designs by stiffness. The status column is essential — a run that failed to converge must not contribute a misleadingly low stress value to the aggregate.

Failure Symptoms and Diagnostics

Three failure modes appear repeatedly in CAE automation. First, unit mismatch: a solver writes displacement in meters while the platform assumes millimeters, producing a 1000× error that passes every shape check. The fix is to read and assert units from HDF5 attributes at load time. Second, memory exhaustion: reading a full transient field into RAM when only peak values are needed. The fix is hyperslab reads and streaming aggregation. Third, silent parse drift: a solver version changes its output column order, and a hard-coded index now reads the wrong component. The fix is to validate against a known reference case after every solver upgrade and to store the solver version as an HDF5 attribute alongside the data.

These four libraries are not interchangeable alternatives; they are layers. NumPy holds the numbers, SciPy transforms them, pandas organizes the metadata and summaries, and h5py persists the large arrays between runs. A platform that respects those boundaries stays fast, debuggable, and correct as solver versions and dataset sizes grow.