Neural Operators and PDE Surrogates

A standard neural network maps vectors to vectors. Train it on 64×64 fluid simulations and it is a 4,096-dimensional function; ask it for 128×128 and it has nothing to say. A neural operator is trained to approximate a mapping between function spaces, so the same learned weights apply at any discretisation. That property, called discretisation invariance, is the reason the idea is more than a faster regressor.

The target is an operator \(\mathcal{G}: \mathcal{A} \to \mathcal{U}\) taking a coefficient or initial-condition function \(a\) to the PDE solution \(u\). Learn it once from data and you replace a solver call with a forward pass, which is worth doing when you need the solution for thousands of different inputs: inverse problems, design optimisation, uncertainty quantification, real-time control.

The Fourier neural operator

The key layer generalises a convolution to a kernel integral operator:

Evaluating this directly is quadratic in grid points. If the kernel is translation-invariant, \(\kappa(x,y) = \kappa(x-y)\), the integral is a convolution, and the convolution theorem turns it into a pointwise product in Fourier space. FNO does exactly that: FFT, multiply by a learned complex weight tensor on the lowest \(k\) modes, inverse FFT, add a pointwise linear term, apply a nonlinearity (Li et al., ICLR 2021, arXiv:2010.08895).

Three consequences follow. The layer is \(O(n \log n)\) rather than \(O(n^2)\). Truncating to \(k\) modes acts as a global low-pass filter, so every layer sees the whole domain, unlike a CNN's local receptive field. And because the parameters live in mode space rather than grid space, the same weights evaluate on a finer grid, which is where zero-shot super-resolution comes from.

The other main family, DeepONet, splits the problem instead: a branch network encodes the input function sampled at fixed sensor locations, a trunk network encodes the query coordinate, and the output is their inner product (Lu, Jin, Pang, Zhang and Karniadakis, Nature Machine Intelligence 3(3), 218–229, 2021). It rests on a universal approximation theorem for operators and handles arbitrary query points naturally, at the cost of a fixed sensor layout for the input.

How this differs from a PINN

Physics-informed neural networks put the PDE residual in the loss and fit a network to one solution instance (Raissi, Perdikaris and Karniadakis, Journal of Computational Physics 378, 2019). Change the boundary condition and you retrain.

A neural operator amortises across instances: expensive training once, then near-free inference for any input in the training distribution. The two are complementary rather than competing, and the physics-informed loss can be applied to an operator to reduce its data requirement.

Speedups, honestly stated

Reported inference speedups over classical solvers are large, often three orders of magnitude or more on the benchmark problems these papers use. Three caveats belong with every such number.

The comparison is usually against a general-purpose solver at the same resolution, not against a tuned production code with adaptive meshing. Training cost is excluded, and generating the training set requires running the solver thousands of times, so the break-even point is real. And the surrogate provides no error bound: a solver's error can be controlled by refining the mesh, a surrogate's cannot be controlled at all at inference time.

When it breaks

Distribution shift is silent and total. Evaluate on a Reynolds number outside the training range and the output is smooth, plausible, and wrong. There is no residual to check unless you compute one explicitly, which many pipelines never do.

Long rollouts drift. Autoregressive application of a learned operator accumulates error and can violate conservation laws that the training data satisfied exactly. Stability is not guaranteed by the architecture; it is coaxed with rollout training and noise injection.

Spectral bias against sharp features. Truncating to low Fourier modes is what makes FNO efficient and what makes it poor at shocks, discontinuities and boundary layers. Complex or non-periodic geometry needs variants (geometry-aware or graph-based operators), not the vanilla layer.

Benchmark saturation. Much of the literature reports on the same handful of problems (Burgers, Darcy flow, 2D Navier-Stokes at modest Reynolds number). Strong numbers there transfer poorly to industrial CFD, and the field has been explicit that harder, standardised benchmarks are needed.