From runs to arrays.

A model reads numbers, not result files. The first job is to turn every run into rows of numbers with a fixed layout.

A run has inputs (the design parameters and operating conditions you chose) and outputs (what the solver returned). How you lay them out depends on whether the outputs are a few numbers or a field. NumPy arrays are the common currency: scikit-learn and PyTorch both read them directly.

Scalars: one row per run

For scalar outputs, stack the runs as rows. X holds the inputs, one column per parameter; Y holds the outputs, one column per quantity. Row i of X and row i of Y are the same run.

X: (n_runs, n_params)

The aircraft: (44, 2), the angle of attack and the sweep of each run.

Y: (n_runs, n_outputs)

The aircraft: (44, 2) if you keep the lift and induced-drag coefficients.

Fields: runs × points × channels

When every run uses the same mesh (same points, in the same order), a field fits in a three-dimensional array: one slice per run, one row per point, one column per channel (a physical quantity at that point, such as Cp, or one component of the wall shear stress).

Y: (n_runs, n_points, n_channels)

The aircraft: (44, 7,952, 1), Cp at every surface point of every run. With Cp and three wall-shear components it would be (44, 7,952, 4).

A list of (n_points_i, n_channels)

When the geometry changes from run to run, so does the mesh, and the arrays cannot be stacked. Keep one array per run. How networks deal with that is the subject of grids, meshes and point clouds.

The point-wise table

Some models, such as the neural field on this site, learn from one point at a time. For them, flatten the fields into a long table with one row per point per run: where the point is, which design it belongs to, and the value there. The aircraft gives 44 × 7,952 = 349,888 rows.

The design columns are the site’s first two runs; the point values are for illustration.

import numpy as np
from pathlib import Path

# Stand-in for solver output: params.csv (run, alpha, sweep)
# and one CSV per run with one row per surface point
# (x, y, z, cp). The mesh is the same in every run.
rng = np.random.default_rng(0)
root = Path("runs")
root.mkdir(exist_ok=True)
params = np.column_stack([np.arange(44),
                          rng.uniform(-2, 8, 44),
                          rng.uniform(20, 36, 44)])
np.savetxt(root / "params.csv", params, delimiter=",",
           header="run,alpha,sweep", comments="")
xyz = rng.uniform(-1, 1, (500, 3))
for run, a, s in params:
    cp = 0.1 * a * np.sin(3 * xyz[:, 0]) - 0.01 * s * xyz[:, 1]
    out = np.column_stack([xyz, cp])
    np.savetxt(root / f"run_{int(run):03d}.csv", out,
               delimiter=",")

# Read it back as arrays. X: (runs, params).
P = np.loadtxt(root / "params.csv", delimiter=",", skiprows=1)
run_id, X = P[:, 0].astype(int), P[:, 1:]
F = np.stack([np.loadtxt(root / f"run_{r:03d}.csv",
                         delimiter=",") for r in run_id])
# Y: (runs, points, channels).
Y = F[:, :, 3:]
print(X.shape, F.shape, Y.shape)
# (44, 2) (44, 500, 4) (44, 500, 1)

# The point-wise table: one row per point per run.
n_runs, n_pts, _ = F.shape
table = np.column_stack([
    np.repeat(run_id, n_pts),       # run id
    F[:, :, :3].reshape(-1, 3),     # x, y, z
    np.repeat(X, n_pts, axis=0),    # alpha, sweep
    Y.reshape(-1, 1),               # cp
])
print(table.shape)  # (22000, 7)

With the data laid out, the next question is which runs to make in the first place: sampling the design space.

Simulations are expensive, so each one should teach the surrogate something new.

First choose the design variables (what changes from run to run) and their ranges. The ranges matter: the surrogate will be trusted only inside them. Then decide where in that box to put a limited number of runs. This choice is called the design of experiments.

A regular grid looks natural but scales badly: 10 levels for each of 5 variables is 10⁵ runs, and each level of each variable is repeated thousands of times. Purely random points are cheap but leave clumps and gaps. Space-filling designs spread a fixed number of runs evenly:

Latin hypercube

Cut each variable’s range into N equal slices, then place N samples so that every slice of every variable holds exactly one of them, with the slices paired at random. Seen along any one axis, the samples are perfectly spread. A maximin variant draws many such designs and keeps the one whose two closest samples are farthest apart.

Sobol sequence

A deterministic low-discrepancy sequence, built digit by digit in base 2: every block of 2^m points puts the right number of points in each binary sub-box of the space, so any prefix of the sequence stays evenly spread. You can add points later, which suits budgets that grow; the balance is exact only at powers of two (32, 64, 128 runs). Sobol indices, a sensitivity measure met in the design unit, are a different tool by the same mathematician (and are often estimated with Sobol samples).

How many runs? A common starting point for smooth scalar outputs is about ten runs per design variable. Rough responses, many variables or full fields need more, and the only reliable check is the error on held-out runs (see train, validation and test). The code below compares designs by their discrepancy, a measure of how unevenly the points fill the box (lower is more even).

import numpy as np
from scipy.stats import qmc

# Ranges of alpha and sweep, in degrees.
lo, hi = [-2.0, 20.0], [8.0, 36.0]

# 44 Latin hypercube points in the unit square, then in degrees.
u = qmc.LatinHypercube(d=2, rng=0).random(n=44)
designs = qmc.scale(u, lo, hi)

# 2^5 = 32 Sobol points, and plain random ones to compare.
sobol = qmc.Sobol(d=2, rng=0).random_base2(m=5)
plain = np.random.default_rng(0).random((44, 2))

# Discrepancy: how unevenly the points fill the box.
for name, pts in [("LHS", u), ("Sobol", sobol),
                  ("random", plain)]:
    print(name, round(qmc.discrepancy(pts), 4))

# Hold out validation designs now, before any training.
order = np.random.default_rng(1).permutation(len(designs))
val, train = designs[order[:8]], designs[order[8:]]
print(train.shape, val.shape)  # (36, 2) (8, 2)

The rng= argument needs SciPy 1.15 or later; older versions call it seed=.

Not every run has to be expensive. Many cheap runs (a coarse mesh, a simpler physics model) can be mixed with a few accurate ones; see cheap and expensive data. And instead of fixing all runs up front, the surrogate itself can choose where to simulate next: choosing the next simulation. Before any of the runs train a model, check that they are right: the next lesson, is the training data right?

A surrogate learns whatever the solver says, including its errors. Check the runs before you train on them.

Every lesson after this one compares the surrogate with the solver and calls the solver the truth. It is not: it is a numerical approximation of a physical model, which is itself an approximation of reality. Three checks, familiar from any careful simulation campaign, decide whether the runs deserve to be training data.

1. Is each run converged?

Mesh independence. Run a few representative designs on three meshes, each refined by the same ratio, and watch the quantities you care about. If they still move between the medium and the fine mesh, the dataset will teach the surrogate a mesh error. Richardson extrapolation estimates the mesh-independent value, and the grid convergence index (GCI, Celik et al. 2008) turns the change into an error band for the fine mesh.

Iterative convergence. For each run, check that the residuals have dropped by several orders of magnitude and, more importantly, that the monitored forces have stopped drifting. A run stopped early looks like a valid result and is not. Runs that never settle may be physically unsteady (vortex shedding, buffet), and a steady-state value for them is an arbitrary snapshot.

2. Did any run fail quietly?

In a batch of hundreds of runs, a few always go wrong: a bad mesh, a diverged solution, a setup typo, a different solver version. They rarely crash; they just come out wrong. Look for outliers: plot each integrated quantity against each input, and compare every run with a smooth fit through all the others. A run far from its neighbours is either a real physical change, which you should understand, or an error, which you should fix or drop. Also check that every run used the same setup: the same turbulence model, boundary conditions, reference values and post-processing.

import numpy as np

# 1. Mesh independence: drag coefficient on three meshes,
# each refined by r = 2 in every direction.
cd = {"coarse": 0.02910, "medium": 0.02840, "fine": 0.02820}
f3, f2, f1 = cd["coarse"], cd["medium"], cd["fine"]
r = 2.0
# Observed order of convergence, Richardson extrapolation,
# and the grid convergence index (GCI) of the fine mesh.
p = np.log((f3 - f2) / (f2 - f1)) / np.log(r)
f_exact = f1 + (f1 - f2) / (r ** p - 1)
gci = 1.25 * abs((f2 - f1) / f1) / (r ** p - 1)
print(f"order {p:.2f}, extrapolated CD {f_exact:.5f},"
      f" GCI {gci:.2%}")

# 2. Convergence: final residual and drift of the
# monitored lift over the last 200 iterations, per run.
rng = np.random.default_rng(0)
n_runs = 44
final_residual = 10.0 ** rng.uniform(-7, -5, n_runs)
final_residual[[7, 30]] = [3e-3, 1e-4]
cl_drift = np.abs(rng.normal(0, 2e-4, n_runs))
cl_drift[30] = 0.02
unconverged = (final_residual > 1e-5) | (cl_drift > 1e-3)
print("not converged:", np.flatnonzero(unconverged))

# 3. Outliers: compare each run's lift with a smooth fit
# of all the others (leave one out).
alpha = rng.uniform(-2, 8, n_runs)
cl = 0.14 + 0.08 * alpha + rng.normal(0, 0.005, n_runs)
cl[12] += 0.15
res = np.empty(n_runs)
for i in range(n_runs):
    keep = np.arange(n_runs) != i
    w = np.polyfit(alpha[keep], cl[keep], deg=2)
    res[i] = cl[i] - np.polyval(w, alpha[i])
# A robust spread: 1.4826 x the median absolute deviation.
sigma = 1.4826 * np.median(np.abs(res - np.median(res)))
print("outliers:", np.flatnonzero(np.abs(res) > 5 * sigma))

3. How large is each error?

The error that matters to the design is the surrogate’s error plus the solver’s error against reality. Put numbers on both before tuning the surrogate. Solver-against-test error comes from validation cases (wind tunnel, rig or flight data); mesh error from the GCI; surrogate error from held-out runs (see measuring error).