A neural network is a model with many weights, built from one simple unit repeated. Fitting it is still regression.

In fitting a response surface we chose the features ourselves (α², αΛ…) and fitted their weights. A neural network also learns the features. Its building block is the neuron: a weighted sum of its inputs plus a constant (the bias), passed through a fixed nonlinear function called the activation.

Many neurons side by side, all reading the same inputs, form a layer; in matrix form, h = φ(Wx + b), with one row of W per neuron. The outputs of one layer are the inputs of the next. A stack of such layers is a multilayer perceptron (MLP). The number of neurons in a layer is its width, the number of layers its depth. The inner layers are called hidden layers; the last layer usually has no activation, so it can output any value.

Why the activation matters

Without activations, each layer is a linear map, and a stack of linear maps is just one linear map: a deep network would be no more capable than the plane of the regression lesson. The bend in the activation is what lets layers combine into curved surfaces. Common choices:

ReLU is the default in many fields. tanh is smooth and bounded. SiLU (also called swish) is smooth like tanh but unbounded like ReLU, a good fit for smooth physical fields; the site’s model uses it.

How capable is an MLP? A theorem (the universal approximation theorem) says that with enough neurons it can approximate any continuous function on a bounded range as closely as you like. It says nothing about how many neurons, how much data, or whether training will find those weights. In practice the limits are data and training, which the rest of this unit is about.

Counting weights

A layer from n inputs to m neurons has n × m weights and m biases. Counting them tells you the size of the model, its memory, and roughly how much data it may need.

In PyTorch a model is a nn.Module: layers are created in __init__ and applied in forward. nn.Linear(n, m) is one layer of weighted sums; nn.Sequential chains layers when the model is a simple stack.

import torch
from torch import nn

torch.manual_seed(0)
# A small MLP: 2 inputs (alpha, sweep), two hidden layers
# of 64 neurons, 1 output (CL).
model = nn.Sequential(
    nn.Linear(2, 64), nn.SiLU(),
    nn.Linear(64, 64), nn.SiLU(),
    nn.Linear(64, 1),
)
# One design, already normalised.
x = torch.tensor([[0.2, -0.5]])
# Untrained, so the number means nothing yet.
print(model(x))
# 4417 weights and biases.
print(sum(p.numel() for p in model.parameters()))

# The field model on this site, written as an nn.Module.
class NeuralField(nn.Module):
    def __init__(self, n_in=32, width=48):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(n_in, width), nn.SiLU(),
            nn.Linear(width, width), nn.SiLU(),
            nn.Linear(width, width), nn.SiLU(),
            nn.Linear(width, 1),
        )

    def forward(self, features):
        return self.net(features)

# 6337
print(sum(p.numel() for p in NeuralField().parameters()))

Training a network means searching for the weights with the lowest loss, one small step at a time.

The loss is one number that measures how far the predictions are from the data. For surrogates it is usually the mean squared error of fitting a response surface. For a polynomial the best weights could be computed in one go; for a network they cannot, because the prediction is not linear in the weights. Instead we start from random weights and improve them repeatedly.

Picture the loss as a landscape over the weights: one dimension per weight, height = loss. The gradient ∇L is the list of slopes of the loss along every weight; it points uphill. Gradient descent takes a small step the other way:

The learning rate η is the most important setting in training. Too small, and progress is slow. Too large, and each step overshoots the valley floor: the weights zigzag from side to side, and beyond a limit they fly off. The figure is the loss of a straight-line fit, y ≈ w x + b, over six points: only two weights, so the landscape can be drawn as contours.

The valley is long and narrow because the two weights affect the loss very differently. Small steps crawl along the floor; large steps bounce between the walls, and above η ≈ 0.44 they diverge. Real networks have thousands to billions of weights and many such valleys at once, which is why the optimisers below adapt the step.

In practice

Minibatches

Computing the gradient over all the data for every step is slow. Stochastic gradient descent (SGD) uses a small random minibatch instead, for example 512 points. Each step is cheaper and slightly noisy, and many more steps fit in the same time. One pass over all the data is an epoch.

Adam

The usual default optimiser. It keeps running averages of each weight’s gradient and squared gradient and scales each weight’s step accordingly, which handles narrow valleys much better than plain SGD. A learning rate around 10⁻³ is a common start.

Schedules

Lower the learning rate during training: large steps early for fast progress, small steps late to settle into the minimum. A cosine decay is a common choice.

import numpy as np
import torch

# Six points on a line, with noise. Find w and b by
# gradient descent on the mean squared error.
x = np.array([0.0, 0.4, 0.8, 1.2, 1.6, 2.0])
y = np.array([0.52, 1.07, 1.73, 2.31, 2.86, 3.55])

w, b, lr = 0.0, 0.0, 0.1
for step in range(500):
    err = (w * x + b) - y
    loss = np.mean(err ** 2)
    # The gradient of the loss...
    dw, db = 2 * np.mean(err * x), 2 * np.mean(err)
    # ...and one step downhill.
    w, b = w - lr * dw, b - lr * db
print(f"w = {w:.3f}, b = {b:.3f}, loss = {loss:.5f}")

# The same in PyTorch: autograd finds the gradient,
# the optimiser takes the step.
xt, yt = torch.tensor(x), torch.tensor(y)
params = torch.zeros(2, dtype=torch.float64,
                     requires_grad=True)
opt = torch.optim.SGD([params], lr=0.1)
for step in range(500):
    pred = params[0] * xt + params[1]
    loss = torch.mean((pred - yt) ** 2)
    opt.zero_grad()
    loss.backward()
    opt.step()
print(params.detach().numpy().round(3))

The PyTorch version never writes the gradient by hand: loss.backward() computes it. How is the next lesson.

Gradient descent needs the slope of the loss for every weight. Backpropagation gets all of them for about the cost of three predictions.

The obvious way to find a slope is to nudge one weight, recompute the loss, and divide. With 6,337 weights that is 6,337 extra predictions per step, and the result is only approximate. Networks in papers have millions of weights. Something better is needed.

A network is a chain of simple operations: multiply, add, activate, repeat, then compare with the data. The chain rule of calculus says that the slope of a chain is the product of the slopes of its links. For a weight w that affects the loss L only through the prediction ŷ:

Backpropagation organises this. The forward pass computes the prediction and stores every intermediate value. The backward pass starts from the loss and walks back through the layers, multiplying local slopes as it goes; at each layer it reads off the gradient for that layer’s weights. One backward pass gives every gradient, at a cost of about two forward passes (three in all, counting the forward pass itself), however many weights there are. Engineers who know adjoint solvers will recognise the idea: backpropagation is reverse-mode differentiation, the same principle as the discrete adjoint.

Autograd in PyTorch

PyTorch records the operations applied to any tensor created with requires_grad=True (model weights have it by default). Calling loss.backward() runs the backward pass and stores each gradient in the tensor’s .grad. Two habits follow. Gradients accumulate, so the optimiser’s zero_grad() clears them before each step. And when you only predict, wrap the code in torch.no_grad() so nothing is recorded.

The same machinery gives gradients with respect to the inputs: how the predicted lift changes with the sweep, for example. These design sensitivities drive optimising on a surrogate, and derivatives of the predicted field with respect to the coordinates are what physics losses penalise.

import torch

# Autograd records every operation, then applies the
# chain rule backwards.
x = torch.tensor([0.0, 0.4, 0.8, 1.2, 1.6, 2.0])
y = torch.tensor([0.52, 1.07, 1.73, 2.31, 2.86, 3.55])
w = torch.tensor(1.0, requires_grad=True)
b = torch.tensor(0.0, requires_grad=True)

loss = torch.mean((w * x + b - y) ** 2)
# Fills w.grad and b.grad.
loss.backward()
err = (w * x + b - y).detach()
# Autograd and the hand-derived gradient agree.
print(w.grad.item(), (2 * err * x).mean().item())
print(b.grad.item(), (2 * err).mean().item())

# Gradients with respect to the INPUT: how the output
# changes with the design.
torch.manual_seed(0)
net = torch.nn.Sequential(torch.nn.Linear(2, 16),
                          torch.nn.Tanh(),
                          torch.nn.Linear(16, 1))
design = torch.tensor([[0.2, -0.5]], requires_grad=True)
out = net(design).sum()
(dout_ddesign,) = torch.autograd.grad(out, design)
# One slope per design variable.
print(dout_ddesign)

Networks train best when every input and output is a number of order one. Scale them before training, and undo it after.

Raw engineering inputs come in very different units and sizes: an angle of attack of a few degrees, a sweep of 20 to 36 degrees, a Reynolds number of millions, a pressure of 10⁵ Pa. A freshly initialised network assumes inputs of order one. Large inputs push its neurons into the flat parts of their activations, and the gradients of different weights end up with wildly different sizes, so no single learning rate suits them all. The cure is cheap: rescale.

Standardise

Subtract the mean and divide by the standard deviation, as above. The usual default.

Min-max

Map the known range onto [−1, 1] or [0, 1]. Natural for design variables with fixed bounds.

Log first

For positive quantities that span decades (Reynolds number, viscosity, a heat flux), take the logarithm, then standardise.

Outputs need the same care. Train on standardised targets and convert predictions back to physical units before reporting any error. Predicting non-dimensional quantities such as Cp (see physical fields) already does much of the work. For fields with several channels (pressure, three velocity components, temperature), use one mean and one standard deviation per channel, computed over all points and runs.

These global statistics, one set for the whole training set, are the default. Some datasets instead scale each run by its own reference values, such as the free-stream speed or the inlet pressure; that is exactly what Cp does, and it is fine because those values are known inputs of the run. What you must not do is scale a run by statistics of its own output field (its mean or maximum Cp): at prediction time that field is the unknown, so the model would be leaning on the answer.

With several output channels, the loss adds up their errors, so the scaling also decides how much each channel counts. Equal weights on standardised channels are a sensible start; if one channel matters more to the design (surface pressure over the volume, say), weight it explicitly in the loss and treat the weights as one more setting to choose on the validation runs.

import numpy as np
import torch
from sklearn.preprocessing import StandardScaler

# alpha, sweep and CL for 44 runs; 36 train, 8 validation.
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]).reshape(-1, 1)
X_tr, X_val, y_tr, y_val = X[:36], X[36:], y[:36], y[36:]

# Statistics from the training runs only.
x_scaler = StandardScaler().fit(X_tr)
y_scaler = StandardScaler().fit(y_tr)
Xn_tr = x_scaler.transform(X_tr)
Xn_val = x_scaler.transform(X_val)
yn_tr = y_scaler.transform(y_tr)
# [0 0] [1 1]
print(Xn_tr.mean(axis=0).round(3), Xn_tr.std(axis=0).round(3))

# After predicting in normalised units, undo the scaling
# to report CL. (yn_pred stands in for the model output.)
yn_pred = yn_tr[:3]
print(y_scaler.inverse_transform(yn_pred).ravel())

# In PyTorch, keep the statistics with the model, so they
# are saved together.
mean = torch.tensor(x_scaler.mean_, dtype=torch.float32)
std = torch.tensor(x_scaler.scale_, dtype=torch.float32)
x_new = (torch.tensor([[4.0, 28.0]]) - mean) / std
print(x_new)

A network can fit its training runs ever more closely while getting worse on new designs. The learning curves show when.

Plot the training loss and the validation loss against the training step: these are the learning curves. Three patterns cover most cases. Both fall together: the network is learning, keep going. Both stay high: it is underfitting, so train longer or use a larger model. Training keeps falling while validation flattens and rises: it is overfitting, fitting the noise and the peculiarities of its training runs, as the high-degree polynomials did in train, validation and test.

Remedies

Early stopping

Keep a copy of the weights whenever the validation loss reaches a new low; stop when it has not improved for a while, and roll back to that copy.

Weight decay

Shrink every weight slightly at each step, which keeps the weights small and favours smoother functions. (AdamW applies the shrinking directly to the weights, separately from the gradient step; it is not a term added to the loss.) Use AdamW with a small weight_decay.

Dropout

During training, set a random fraction of a layer’s outputs to zero, so the network cannot rely on any one neuron. Common in classification, used more sparingly in regression surrogates, where it can blur sharp features.

More data

The surest cure. For a surrogate that means more designs, not more points per design. Where the physics has a symmetry (a mirror plane, a periodic sector), mirrored or rotated copies of runs are free extra data; only use symmetries that truly hold.

A smaller model

Fewer or narrower layers. Parameter counts are only a rough guide, though: large networks often generalise well. Trust the validation curve, not the count.

import copy
import torch
from torch import nn

torch.manual_seed(0)
X_tr = torch.rand(36, 2) * 2 - 1
X_val = torch.rand(8, 2) * 2 - 1
f = lambda X: torch.sin(2 * X[:, :1]) + 0.3 * X[:, 1:] ** 2
# Noisy training labels, clean validation labels.
y_tr = f(X_tr) + 0.05 * torch.randn(36, 1)
y_val = f(X_val)

model = nn.Sequential(
    nn.Linear(2, 256), nn.SiLU(), nn.Dropout(0.1),
    nn.Linear(256, 256), nn.SiLU(), nn.Linear(256, 1))
opt = torch.optim.AdamW(model.parameters(), lr=1e-3,
                        weight_decay=1e-4)
loss_fn = nn.MSELoss()

best, best_state = float("inf"), None
patience, waited = 200, 0
for epoch in range(5000):
    model.train()
    loss = loss_fn(model(X_tr), y_tr)
    opt.zero_grad()
    loss.backward()
    opt.step()
    model.eval()
    with torch.no_grad():
        val = loss_fn(model(X_val), y_val).item()
    if val < best:
        # A new best: keep a copy of the weights.
        best, waited = val, 0
        best_state = copy.deepcopy(model.state_dict())
    else:
        waited += 1
        # No progress for a while: stop.
        if waited > patience:
            break
# Roll back to the best epoch.
model.load_state_dict(best_state)
print(f"stopped at epoch {epoch}, best val MSE {best:.5f}")

The simplest full-field surrogate is an MLP asked about one point at a time. Ask it about every point and you have the field.

So far the networks mapped a design to a few numbers. To predict a field, add the location to the inputs. A neural field (or coordinate network) takes the coordinates of one point, here also the surface normal there, together with the design parameters, and returns the field value at that point: Cp(x, n, α, Λ). Evaluate it at all 8,000 surface points with the same α and Λ and you get the pressure on the whole aircraft; evaluate it anywhere else on the surface and you get a value there too, at any resolution.

It trains on the point-wise table: each minibatch is a random set of points from random training designs, and the loss is the mean squared error on Cp. Everything from the previous lessons applies unchanged.

Fourier features

A plain MLP learns smooth, slowly varying functions first and struggles with sharp features, such as the suction peak at a wing’s leading edge. The remedy is to feed it sines and cosines of the coordinates at several frequencies, called Fourier features. The network then combines ready-made oscillations instead of building them from scratch. The site uses fixed frequencies doubling at each step (octaves, as in NeRF); Tancik et al. also study random frequencies, which the tokens lesson uses.

Strengths and limits

Strengths

Tiny and fast; any resolution, since you can ask about any point; points are independent, so evaluation parallelises perfectly on a GPU.

Limits

Each point is predicted on its own, without looking at its neighbours. The model knows the geometry only through coordinates and a few design parameters, so it covers one parametric family of shapes and generalises poorly to genuinely new ones. Models that read the geometry itself come next: grids, meshes and point clouds, then neural operators, which generalise the idea of querying any point.

import math
import torch
from torch import nn

def fourier(x, octaves=4):
    # sin and cos of each coordinate at π, 2π, 4π, 8π.
    k = 2.0 ** torch.arange(octaves) * math.pi
    # (..., 3, octaves)
    xk = x[..., None] * k
    sc = torch.cat([torch.sin(xk), torch.cos(xk)], dim=-1)
    return sc.flatten(-2)

class NeuralField(nn.Module):
    def __init__(self, width=48):
        super().__init__()
        # 3 + 3 + 2 + 24 = 32 inputs.
        n_in = 3 + 3 + 2 + 3 * 2 * 4
        self.net = nn.Sequential(
            nn.Linear(n_in, width), nn.SiLU(),
            nn.Linear(width, width), nn.SiLU(),
            nn.Linear(width, width), nn.SiLU(),
            nn.Linear(width, 1))

    def forward(self, xyz, normal, design):
        # One row per point.
        f = [xyz, normal, design, fourier(xyz)]
        return self.net(torch.cat(f, dim=-1))

# Stand-in rows of a point-wise table:
# (point, normal, design) -> Cp.
torch.manual_seed(0)
n = 20000
xyz = torch.rand(n, 3) * 2 - 1
normal = nn.functional.normalize(torch.randn(n, 3), dim=-1)
# Normalised alpha and sweep.
design = torch.rand(n, 2) * 2 - 1
peak = torch.exp(-20 * (xyz[:, :1] + 0.8) ** 2)
cp = -(1 + design[:, :1]) * peak + 0.2 * normal[:, 2:]

model = NeuralField()
# 6337
print(sum(p.numel() for p in model.parameters()))
opt = torch.optim.Adam(model.parameters(), lr=3e-3)
for step in range(3000):
    # 512 random points per step.
    i = torch.randint(0, n, (512,))
    pred = model(xyz[i], normal[i], design[i])
    loss = nn.functional.mse_loss(pred, cp[i])
    opt.zero_grad()
    loss.backward()
    opt.step()
print(f"training RMSE {loss.sqrt().item():.3f}")

Fourier features: Tancik et al., NeurIPS 2020 (project page ↗). A close relative that encodes the input function with a second network is DeepONet, met again in neural operators.

Everything so far, in one script, on real solver runs: load, split by design, normalise, train, test, and compare with a baseline.

The data are the site’s own runs, exported from the same solver that drives the 3D aircraft. You need Python 3.10 or later and four packages: pip install torch numpy scipy scikit-learn (PyTorch 2, scikit-learn 1.4 or later). The script runs in about ten seconds on a laptop CPU, with no GPU.

The dataset

aircraft-runs.npz (about 1 MB) holds 56 runs of the aircraft, each with the same 2,000 points of its surface, as NumPy arrays in the layout of from runs to arrays:

alpha, sweep, split, CL: (56,)

The design of each run, its role (0 training, 1 validation, 2 test) and the solver’s lift coefficient. The first 44 runs are the site’s Latin hypercube (36 training, 8 validation); the last 12 are test designs, 4 inside the training box and 8 outside it.

xyz, normal: (56, 2000, 3)

Position of each point and its outward unit normal. The geometry changes with the sweep.

cp, area: (56, 2000)

The pressure coefficient, and the surface area each point stands for, scaled so that Σ area × (…) over the 2,000 points estimates an integral over the whole aircraft. That estimate runs below the stored CL (by about 0.06 on average): this simple solver computes the lift with a vortex lattice and the surface Cp separately, strip by strip, with a 2D panel method, so the two are not one consistent solution as they would be in a single CFD run. Use the integrated lift to compare models with each other, not with the stored CL.

S_ref: (1,)

The reference area for force coefficients.

The recipe

1. Load the arrays, downloading the file once.
2. Split by run with the stored flag, and mark which test designs lie inside the training box (splits).
3. Build the point-wise table and normalise it with statistics of the training runs only (normalising).
4. Train a neural field with Fourier features, Adam and a cosine schedule (neural field, training loop).
5. Build a baseline: four POD modes of the training fields and a Gaussian process on their coefficients (POD).
6. Report per design the RMSE, the relative L2 error and the error in a lift coefficient integrated from the predicted Cp, on the validation designs and on the test designs inside and outside the box (measuring error).

import math
import os
import urllib.request

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

# 1. Data: 56 solver runs x 2,000 surface points (980 KB).
URL = "https://surrogatemodels.com/downloads/aircraft-runs.npz"
FILE = "aircraft-runs.npz"
if not os.path.exists(FILE):
    urllib.request.urlretrieve(URL, FILE)
d = np.load(FILE)
alpha, sweep, split = d["alpha"], d["sweep"], d["split"]
xyz, normal, cp = d["xyz"], d["normal"], d["cp"]
area, S_ref = d["area"], d["S_ref"][0]
n_runs, n_pts = cp.shape

# 2. Split by run: 0 train, 1 validation, 2 test.
train = np.flatnonzero(split == 0)
val = np.flatnonzero(split == 1)
test = np.flatnonzero(split == 2)
in_box = (alpha >= -2) & (alpha <= 8)
in_box &= (sweep >= 20) & (sweep <= 36)
print(len(train), len(val), len(test), "runs")

# 3. Point-wise table and normalisation (training runs only).
D = np.stack([alpha, sweep], axis=1)
design = np.repeat(D[:, None, :], n_pts, axis=1)
X = np.concatenate([xyz, normal, design], axis=-1)
mu = X[train].reshape(-1, 8).mean(axis=0)
sd = X[train].reshape(-1, 8).std(axis=0)
# Normals are unit vectors already: leave them as they are.
mu[3:6], sd[3:6] = 0.0, 1.0
X = torch.tensor((X - mu) / sd, dtype=torch.float32)
Y = torch.tensor(cp, dtype=torch.float32)[..., None]

# 4. The neural field: Fourier features + MLP.
def fourier(p, octaves=4):
    k = 2.0 ** torch.arange(octaves) * math.pi
    pk = p[..., None] * k
    sc = torch.cat([torch.sin(pk), torch.cos(pk)], dim=-1)
    return sc.flatten(-2)

class NeuralField(nn.Module):
    def __init__(self, width=64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(8 + 24, width), nn.SiLU(),
            nn.Linear(width, width), nn.SiLU(),
            nn.Linear(width, width), nn.SiLU(),
            nn.Linear(width, 1))

    def forward(self, x):
        f = torch.cat([x, fourier(x[..., :3])], dim=-1)
        return self.net(f)

torch.manual_seed(0)
model = NeuralField()
steps = 4000
opt = torch.optim.Adam(model.parameters(), lr=3e-3)
sched = torch.optim.lr_scheduler.CosineAnnealingLR(
    opt, T_max=steps, eta_min=1e-4)
Xtr = X[train].reshape(-1, 8)
Ytr = Y[train].reshape(-1, 1)
for step in range(steps):
    i = torch.randint(0, len(Xtr), (1024,))
    loss = nn.functional.mse_loss(model(Xtr[i]), Ytr[i])
    opt.zero_grad()
    loss.backward()
    opt.step()
    sched.step()

with torch.no_grad():
    nf = model(X).squeeze(-1).numpy()

# 5. Baseline: POD of the training fields + a GP on the modes.
pca = PCA(n_components=4).fit(cp[train])
Dn = (D - D[train].mean(axis=0)) / D[train].std(axis=0)
gp = GaussianProcessRegressor(
    kernel=RBF([1.0, 1.0]) + WhiteKernel(1e-4),
    normalize_y=True)
gp.fit(Dn[train], pca.transform(cp[train]))
pod = pca.inverse_transform(gp.predict(Dn))

# 6. Errors per design: RMSE, relative L2 and lift.
def lift(c):
    a = np.radians(alpha)[:, None]
    nL = -normal[..., 0] * np.sin(a)
    nL += normal[..., 1] * np.cos(a)
    return -(c * nL * area).sum(axis=1) / S_ref

def report(name, pred):
    err = pred - cp
    rmse = np.sqrt((err ** 2).mean(axis=1))
    rel = np.linalg.norm(err, axis=1)
    rel /= np.linalg.norm(cp, axis=1)
    dCL = np.abs(lift(pred) - lift(cp))
    print(name)
    groups = [("validation", val),
              ("test, inside box", test[in_box[test]]),
              ("test, outside box", test[~in_box[test]])]
    for label, idx in groups:
        print(f"  {label:18s} RMSE {rmse[idx].mean():.3f}"
              f"  rel. L2 {rel[idx].mean():.1%}"
              f"  |dCL| {dCL[idx].mean():.3f}")

report("neural field", nf)
report("POD + GP baseline", pod)

The script downloads the file from surrogatemodels.com the first time; if you are offline, put aircraft-runs.npz next to the script and it is used directly. Numbers can differ slightly between machines: PyTorch’s random numbers are reproducible on one machine, not across all of them.

What it prints

Three lessons in one table.

• The baseline wins where there is data. The POD + GP is about 36 times more accurate on the validation designs, measured by relative L2. This solver’s pressure distributions are almost exactly a mean plus a few fixed shapes (as the POD lesson showed), which is the ideal case for POD. On RANS data with separation, the gap narrows or reverses; you only know by measuring.
• Both fail outside the box. 18 to 20 % relative L2, with lift errors several times larger than inside. Neither model warns you: the warning has to come from you, by checking each new design against the training range.
• Field and integral errors differ. Outside the box the network’s pressure field is no better than the baseline’s, yet its lift error is about 3 times smaller. Report both, as measuring error recommends.

Next: real RANS data

This dataset is small and its solver idealised. The natural next step is AirfRANS ↗: about a thousand RANS simulations of 2D airfoils at different Reynolds numbers and angles of attack, with full volume fields, a Python loader and published splits for interpolation and for Reynolds and angle-of-attack extrapolation. The same recipe applies: split by simulation, normalise with training statistics, train, and compare with a simple baseline before believing a network. Its meshes differ from run to run, which is where the next unit, grids, meshes and point clouds, begins.