2. Reading solver output files

Once a solver has finished a run, the results live in files whose formats were chosen by the solver vendor, not by you. Before any contour plot, probe table, or automated report can be produced, those files must be opened, decoded, and converted into Python objects that preserve the mesh topology, the field values, and the metadata that describes them. This chapter focuses on the four formats that dominate CAE postprocessing workflows: VTK (and its XML variants), CGNS, HDF5, and CSV. Each has a different philosophy, and each demands a different parsing strategy.

Choosing the right reader for each format

The first decision is always whether a mature library already exists. Writing a binary parser by hand is rarely justified when meshio, h5py, pyvista, or pandas can do the job. The table below summarizes the practical mapping between format, library, and typical use.

The choice is not merely about convenience. VTK and CGNS carry mesh connectivity, so a reader must reconstruct cells and their point ordering. HDF5 and CSV are container formats: they store arrays and tables but say nothing about geometry. That distinction determines whether you can plot a field directly or must first map it onto a mesh you already have.

Reading VTK files

Legacy .vtk files begin with a header line such as # vtk DataFile Version 4.2, followed by a title, an encoding keyword (ASCII or BINARY), and a DATASET declaration. The dataset type — UNSTRUCTURED_GRID, STRUCTURED_GRID, POLYDATA — tells the reader how points and cells are laid out. With meshio, the entire structure is exposed in a few lines:

import meshio
mesh = meshio.read("case_042.vtu")
print(mesh.points.shape)          # (n_points, 3)
print(mesh.point_data.keys())     # dict_keys(['pressure', 'velocity'])
print(mesh.cells[0].type)         # 'tetra'
pressure = mesh.point_data["pressure"]

For XML variants such as .vtu, the data may be stored in an appended binary block referenced by offset. meshio and pyvista handle this transparently, but if you inspect the file with a text editor you will see <AppendedData encoding="raw"> followed by unreadable bytes. Do not attempt to parse those bytes manually; the offset table in the <DataArray> elements is the only reliable index.

A common pitfall is confusing point data with cell data. Point data is defined at nodes and is interpolated across cells; cell data is constant per element. If you read mesh.cell_data["wall_shear"] and try to plot it as a nodal field, the result will be misleading. Always check which container the variable lives in before visualization.

Reading CGNS files

CGNS organizes results as a tree of nodes: a CGNSBase_t contains one or more Zone_t nodes, each with a GridCoordinates_t node and one or more FlowSolution_t nodes. The backend may be ADF (the original) or HDF5 (CGNS 3.x and later). meshio can read CGNS directly:

import meshio
cgns = meshio.read("wing_steady.cgns")
print(cgns.points.shape)
for name, arr in cgns.point_data.items():
    print(name, arr.shape)

If meshio cannot resolve a particular node, pyCGNS or the h5py route through the HDF5 backend gives finer control. When the CGNS file uses the HDF5 backend, you can open it with h5py and walk the tree manually, but you must respect the CGNS naming conventions (/Base/Zone/FlowSolution) or downstream tools will reject the data.

Reading HDF5 files

HDF5 is a general-purpose hierarchical container. A solver may write a single .h5 file containing mesh coordinates, connectivity, and dozens of time steps. The h5py library exposes groups like directories and datasets like NumPy arrays:

import h5py
with h5py.File("transient.h5", "r") as f:
    def show(name, obj):
        if isinstance(obj, h5py.Dataset):
            print(name, obj.shape, obj.dtype)
    f.visititems(show)
    coords = f["/mesh/coordinates"][:]
    temp_t0 = f["/results/temperature/step_0000"][:]

Two details matter in practice. First, HDF5 datasets are lazy: slicing f["/results/temperature/step_0000"][:] loads the whole array into memory, while f[...][0:1000] reads only a slab. For large transient cases, always slice. Second, attributes attached to datasets (for example units or time) are accessed with obj.attrs["units"] and should be preserved when you export to another format.

Reading CSV result tables

CSV files rarely contain mesh data; they usually hold scalar monitors such as residual histories, probe values, or force coefficients. pandas.read_csv is the standard entry point, but CAE exports often include comment lines, inconsistent delimiters, or units embedded in the header:

import pandas as pd
df = pd.read_csv("residuals.csv", comment="#", sep=",",
                 skipinitialspace=True)
df.columns = [c.strip().lower() for c in df.columns]
print(df.head())
print(df[["iteration", "continuity"]].tail())

If the file uses semicolons or tabs, set sep accordingly. If the first row contains units rather than names, read with header=None and assign column names explicitly. Always verify that the iteration column is monotonically increasing; a duplicated or reset iteration index usually indicates a restarted run that was concatenated incorrectly.

Format selection and data flow

The diagram below summarizes how the four formats relate to the two questions that drive every postprocessing task: does the file carry mesh topology, and does it carry time history?

A reproducible parsing workflow

Consider a transient CFD case delivered as transient.h5 plus a residuals.csv. The goal is to extract the temperature field at the final step and plot the continuity residual. A minimal, reproducible script looks like this:

import h5py, pandas as pd, numpy as np

# 1. Inspect the HDF5 tree before loading anything
with h5py.File("transient.h5", "r") as f:
    steps = sorted(f["/results/temperature"].keys())
    last = steps[-1]
    temp = f[f"/results/temperature/{last}"][:]
    coords = f["/mesh/coordinates"][:]
    time = f[f"/results/temperature/{last}"].attrs["time"]

# 2. Load the residual monitor
res = pd.read_csv("residuals.csv", comment="#")
res.columns = [c.strip().lower() for c in res.columns]

# 3. Verify before plotting
assert temp.shape[0] == coords.shape[0], "field/point count mismatch"
assert res["iteration"].is_monotonic_increasing, "iteration index not sorted"
print(f"final time = {time:.3f} s, points = {temp.shape[0]}")

The assertions are not decoration. A field/point count mismatch almost always means the HDF5 file stores cell-centered data while the coordinates are nodal, or that a parallel partition was concatenated without reordering. A non-monotonic iteration column usually means two runs were appended. Catching these before plotting saves hours of chasing a contour that looks plausible but is wrong.

Practical checklist

• Identify whether the file carries mesh topology (VTK, CGNS) or only arrays/tables (HDF5, CSV).
• Inspect the HDF5 or CGNS tree before loading full datasets; slice large arrays instead of reading them whole.
• Confirm whether a field is point data or cell data before visualization.
• Preserve units and time attributes when converting between formats.
• Validate array shapes and iteration ordering with assertions before generating any report.

With these readers in place, the next step is to move from raw arrays to structured postprocessing objects — meshes with attached fields, time series with consistent indexing, and tables ready for aggregation. The parsing decisions made here determine how much cleanup is required later.