A neural field looks at one point at a time and never sees its neighbours. When the geometry and the field live on a regular grid, a network can do what every finite-difference solver does: combine each cell with the cells around it.

Convolution

A convolution slides a small window, the kernel, over the grid. At every cell it takes a weighted sum of the cell and its neighbours, with the same weights everywhere. A 3×3 kernel has nine weights; a finite-difference Laplacian is one particular choice of them. A convolutional layer learns its weights from data, and it learns many kernels at once: each produces one output channel, a new grid of numbers. The input channels are the grids you start with (geometry, inflow velocity); the hidden channels are features the network invents.

Because the weights are shared across the grid, a convolution has few parameters (a 3×3 layer from 32 to 64 channels has 18,432 weights, whatever the grid size) and it treats the same flow pattern the same way wherever it appears.

Seeing further: pooling and the receptive field

The receptive field of a cell is the patch of the input it can be influenced by. One 3×3 layer sees 3×3 cells; two see 5×5; ten see 21×21. That grows too slowly to connect a leading edge to a wake 100 cells away. Pooling fixes this: a 2×2 max-pool keeps the largest of every four cells and halves the grid in each direction. After three poolings, a 3×3 kernel on the coarse grid spans about 24×24 cells of the original one.

The U-Net

A U-Net is an encoder that convolves and pools down to a coarse grid, followed by a decoder that upsamples back to full resolution. The coarse levels see the whole domain but have lost fine detail, so skip connections copy each encoder level across to the decoder level of the same size, where the two are concatenated. The network combines global context with local detail, much like a multigrid solver. On the way up, a transposed convolution (ConvTranspose2d) does the reverse of pooling: each coarse cell spreads into a 2×2 block through learned weights, doubling the grid in each direction.

import torch
from torch import nn

def block(c_in, c_out):
    # Two 3x3 convolutions: a learned stencil, twice.
    return nn.Sequential(
        nn.Conv2d(c_in, c_out, 3, padding=1), nn.SiLU(),
        nn.Conv2d(c_out, c_out, 3, padding=1), nn.SiLU())

class UNet(nn.Module):
    def __init__(self, c_in=3, c_out=3, w=32):
        super().__init__()
        self.enc1 = block(c_in, w)
        self.enc2 = block(w, 2 * w)
        self.mid = block(2 * w, 4 * w)   # coarsest level
        self.pool = nn.MaxPool2d(2)      # halve H and W
        # Learned upsampling: double H and W.
        self.up2 = nn.ConvTranspose2d(4 * w, 2 * w, 2, stride=2)
        # 4w channels in = 2w (upsampled) + 2w (skip)
        self.dec2 = block(4 * w, 2 * w)
        self.up1 = nn.ConvTranspose2d(2 * w, w, 2, stride=2)
        self.dec1 = block(2 * w, w)
        # A 1x1 convolution: channels to output fields.
        self.out = nn.Conv2d(w, c_out, 1)

    def forward(self, x):                # (B, c_in, H, W)
        e1 = self.enc1(x)                # (B, w, H, W)
        e2 = self.enc2(self.pool(e1))    # (B, 2w, H/2, W/2)
        m = self.mid(self.pool(e2))      # (B, 4w, H/4, W/4)
        # Skip connections: concatenate e2, then e1.
        d2 = self.dec2(torch.cat([self.up2(m), e2], 1))
        d1 = self.dec1(torch.cat([self.up1(d2), e1], 1))
        return self.out(d1)              # (B, c_out, H, W)

# Channels: inflow u, inflow v, airfoil mask.
x = torch.randn(8, 3, 128, 128)
print(UNet()(x).shape)                   # (8, 3, 128, 128)

Train it exactly as in the training loop: the loss is the mean squared error over all cells, ideally masked so that cells inside the body do not count. Grid models overfit like any other network, so watch the validation curve (overfitting).

Limits

• Voxelisation. An unstructured CFD mesh must be resampled onto the grid; anything thinner than a cell is lost.
• Memory grows as N³ in 3D. Halving the cell size multiplies memory by eight, so industrial 3D grids are coarse.
• Surface values are read back by interpolation from cells cut by the wall, exactly where accuracy matters most.
• One resolution. A U-Net is trained for one grid spacing. Neural operators relax this.

A CFD or FE mesh is already a graph: nodes (mesh points) joined by edges (mesh connections). A graph network keeps it as it is, with no grid and no resampling, and lets every node exchange information with the nodes it is connected to.

Features on nodes and edges

Each node carries a vector of node features: a flag saying wall, inlet or fluid, the flow condition, and for time-dependent problems the current field values, but not its absolute position. Each edge carries edge features: the relative position of its two nodes and their distance. Because position enters only through these differences, the model gives the same answer wherever the part sits in space (translation invariance). Rotating the part is another matter, taken up in physics in the loss and in the model.

One round of message passing

In one message-passing layer, every edge computes a message with a small MLP from the features of its two nodes and of the edge. Each node then sums the messages arriving from its neighbours (a sum does not care how many neighbours there are or in which order) and updates its own features with a second MLP.

The update is residual: the layer computes a correction and adds it to what the node already had, rather than replacing it. This makes deep stacks of layers much easier to train. This is a simplified block: MeshGraphNets also updates the edge features in every block (eᵢⱼ ← eᵢⱼ + MLP(eᵢⱼ, hᵢ, hⱼ)) and uses those updated edges as the messages; the reference code linked below shows it.

Encode, process, decode

MeshGraphNets (Pfaff and colleagues, DeepMind) made this recipe standard. An encoder MLP lifts the raw node and edge features to, say, 128 numbers each. A processor of K message-passing blocks, typically 10 to 15, does the work. A decoder MLP turns each node's final features into the output field. For moving structures such as cloth, the paper adds a second set of edges between nodes that are close in space but not connected in the mesh (world-space edges, for contact), and it can refine the mesh as it goes. Every block lets information travel one edge further, so after K blocks a node has heard from its K-hop neighbourhood, and from nothing beyond it.

On a fine mesh, 15 hops may cover a few per cent of the domain: too little to connect a wing tip to its root. Remedies are more blocks, or multiscale graphs that add coarser levels with long edges, the graph counterpart of a U-Net's pooling.

import torch
from torch import nn
from torch_geometric.nn import MessagePassing

def mlp(i, o, h=64):
    return nn.Sequential(nn.Linear(i, h), nn.SiLU(),
                         nn.Linear(h, o))

class Block(MessagePassing):
    def __init__(self, d):
        # Sum the messages arriving at each node.
        super().__init__(aggr="sum")
        self.msg = mlp(3 * d, d)    # from (h_i, h_j, e_ij)
        self.upd = mlp(2 * d, d)    # from (h_i, sum of msgs)

    def forward(self, h, edge_index, e):
        # propagate() calls message(), then sums per node.
        m = self.propagate(edge_index, h=h, e=e)
        # Residual: add a correction to h.
        return h + self.upd(torch.cat([h, m], dim=-1))

    def message(self, h_i, h_j, e):
        # h_i: the receiving node, h_j: the sender
        return self.msg(torch.cat([h_i, h_j, e], dim=-1))

class MeshGNN(nn.Module):
    # Encode, process (K blocks), decode.
    def __init__(self, node_in, edge_in, out, d=64, K=8):
        super().__init__()
        self.enc_n = mlp(node_in, d)
        self.enc_e = mlp(edge_in, d)
        self.blocks = nn.ModuleList(
            Block(d) for _ in range(K))
        self.dec = mlp(d, out)

    def forward(self, x, edge_index, edge_attr):
        h, e = self.enc_n(x), self.enc_e(edge_attr)
        for blk in self.blocks:
            h = blk(h, edge_index, e)
        return self.dec(h)

pos = torch.rand(500, 2)
# In practice: the mesh's own edges.
edge_index = torch.randint(0, 500, (2, 3000))
# Edges carry relative position and distance.
rel = pos[edge_index[1]] - pos[edge_index[0]]
edge_attr = torch.cat(
    [rel, rel.norm(dim=-1, keepdim=True)], dim=-1)
# Nodes carry their type (fluid, wall, inlet), one-hot,
# not their absolute position.
node_type = nn.functional.one_hot(
    torch.randint(0, 3, (500,)), 3).float()
model = MeshGNN(node_in=3, edge_in=3, out=1)
print(model(node_type, edge_index, edge_attr).shape)
# torch.Size([500, 1]): one value per node

On recent PyTorch versions, PyG may print a FutureWarning about torch.jit.script when MessagePassing is imported. It comes from inside the library and is harmless.

For transient problems (vortex shedding, a crash, a cloth), MeshGraphNets predicts the state at the next time step from the current one and is fed its own predictions to go further. How to train such a model so that errors do not run away is the next lesson, time stepping and rollouts.

Watch out: cost grows with the number of edges, the reach is K hops, and on very large 3D meshes the graph does not fit in memory without partitioning. Attention, in the next unit, reaches everywhere in one layer, at a different price.

Many surrogates predict how a field changes in time: vortex shedding behind a cylinder, a crash, the flow in a combustor. A solver marches in small time steps; a surrogate usually does the same with a larger step. It is trained to map the state at time t to the state at t + Δt, and to predict further it is fed its own output: the prediction for t + Δt becomes the input for t + 2Δt. This is an autoregressive rollout, and it is how MeshGraphNets, GNS, the latent rollouts of UPT and most models trained on The Well produce a trajectory.

Why errors grow

A model trained only to predict one step ahead is never shown its own mistakes. At test time it is fed them, and small errors compound: a state slightly off becomes the input for a step the model was never trained on, which drifts further. Errors often grow fastest at high frequencies, where a small one-step error is large relative to the signal. Past a certain number of steps the rollout can leave the range of states seen in training altogether.

The standard remedies

• Predict the change. Learn Δu = u(t + Δt) − u(t) and add it to the current state, rather than predicting the whole next state. The change is small and the network does not have to copy the state through.
• Condition on a short history. Give the model the last two to four states, so it can estimate rates of change, as a multistep time integrator does.
• Noise injection. Add a little random noise to the input states during training (MeshGraphNets, GNS). The model learns to pull a slightly wrong state back towards the right one.
• Unrolled training. Roll the model out for a few steps during training and put the loss on all of them, so it sees its own errors. Backpropagating through every step is expensive; the pushforward trick of Brandstetter and colleagues rolls out without gradients and trains only the last step.
• Refinement and generative models. PDE-Refiner and diffusion surrogates correct the high-frequency part of each step in several passes (pretrained physics models), at the price of several network evaluations per step.

The code below trains a small 1D convolutional model on Burgers' equation (a 1D model of convection and diffusion that forms steep fronts) with a one-step loss and with a four-step unrolled loss, then rolls both out for 100 steps. It runs in a minute or two on a laptop CPU.

import torch
from torch import nn

torch.manual_seed(0)
torch.set_num_threads(1)
n, nu, h = 64, 0.01, 0.002
dx = 1.0 / n

def solver(u, substeps=10):
    # Burgers' equation, explicit, periodic in x.
    # One saved time step = 10 small solver steps.
    for _ in range(substeps):
        up, um = u.roll(-1, -1), u.roll(1, -1)
        ux = (up - um) / (2 * dx)
        uxx = (up - 2 * u + um) / dx**2
        u = u + h * (-u * ux + nu * uxx)
    return u

def trajectories(m, steps):
    x = torch.arange(n) * dx * 2 * torch.pi
    k = torch.randint(1, 4, (m, 1)).float()
    amp = 0.2 + 0.5 * torch.rand(m, 1)
    u = amp * torch.sin(k * x + 6.28 * torch.rand(m, 1))
    out = [u]
    for _ in range(steps):
        u = solver(u)
        out.append(u)
    # (m, steps + 1, n): one row per saved time step
    return torch.stack(out, 1)

def step(model, u):
    # Predict the change du, not the next state.
    return u + model(u[:, None])[:, 0]

def train(unroll, iters=1500):
    c = dict(padding=2, padding_mode="circular")
    model = nn.Sequential(
        nn.Conv1d(1, 32, 5, **c), nn.SiLU(),
        nn.Conv1d(32, 32, 5, **c), nn.SiLU(),
        nn.Conv1d(32, 1, 5, **c))
    opt = torch.optim.Adam(model.parameters(), lr=1e-3)
    T = data.shape[1]
    for _ in range(iters):
        # 64 random (run, start time) pairs
        r = torch.randint(0, len(data), (64,))
        t = torch.randint(0, T - unroll, (64,))
        u, loss = data[r, t], 0.0
        for j in range(1, unroll + 1):
            # unroll > 1: the model is fed its own output
            u = step(model, u)
            loss = loss + ((u - data[r, t + j]) ** 2).mean()
        opt.zero_grad()
        (loss / unroll).backward()
        opt.step()
    return model

def rollout_error(model, runs, steps=100):
    u, errs = runs[:, 0], []
    with torch.no_grad():
        for t in range(1, steps + 1):
            u = step(model, u)
            e = (u - runs[:, t]).norm(dim=1)
            rel = e / runs[:, t].norm(dim=1)
            errs.append(rel.mean().item())
    return errs

data = trajectories(64, 60)
test = trajectories(16, 100)
for unroll in (1, 4):
    errs = rollout_error(train(unroll), test)
    # relative L2 error after 1, 10, 30, 60, 100 steps
    print(unroll, [round(errs[t - 1], 3)
                   for t in (1, 10, 30, 60, 100)])

On this run the one-step model's error goes from 0.4 % after one step to 21 % after 60 and 700 % after 100; the four-step model ends at 17 % and 480 %. Change the seed and the ranking can flip: on this small problem neither training mode wins every time, and the figure shows how much one training run differs from the next. That is the honest state of the art. The remedies help on average and need tuning, and a paper that reports one seed and one horizon has not shown much.

How papers report it

Look for error as a function of rollout time, not only the error after one step. Common forms are the relative L2 error at several horizons, the VRMSE of The Well (the RMSE divided by the standard deviation of the true field, so 1.0 means no better than predicting the mean), and the time until the correlation between prediction and truth drops below a threshold. For chaotic flows a pointwise error must grow eventually; there, compare statistics such as spectra or mean fields instead.

Further on this site: MeshGraphNets (noise injection), UPT (rollout in latent space), The Well (VRMSE over long horizons) and PDE-Transformer (autoregressive and diffusion training on grids).

A U-Net maps one grid to another grid of the same size. A neural operator aims higher: it learns a map from one function to another: from the permeability field of a porous plate to its pressure field (Darcy flow, slow flow through a porous medium and the standard first benchmark), or from a wing's shape to its surface pressure. The functions are sampled on a grid or mesh for the computer, but the learned map is meant not to depend on how they are sampled.

Engineers know one family of operators already. For a linear PDE, the solution is the source integrated against a Green's function: u(x) = ∫ G(x, y) f(y) dy. A neural operator stacks layers of this form, with a learned kernel in place of G and a nonlinearity between layers. The two designs below differ in how they make that integral cheap.

DeepONet: a neural field with a branch

DeepONet (Lu, Karniadakis and colleagues) has two MLPs. The branch net reads the input function at m fixed sensor points and returns p numbers. The trunk net reads a query point y and also returns p numbers. The prediction at y is their dot product. The trunk is exactly a neural field; the branch decides how to weight its p outputs for this particular input.

import torch
from torch import nn

def mlp(i, o, h=128):
    return nn.Sequential(
        nn.Linear(i, h), nn.SiLU(),
        nn.Linear(h, h), nn.SiLU(), nn.Linear(h, o))

class DeepONet(nn.Module):
    def __init__(self, m_sensors, dim_y, p=64):
        super().__init__()
        # Branch: the input function at m fixed sensors.
        self.branch = mlp(m_sensors, p)
        # Trunk: the query point y.
        self.trunk = mlp(dim_y, p)
        self.bias = nn.Parameter(torch.zeros(1))

    def forward(self, u, y):         # u: (B, m)
        b = self.branch(u)[:, None]  # (B, 1, p)
        t = self.trunk(y)            # y: (B, Q, dim_y)
        # One value per query point: (B, Q)
        return (b * t).sum(-1) + self.bias

# 16 input functions, each sampled at 100 sensors,
# and 500 query points (x, y) per case.
u = torch.randn(16, 100)
y = torch.rand(16, 500, 2)
print(DeepONet(100, 2)(u, y).shape)  # (16, 500)

FNO: mixing in Fourier space

Any periodic signal is a sum of sines and cosines, its Fourier modes: mode 1 is one wave across the domain, mode 10 is ten waves. The fast Fourier transform (FFT) computes the modes from the samples, and the inverse FFT goes back. A Fourier Neural Operator layer takes the FFT of its input channels, keeps only the lowest k modes, multiplies them by learned complex weights (which mixes channels), and transforms back. A pointwise linear path is added alongside, then a nonlinearity. Because every kept mode spans the whole domain, one layer is global, unlike a 3×3 convolution.

The slider shows the cost of truncation. Smooth parts of the signal need few modes; the narrow peak needs many. An FNO with 16 modes per direction captures large-scale structure well and struggles with shocks and thin boundary layers, which the pointwise path and the nonlinearities have to rebuild.

import torch
from torch import nn

class SpectralConv1d(nn.Module):
    def __init__(self, c, k):
        super().__init__()
        self.k = k                     # Fourier modes kept
        w = torch.randn(c, c, k, dtype=torch.cfloat)
        self.w = nn.Parameter(w / c)

    def forward(self, x):              # x: (B, c, n)
        X = torch.fft.rfft(x)          # (B, c, n/2 + 1)
        Y = torch.zeros_like(X)
        # Mix channels mode by mode; drop the high modes.
        Y[..., :self.k] = torch.einsum(
            "bik,iok->bok", X[..., :self.k], self.w)
        return torch.fft.irfft(Y, n=x.shape[-1])

# In 2D and 3D (rfftn), the low modes sit at both ends of
# every axis except the last, so implementations keep two
# or four corner blocks of modes.
layer = SpectralConv1d(c=8, k=12)
print(layer(torch.randn(2, 8, 256)).shape)
print(layer(torch.randn(2, 8, 1024)).shape)

In practice, use the neuraloperator ↗ library. It adds a lifting layer (a pointwise MLP that widens the input channels to the hidden width) and a projection layer (the reverse, down to the output fields), pads non-periodic domains so the FFT does not wrap one edge onto the other, and ships tested training code. Its Darcy example ↗ trains the model below on a laptop CPU.

import torch
from neuralop.models import FNO

model = FNO(n_modes=(16, 16), hidden_channels=32,
            in_channels=1, out_channels=1)
# e.g. permeability on a 64 x 64 grid (Darcy flow)
a = torch.randn(4, 1, 64, 64)
print(model(a).shape)                  # (4, 1, 64, 64)
# The same weights run on a finer grid.
print(model(torch.randn(4, 1, 128, 128)).shape)

What discretisation-invariant means

An operator is called discretisation-invariant when the same weights can be applied to inputs sampled at different resolutions and give consistent answers. DeepONet's trunk accepts any query point; the FNO's weights live on Fourier modes, not on grid cells, so the last line above runs on a 128×128 grid. That is useful: train on coarse data, evaluate where you need.

GINO: operators on real 3D geometry

An FFT needs a regular grid, and a car is not one. GINO (Li and colleagues, 2023) wraps an FNO between two graph operators: the first aggregates the surface points onto a regular latent grid around the car (each grid point takes a learned, kernel-weighted sum over the surface points within a radius), the FNO works on that grid, and the second graph operator brings the result back to any query point on the surface. It was one of the first neural operators trained on full 3D car geometries for surface pressure.

Further on this site: DeepONet, FNO and GINO on Papers. The decoders lesson returns to DeepONet's trunk as the first query decoder.