The aircraft on this site has 8,000 surface points per run and trains on a laptop. An industrial car run has millions of surface points and up to hundreds of millions of volume cells, and a dataset has hundreds of runs. One run no longer fits in GPU memory, let alone a batch. The models of the previous unit keep the cost per point linear; this lesson is about the training around them.

Sample points, not runs

At each training step, draw a random subset of points from each run, typically 10,000 to 100,000, and compute the loss on those. Over many steps every point is seen, and each step fits in memory. This works because models with query decoders do not need the whole mesh to make a prediction at a point. Two things need care:

• Sampling weights. Uniform sampling of mesh nodes follows the mesh density, high near walls and low in the far field, so the loss is weighted towards refined regions. That often matches what matters; if not, weight each point by its cell volume or area, or sample the regions you care about more.
• Validation on the full mesh. Evaluate held-out runs on all points, in chunks, so the metric matches what you will report (measuring error).

Batching runs of different sizes

A PyTorch DataLoader stacks the items of a batch into one tensor, which fails when runs have different numbers of points. A collate_fn decides how a batch is assembled; here it draws the same number of points from every run.

import torch
from torch import nn
from torch.utils.data import DataLoader

# Stand-in for 20 runs with different point counts.
# Columns: x, y, z, nx, ny, nz, Cp. In practice a Dataset
# would read each run from disk (Zarr, HDF5).
runs = [torch.randn(200_000 + 10_000 * i, 7)
        for i in range(20)]

def sample_points(batch, k=16_384):
    # The same number k of random points from every run.
    picked = torch.stack(
        [run[torch.randint(len(run), (k,))] for run in batch])
    # Inputs (B, k, 6), targets (B, k, 1).
    return picked[..., :6], picked[..., 6:]

loader = DataLoader(runs, batch_size=4, shuffle=True,
                    collate_fn=sample_points)
dev = "cuda" if torch.cuda.is_available() else "cpu"
model = nn.Sequential(nn.Linear(6, 128), nn.SiLU(),
                      nn.Linear(128, 1)).to(dev)
opt = torch.optim.AdamW(model.parameters(), lr=1e-3)

for x, y in loader:                      # one epoch
    x, y = x.to(dev), y.to(dev)
    # Mixed precision: matrix products in bfloat16.
    with torch.autocast(device_type=dev, dtype=torch.bfloat16):
        pred = model(x)
    loss = nn.functional.mse_loss(pred.float(), y)
    opt.zero_grad()
    loss.backward()
    opt.step()
print(x.shape)                           # (4, 16384, 6)

Memory

• Mixed precision. Inside torch.autocast, matrix products run in bfloat16, a 16-bit format with the same range as 32-bit floats. Activations take half the memory and modern GPUs run faster; the weights and optimiser state stay in 32-bit.
• Activation checkpointing. Backpropagation stores every layer's intermediate results. torch.utils.checkpoint keeps only some and recomputes the rest during the backward pass: less memory for about a third more compute.

Several GPUs, and data on disk

Distributed data parallel (DDP) puts a copy of the model on each GPU, gives each a different batch, and averages the gradients after every step: more runs per second, same model. When a single run is too large even for sampling, the points of one mesh can be split across GPUs; Transolver++ does this and only has to share its small slice tokens. On the data side, store runs in chunked formats such as Zarr or HDF5, so a loader can read a random subset of points without loading the whole file.

Budget the compute as you would a CFD campaign: generating the data usually costs more than training, and a training run that does not improve on the validation set after a fixed number of GPU hours is a signal to revisit the data, not to train longer. Companies now sell both halves: GPU-native solvers that produce training data and the surrogates trained on it (Luminary Cloud's Series B).

Simulations come in grades. A coarse mesh, a loose solver tolerance, a RANS model instead of LES, a panel method instead of RANS: each step down costs less and is less accurate. Multi-fidelity training uses many cheap runs to learn the overall shape of the response and a few expensive ones to correct it.

Learn the difference

The simplest scheme is residual or delta learning. Fit a model f_lo to the cheap data. At the few points where you also have an expensive run, compute the difference between the two, and fit a second, small model g to it. The difference is usually smoother and smaller than the field itself, so a few runs pin it down.

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

def lo(a):
    # Cheap model: the right shape, biased values.
    return 0.10 * a + 0.08 * np.sin(1.5 * a)
def hi(a):
    # Expensive solver (a toy stand-in).
    return 1.1 * lo(a) + 0.02 * a - 0.03

def gp():
    return GaussianProcessRegressor(
        RBF(1.0), alpha=1e-6, normalize_y=True)

rng = np.random.default_rng(0)
a_lo = np.linspace(-2, 8, 40)[:, None]   # 40 cheap runs
a_hi = rng.uniform(-2, 8, (5, 1))        # 5 expensive runs

f_lo = gp().fit(a_lo, lo(a_lo[:, 0]))
# Learn only the difference at the expensive runs.
diff = hi(a_hi[:, 0]) - f_lo.predict(a_hi)
g = gp().fit(a_hi, diff)

a = np.linspace(-2, 8, 200)[:, None]
delta = f_lo.predict(a) + g.predict(a)   # y_lo + g(x)
# For comparison: the 5 expensive runs alone.
only_hi = gp().fit(a_hi, hi(a_hi[:, 0])).predict(a)

def rmse(p):
    return np.sqrt(np.mean((p - hi(a[:, 0])) ** 2))
print(f"lo + delta {rmse(delta):.4f}, "
      f"hi only {rmse(only_hi):.4f}")    # about 10x lower

Co-kriging is the Gaussian-process version that fits both levels together and also learns a scale factor ρ between them (y_hi ≈ ρ · y_lo + δ); the SMT toolbox ↗ implements it. For neural surrogates the same idea appears as pretraining on cheap data and fine-tuning on expensive data, or as a fidelity flag in the input.

How much of each

A recent study treated fidelity as a spectrum and asked, for a fixed budget of simulation compute, which mix of cheap and accurate runs gives the lowest error. Across six PDEs and four architectures, the best mix matched the accuracy of high-fidelity-only training for about three times less compute on average, up to 6.3 times; loosening the solver tolerance was the easiest saving to start with. Their plots are error against data-generation cost, one curve per mix: the curve that is lowest at your budget is the one to pick.

Further on this site: Compute-optimal multi-fidelity training on Papers; Quanscient builds a cloud multiphysics platform aimed at generating such training data.

A surrogate trained only on data learns physics approximately: mass is nearly conserved, forces nearly balance. There are two ways to do better. A soft constraint adds a penalty to the loss when a law is broken; the network is encouraged, not forced. A hard constraint builds the law into the model's output, so it holds for any weights.

Soft: a PDE residual in the loss

The network's output is a function of its input coordinates, and autograd (backpropagation) gives its derivatives with respect to them. Plug those derivatives into the governing equation, and its residual (how far the equation is from being satisfied) can be added to the loss at a set of points in the domain. For steady heat conduction, ∇²T = 0:

This is the idea of physics-informed neural networks (PINNs). It helps most when data is scarce and the equation is known exactly. It has costs: second derivatives make each step several times slower, λ needs tuning, and for turbulent flows the averaged equations a surrogate would need are themselves approximate.

Hard: a divergence-free velocity

In incompressible flow, mass conservation says the divergence ∂u/∂x + ∂v/∂y is zero. Instead of predicting u and v, predict a stream function ψ and define u = ∂ψ/∂y and v = −∂ψ/∂x. Then the divergence is ∂²ψ/∂x∂y − ∂²ψ/∂y∂x, which is zero whatever the network does. In 3D the same trick uses a vector potential. AB-UPT offers this for its volume velocity.

import torch
from torch import nn

psi = nn.Sequential(nn.Linear(2, 64), nn.Tanh(),
                    nn.Linear(64, 64), nn.Tanh(),
                    nn.Linear(64, 1))

xy = torch.rand(1000, 2, requires_grad=True)
grad = torch.autograd.grad
dpsi = grad(psi(xy).sum(), xy, create_graph=True)[0]
# u = d psi / dy,  v = - d psi / dx
u, v = dpsi[:, 1], -dpsi[:, 0]

# Divergence du/dx + dv/dy: zero for any weights.
du = grad(u.sum(), xy, create_graph=True)[0]
dv = grad(v.sum(), xy, create_graph=True)[0]
div = du[:, 0] + dv[:, 1]
# Below 1e-6: float32 round-off only.
print(div.abs().max().item())

Symmetries are a related kind of hard constraint. Feeding only relative positions, as in graph networks, makes a model translation-invariant: moving the part changes nothing. Rotation needs more: use only inputs that do not change under rotation, such as distances and angles, or an equivariant network, whose output vectors rotate with its input (for example E(n)-equivariant graph networks). A model with the symmetry built in cannot violate it and needs less data. But the freestream direction and gravity break rotational symmetry in most aerodynamic problems, so build in only the symmetries that truly hold.

Check the physics, not only the error

A small field error can hide a large error in what the designer uses. Small pressure errors of the same sign add up over a wing, so a field with a low RMSE can still be several per cent off in lift. Alongside field metrics, report integrated quantities (forces, moments, heat flux, mass flow through a section) and check them against the solver on held-out runs. AirfRANS was designed with this in mind: its benchmark scores surrogates on lift and drag as well as on the fields.

A prediction without an error bar invites over-trust; an error bar that has never been checked is worse. This lesson covers two practical ways to get one for a neural surrogate, and the test that says whether to believe it.

Two kinds of uncertainty

• Epistemic: the model does not know, because there was little or no data near this input. It shrinks with more runs. This is what matters most for surrogates.
• Aleatoric: the data itself scatters, from turbulence statistics, averaging windows or solver noise. More runs do not remove it.

A Gaussian process gives epistemic uncertainty for free: its σ grows away from the data. Neural networks do not, so it has to be added.

Deep ensembles

Train the same network five times from different random initial weights (and, optionally, different shuffles of the data). Near the data the members agree; away from it they tend to disagree. Use the mean as the prediction and the spread as the uncertainty. It costs five trainings and five evaluations, and it is hard to beat. MC dropout (keep dropout on at prediction time and repeat) is cheaper and usually less reliable.

import torch
from torch import nn

torch.manual_seed(0)
x = torch.rand(36, 2) * 2 - 1     # 36 designs, normalised
y = (x[:, :1] * 3).sin() + 0.5 * x[:, 1:]  # a toy output

def train_one(seed):
    # Different initial weights for each member.
    torch.manual_seed(seed)
    net = nn.Sequential(nn.Linear(2, 64), nn.SiLU(),
                        nn.Linear(64, 64), nn.SiLU(),
                        nn.Linear(64, 1))
    opt = torch.optim.Adam(net.parameters(), lr=1e-2)
    for _ in range(500):
        opt.zero_grad()
        nn.functional.mse_loss(net(x), y).backward()
        opt.step()
    return net

ensemble = [train_one(s) for s in range(5)]
# One design inside the data, one far outside it.
x_new = torch.tensor([[0.2, 0.1], [3.0, 3.0]])
with torch.no_grad():
    preds = torch.stack([net(x_new) for net in ensemble])
print(preds.mean(0).squeeze())    # the predictions
print(preds.std(0).squeeze())     # the spread grows outside

Calibration: do the error bars mean what they say?

An interval is calibrated if a 90 % interval contains the solver's value for about 90 % of held-out runs. Check it: for several levels (50 %, 80 %, 90 %, 95 %), count the share of held-out runs inside. Plotted against the nominal level, the points of a calibrated model lie on the diagonal. Ensembles are often too confident; their spread then needs scaling.

Split conformal prediction

Conformal prediction turns any surrogate into one with calibrated intervals, with a guarantee, using held-out runs. Keep n runs aside as a calibration set. Compute the absolute residual |yᵢ − ŷᵢ| on each. For a target coverage of 1 − α (for example 90 %, α = 0.1), take q̂, a slightly conservative quantile of these residuals. The interval for a new design is ŷ ± q̂.

import numpy as np

rng = np.random.default_rng(1)
# Solver values on 500 calibration runs, and the
# surrogate's predictions for the same runs.
y_cal = rng.normal(size=500)
y_hat_cal = y_cal + rng.normal(0, 0.1, 500)
y_hat_new = 0.37               # prediction for a new design

alpha = 0.1                    # aim for 90 % coverage
r = np.sort(np.abs(y_cal - y_hat_cal))   # residuals
n = len(r)
k = int(np.ceil((n + 1) * (1 - alpha)))  # 451 for n = 500
# The k-th smallest residual; too few runs: no interval.
q = r[k - 1] if k <= n else np.inf
print(f"[{y_hat_new - q:.3f}, {y_hat_new + q:.3f}]")

# Check on fresh runs drawn the same way:
# the share inside should be close to 0.9.
y_test = rng.normal(size=10_000)
y_hat_test = y_test + rng.normal(0, 0.1, 10_000)
print(np.mean(np.abs(y_test - y_hat_test) <= q))

The guarantee: if the new design and the calibration designs are exchangeable (drawn the same way, so no ordering among them is special), the interval contains the true value with probability at least 1 − α, averaged over designs and calibration sets. No assumption about the model or the noise is needed. What it does not promise:

• Nothing outside the calibration distribution. A design beyond the training range is not exchangeable with the calibration runs, and coverage can collapse, as the figure shows.
• Coverage on average, not for every design. The interval has the same width everywhere; it is too wide in easy regions and too narrow in hard ones. Scaling the residuals by an ensemble's spread gives intervals that adapt.
• One quantity at a time. For a whole field, conformalise a summary (the maximum error, or the error in lift), or use methods built for functions.

Coverage is guaranteed to be at least 1 − α, and, with continuous residuals, at most 1 − α + 1/(n + 1), on average over calibration sets. On one calibration set it is itself random: with 40 calibration runs, as in the figure, it typically lands within about ±5 % of the target, which is why the figure shows 96 % for a 90 % target.

Further on this site: Robust predict-then-optimize applies conformal bounds to neural-operator outputs inside a design optimisation; Active operator learning uses cheap predictive uncertainty to choose the next simulations, the subject of the next unit.