Every surrogate is a function with adjustable numbers inside. Fitting means choosing those numbers so the function matches the runs.

Start with the simplest case: a few design parameters in, one number out, such as the lift coefficient CL from the angle of attack α and the sweep Λ. A model is a formula with parameters (also called weights) that we are free to choose. A plane is one such formula: ŷ = w₀ + w₁α + w₂Λ, with three weights. The hat on ŷ marks a prediction, as opposed to the solver’s value y.

To bend the plane into a curved surface, add more features computed from the inputs: α², αΛ, Λ². The model is still a weighted sum, ŷ = w · φ(x), where φ(x) is the list of features. Because it is linear in the weights, the best weights can be computed directly.

Least squares

"Best" needs a definition. The usual one: the weights that make the average squared difference between prediction and data as small as possible. That average is the loss, here the mean squared error:

The differences yᵢ − ŷᵢ are the residuals. Plot them against each input: if they look like noise, the model has captured the trend; if they rise and fall smoothly, the model is too simple for the data. Neural networks minimise exactly this kind of loss, only with a model that is not linear in its weights, so they must search for the minimum step by step (see loss and gradient descent).

A polynomial fitted this way is a response surface, the classical surrogate of engineering design since Box and Wilson’s work on experimental design in 1951. It is still a sensible first model for a handful of parameters and smooth scalar outputs.

In scikit-learn every model follows the same pattern: create it, call fit(X, y) on the training data, then predict(X_new). A pipeline chains the feature step and the regression into one model.

import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# 36 training designs (alpha, sweep in degrees) and their CL.
# The formula stands in for the solver: use your own runs.
rng = np.random.default_rng(0)
X = np.column_stack([rng.uniform(-2, 8, 36),
                     rng.uniform(20, 36, 36)])
y = (0.14 + 0.08 * X[:, 0])
y *= np.cos(np.radians(X[:, 1] - 20)) ** 0.5

model = make_pipeline(PolynomialFeatures(degree=2),
                      LinearRegression())
# Least squares: choose the weights w.
model.fit(X, y)
# CL at a new design.
print(model.predict([[4.0, 28.0]]))
residuals = y - model.predict(X)
print(f"largest residual {abs(residuals).max():.4f}")
# The six fitted weights.
print(model[-1].coef_.round(4))

A surrogate exists to answer for new designs. Its error on the runs it was fitted to says almost nothing about that.

The ability to predict well on new inputs is called generalisation. To measure it, keep some runs out of the fit and compare the model’s predictions with the solver there. The slider below fits polynomials of growing degree to 10 noisy points and measures the error on 40 points it never saw, as the root mean squared error (RMSE): the square root of the mean squared error, in the units of the output. Two lessons on, Measuring error defines it and its relatives properly.

At low degree the model is too simple to follow the curve: both errors are high. This is underfitting. As the degree grows the training error keeps falling, but past some point the validation error rises again: the polynomial bends to pass through the noise of each training point. This is overfitting. The best model sits at the bottom of the validation curve, not the training curve. The same U-shape appears with neural networks, where size and training time play the role of the degree (see overfitting and regularisation).

Three sets, three jobs

Training

Runs used to fit the weights.

Validation

Runs used to make choices: the degree, the network size, when to stop training. Every choice made on them uses them a little.

Test

Runs kept aside until the very end and used once, to report the error. If you go back and change the model after looking, they have become validation runs.

With few runs, holding out a third of them is costly. k-fold cross-validation splits the runs into k groups (often 5), fits k times, each time validating on a different group, and averages the k errors. Every run is used for validation exactly once.

Finally, decide what you are testing. Held-out runs scattered inside the training range measure interpolation. To learn how the model behaves beyond its data, hold out a whole region instead, such as every run above α 6°: an extrapolation test. Benchmarks such as AirfRANS publish both kinds of split.

import numpy as np
from sklearn.model_selection import (GroupShuffleSplit, KFold,
                                     cross_val_score)
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
X = np.column_stack([rng.uniform(-2, 8, 44),
                     rng.uniform(20, 36, 44)])
y = (0.14 + 0.08 * X[:, 0])
y *= np.cos(np.radians(X[:, 1] - 20)) ** 0.5
y += rng.normal(0, 0.01, 44)

# Small data: 5-fold cross-validation on whole designs,
# to choose the polynomial degree.
cv = KFold(5, shuffle=True, random_state=0)
for degree in [1, 2, 3, 5, 8]:
    model = make_pipeline(PolynomialFeatures(degree),
                          LinearRegression())
    scores = cross_val_score(
        model, X, y, cv=cv,
        scoring="neg_root_mean_squared_error")
    print(degree, f"validation RMSE {-scores.mean():.4f}")

# Point-wise tables: split by run id, so that no design
# lands on both sides. One id per row of the table:
run_id = np.repeat(np.arange(44), 500)
split = GroupShuffleSplit(n_splits=1, test_size=0.2,
                          random_state=0)
rows = np.zeros(len(run_id))
tr, te = next(split.split(rows, groups=run_id))
# set(): no run is shared.
print(set(run_id[tr]) & set(run_id[te]))

One number never tells the whole story. Each common metric answers a different question, and hides something.

Compare predictions ŷ with solver values y over the test runs (and, for fields, over all their points):

MAE

Mean absolute error: the typical size of an error, in the units of the output.

RMSE

Root mean squared error: like MAE but weights large errors more. The square root of the loss from fitting.

Max error

The single worst point. For fields it usually sits where the field changes fastest, such as a suction peak.

R²

1 − (squared error) / (variance of the data): the share of the variation the model explains. 1 is perfect; easy to look good when the data varies a lot.

Relative L2, the metric of field papers

Most field-surrogate papers report the relative L2 error: for each test design, the size of the error field divided by the size of the true field, both measured with the L2 norm (the square root of the sum of squares over all points). The per-design values are then averaged over the test set.

It makes designs with large and small fields comparable. Its weak spot: if the field is close to zero, the denominator is small and the relative error explodes. It also depends on what is predicted: the same model scores differently on Cp and on pressure in pascals, where the large free-stream value inflates ‖y‖. When reading a benchmark table, check the metric, the quantity, the split (interpolation or extrapolation), and whether the ± is a spread over test designs or over training seeds.

Two cousins appear in benchmark tables. VRMSE, the headline metric of The Well, divides the RMSE by the standard deviation of the true field, so 1.0 means no better than predicting the mean. And integrated quantities are reported in their own units: drag often in counts (one count is ΔCD = 0.0001), as a mean absolute error and a maximum.

Weight by area, not by node

Averages over mesh nodes inherit the mesh. A mesh refined near the wall or at the leading edge puts many nodes there, so a plain point average weights those regions heavily, and the far field or the coarse middle of a panel barely counts. To measure the error per square metre of surface (or per cubic metre of volume), weight each point by the area or volume it stands for: the cell areas that the solver already uses to integrate forces. Report which one you used; on the same model they can differ noticeably.

Put errors in context

• Against the spread of the data. An RMSE means little alone. Divide it by the standard deviation of the output to get a fraction of the variation.
• On the integrated quantities. Compute lift and drag from the predicted field and compare them with the solver’s. Small field errors can add up, and drag is a small difference of large pressure forces.
• As a ranking. A design study needs the designs in the right order more than exact values. The Spearman rank correlation between predicted and true drag (1 = same order) checks that; AirfRANS scores it. Small field errors can still scramble the order when the designs differ by little, as the code below shows.
• Against the solver’s own error. A 1 % surrogate error is irrelevant next to a 5 % solver error against the tunnel (see is the training data right?).
• By region. Map the error over the geometry. An average can hide a large error exactly where the design is decided.

import numpy as np
from scipy.stats import spearmanr
from sklearn.metrics import (mean_absolute_error, r2_score,
                             root_mean_squared_error)

# true, pred: (designs, points) surface Cp on 8 test designs;
# area: the surface area each point stands for. In this
# stand-in the error is largest where the mesh is finest.
rng = np.random.default_rng(0)
true = rng.normal(-0.3, 0.37, size=(8, 2000))
area = rng.uniform(0.2, 2.0, size=true.shape)
noise = rng.normal(0, 1, size=true.shape) * 0.02 / area
pred = true + noise + rng.normal(0, 0.01, size=(8, 1))

t, p = true.ravel(), pred.ravel()
rmse = root_mean_squared_error(t, p)
mae = mean_absolute_error(t, p)
max_err = np.abs(pred - true).max()
r2 = r2_score(t, p)
print(f"RMSE {rmse:.3f} MAE {mae:.3f} max {max_err:.3f}"
      f" R² {r2:.3f}")

# Relative L2: one number per design, then the mean.
err = np.linalg.norm(pred - true, axis=1)
rel_l2 = err / np.linalg.norm(true, axis=1)
print(f"rel. L2 {rel_l2.mean():.3f}, worst {rel_l2.max():.3f}")

# Area-weighted RMSE: refined regions no longer dominate.
w = area / area.sum(axis=1, keepdims=True)
rmse_w = np.sqrt((w * (pred - true) ** 2).sum(axis=1))
print(f"area-weighted RMSE {rmse_w.mean():.3f}")

# Ranking: does the model order the designs correctly?
# (Spearman rank correlation of an integrated quantity;
# here the mean Cp stands in for drag.)
drag_true = true.mean(axis=1)
drag_pred = pred.mean(axis=1)
print(f"Spearman {spearmanr(drag_true, drag_pred)[0]:.2f}")

A Gaussian process predicts a value and a standard deviation. Near the runs it is confident; far from them it says so.

A polynomial commits to one curve. A Gaussian process (GP) starts from a whole family of smooth curves, a probability distribution over functions, and keeps those consistent with the runs. At every new design it returns the mean of the remaining curves (the prediction) and their standard deviation σ (how much they disagree). In geostatistics, where it began, the same method is called kriging.

What “smooth” means is set by the kernel, a function that says how strongly the outputs at two designs are related given their distance. The most common one, the squared-exponential or RBF kernel, has a length scale ℓ: designs much closer than ℓ have nearly equal outputs, designs much farther apart are unrelated.

Click to add runs and watch three things. The band collapses at each run and widens between them. Beyond the last run, the mean drifts back to the prior and the band widens to its prior size: the GP knows it is extrapolating. And a short length scale makes the curve wiggly and the uncertainty grow quickly away from data; a long one does the opposite. In practice ℓ, σf and a noise term (for solver or mesh noise, so the curve need not pass exactly through every run) are fitted automatically by maximising how likely the data is under the model.

The button that adds a run where σ is largest is the simplest form of active learning: let the surrogate choose the next simulation where it is least sure.

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# 20 designs (alpha, sweep) and a measured quantity with
# a little solver noise.
rng = np.random.default_rng(0)
X = np.column_stack([rng.uniform(-2, 8, 20),
                     rng.uniform(20, 36, 20)])
y = np.sin(0.4 * X[:, 0]) * np.cos(np.radians(X[:, 1]))
y += rng.normal(0, 0.01, 20)

kernel = 1.0 * RBF(length_scale=[3.0, 10.0])
kernel += WhiteKernel(noise_level=1e-4)
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True,
                              n_restarts_optimizer=3,
                              random_state=0)
# fit() also tunes the length scales and the noise.
gp.fit(X, y)
print(gp.kernel_)

# A design inside the data, then one far outside.
X_new = np.array([[3.0, 28.0], [12.0, 40.0]])
mean, std = gp.predict(X_new, return_std=True)
for m, s in zip(mean, std):
    print(f"{m:.3f} ± {2 * s:.3f}")  # mean ± 2σ

The design inside the data gets a band of about ± 0.02; the one far outside, about ± 0.65.

A GP predicts one number at a time. To predict a whole field, combine it with a compressed description of the field: the next lesson.

A pressure field has thousands of values, but across a family of designs they move together. A few shapes describe most of it.

Take the field of every training run on a shared mesh and lay each one out as a row: the snapshot matrix, of shape (runs × points), as in from runs to arrays. Subtract the mean field. Proper orthogonal decomposition (POD) then finds the shapes, called modes, that best describe how the runs differ from the mean: the first mode captures as much of the variation as any single shape can, the second as much of the rest, and so on. The share of variation each mode captures is often called its energy. The maths is the singular value decomposition, known in statistics as principal component analysis (PCA).

From modes to a surrogate

Once the modes are known, each run is described by K coefficients instead of thousands of values. A POD surrogate (a reduced-order model) regresses those K coefficients on the design inputs, with a Gaussian process or a polynomial, then rebuilds the field from them. Field prediction becomes a handful of scalar predictions. (The figure rebuilds each test design from its true coefficients, to show what the modes can represent; a surrogate has to predict them.)

import numpy as np
from sklearn.decomposition import PCA
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# Y: one row per run, one column per mesh point (same mesh
# in every run). A stand-in field: a suction peak whose
# height and position move with the design.
rng = np.random.default_rng(0)
X = np.column_stack([rng.uniform(-2, 8, 44),
                     rng.uniform(20, 36, 44)])
s = np.linspace(0, 1, 300)

def field(a, w):
    width = 0.08 + 0.004 * a
    peak = np.exp(-((s - 0.02 * w / 28) / width) ** 2)
    return -(0.3 + 0.1 * a) * peak + 0.3 * s * np.cos(
        np.radians(w))

Y = np.array([field(a, w) for a, w in X])
# A little solver noise.
Y += rng.normal(0, 0.003, Y.shape)
X_tr, Y_tr, X_val, Y_val = X[:36], Y[:36], X[36:], Y[36:]

# Modes from the training runs only.
pca = PCA(n_components=3).fit(Y_tr)
# Energy captured by 1, 2 and 3 modes.
print(pca.explained_variance_ratio_.cumsum().round(4))
# (36, 3): three numbers per run.
Z_tr = pca.transform(Y_tr)

kernel = RBF([5.0, 10.0]) + WhiteKernel()
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True)
gp.fit(X_tr, Z_tr)
# Predict 3 coefficients, rebuild 300 values per run.
Y_pred = pca.inverse_transform(gp.predict(X_val))
err = np.linalg.norm(Y_pred - Y_val, axis=1)
rel_l2 = err / np.linalg.norm(Y_val, axis=1)
print(f"relative L2 on validation runs: {rel_l2.mean():.4f}")

Brunton and Kutz ↗ treat the SVD and POD in depth for engineers, with code.