Neural Operators vs PINN: FNO and DeepONet as Engineering Surrogate Models (2026)

Most engineers meet physics-informed neural networks first, because the pitch is irresistible: write the governing equation into the loss function and let the network discover the solution without a mesh. Then they try to use one inside a digital twin and hit a wall. A physics-informed network solves one instance of one problem. Change the inflow velocity, the material coefficient or the geometry, and you retrain. A twin that must answer “what happens if” thousands of times per hour cannot retrain for each question.

Neural operators attack exactly that gap. Instead of learning one solution, they learn the map from inputs (coefficients, forcing terms, initial conditions) to solution fields, so a trained model answers a new query in a single forward pass. This article is a practitioner’s comparison of neural operators vs PINN approaches: the maths of the Fourier neural operator and DeepONet, the failure modes of PINNs, the data you actually need, and when a classical reduced-order model still wins.

What this covers: operator learning versus function approximation, FNO spectral convolution, DeepONet branch and trunk networks, PINN loss design and failure modes, resolution invariance, minimal PyTorch code, tooling, validation and uncertainty, and a decision framework for digital twin surrogates.

Context and Background

Engineering simulation spends most of its compute on repeated solves. A finite element or finite volume solver takes a geometry, boundary conditions and material parameters, then marches to a solution. Design exploration, calibration, optimisation and real-time monitoring all want that solve repeated many times with slightly different inputs. A surrogate model replaces the expensive solver with a cheap approximation that is fast enough to sit inside a loop or a live twin.

The classical surrogate toolbox is well established. Response surfaces, Gaussian processes and polynomial chaos handle low-dimensional parameter spaces. Projection-based reduced-order models, built from proper orthogonal decomposition (POD) and Galerkin projection, handle high-dimensional fields when the solution manifold is low-rank. We covered that family, including the PINN connection, in our guide to reduced-order models and surrogate models for real-time digital twins. This post picks up where it stops and focuses on the two deep learning families that have changed the conversation since 2019.

The first family is the physics-informed neural network (PINN), popularised by Raissi, Perdikaris and Karniadakis, which trains a network to satisfy a PDE residual at sampled points. The second is operator learning: networks that approximate mappings between function spaces. The DeepONet paper by Lu, Jin and Karniadakis grounded this in a universal approximation result for operators (DeepONet, arXiv:1910.03193), and the Fourier neural operator paper by Li and colleagues showed that parameterising the integral kernel in Fourier space gave an efficient, expressive architecture. In their abstract, the FNO authors report it as up to three orders of magnitude faster than traditional PDE solvers on the benchmarks they studied, with zero-shot super-resolution on turbulent flows (Fourier Neural Operator, arXiv:2010.08895).

Treat those speed-up claims with discipline. They are measured on specific benchmark PDEs (Burgers, Darcy flow, Navier-Stokes on periodic domains) against specific baseline solvers, usually after the training cost has been paid. For an industrial geometry with complex boundaries and multi-physics coupling, the honest statement is that neural operators are fast at inference once trained, and that the training data bill and the accuracy envelope decide whether they pay off. The rest of this article is about reading that bill.

A final piece of context: the term “surrogate model for a digital twin” hides two distinct jobs. One is offline surrogate-assisted design, where latency is not critical and accuracy over a wide design space matters. The other is online state estimation and what-if analysis, where latency, calibration against sensors and graceful degradation matter more than global accuracy. Neural operators and PINNs are not equally good at both, and much of the confusion in the literature comes from mixing the two.

Operator Learning vs Function Approximation: The Core Distinction

Short answer: A PINN approximates a single function, the solution u(x, t) of one PDE instance, by minimising a physics residual. A neural operator approximates an operator G that maps an input function a (a coefficient field, forcing term or initial condition) to the solution u = G(a). Train a neural operator once on many input-solution pairs and it solves new instances without retraining.

Figure 1: Two workflows. A PINN trains per problem instance against a PDE residual; a neural operator trains once on a dataset of input-output function pairs and then answers new queries by forward pass.

Figure 1 shows the structural difference. On the left, a PINN takes coordinates as input and returns the field value at those coordinates. Its weights encode one solution. On the right, a neural operator takes a discretised input function and returns a discretised output function, and its weights encode the solution operator for a whole family of problems. The first is function approximation; the second is function-space regression.

Why the distinction matters for digital twins

Consider a heat conduction problem on a plate, where the unknown is temperature and the conductivity field k(x) varies from unit to unit. A PINN for a specific k(x) needs retraining whenever the conductivity map changes. A neural operator trained on a few thousand conductivity fields and their solved temperature fields can predict the temperature for a new k(x) in milliseconds. In a twin that ingests as-built variation for each physical asset, that difference is the whole business case.

The trade is data. A PINN can train with no solution data at all, using only the PDE and boundary conditions, which is why it is attractive when simulations are unavailable. A neural operator is usually trained supervised, which means someone has already run the solver many times. Operator learning therefore amortises solver cost: you pay N solver runs up front and then save on every later query. The break-even count is simple arithmetic, and we return to it in the practical section.

The operator view of the maths

Formally, let A and U be function spaces, for example continuous functions on a domain D. The solution operator G: A to U sends a coefficient field a to the solution u of a PDE such as the Darcy flow equation -div(a grad u) = f. The goal is to learn an approximation G_theta from samples (a_i, u_i). Two properties are desired. First, discretisation convergence: as the grid gets finer, the predicted function should converge to the same underlying function instead of a different one. Second, mesh independence: a model trained on one resolution should evaluate sensibly at another.

Ordinary convolutional networks and fully connected networks tie their parameters to a grid size. A U-Net trained on 64 by 64 images sees a different problem at 128 by 128. Neural operators are built so that the learned parameters define an operator on functions, with the grid entering only as a sampling of those functions. That is the conceptual leap, and it is the origin of the resolution-invariance claims we examine later.

The universal approximation backbone

DeepONet’s theoretical motivation is a result due to Chen and Chen from the 1990s: a network with a single hidden layer can approximate a broad class of continuous nonlinear operators. DeepONet turns that statement into an architecture with two subnetworks, a branch net that encodes the input function sampled at fixed sensor locations and a trunk net that encodes the location where the output is evaluated. The authors demonstrated it on explicit operators such as integrals and fractional Laplacians, and on implicit operators defined through differential equations, both deterministic and stochastic. Approximation theory tells you that a good operator exists; it does not tell you how much data or training time you need to find one. That gap is why empirical validation, not theory, drives engineering decisions.

How the Fourier Neural Operator Works

The Fourier neural operator builds on an old idea. Many linear PDE solution operators are convolutions with a Green’s function. A convolution is diagonal in Fourier space: multiply the transform of the input by the transform of the kernel, then invert. FNO takes this literally and learns the kernel’s action directly in the frequency domain.

The spectral convolution layer

An FNO has three stages: a pointwise lifting layer P that raises the input channels to a wider hidden width, a stack of Fourier layers, and a pointwise projection Q back to the output channels. Each Fourier layer computes

v_{t+1}(x) = sigma( W v_t(x) + F^{-1}( R_phi . F(v_t) )(x) )

where F is the fast Fourier transform, R_phi is a learned complex weight tensor applied to a truncated set of low-frequency modes, W is a pointwise linear map (a 1 by 1 convolution acting as a bypass for local and high-frequency content), and sigma is a nonlinearity such as GELU. The truncation is essential: only the lowest k_max modes get learned weights, and higher modes are zeroed in the spectral branch. That makes the parameter count independent of grid resolution, and it acts as a built-in low-pass regulariser.

Figure 2: One Fourier layer. The spectral branch transforms, keeps the lowest k modes, multiplies by learned complex weights, and transforms back; the pointwise bypass preserves local detail. Layers are stacked between a lifting and a projection map.

Figure 2 lays out the data path. The cost of a layer is dominated by the FFT, which is O(N log N) in the number of grid points N, versus O(N^2) for a dense global kernel. That scaling is why FNO handles global interactions cheaply and why it ran on grids that would be prohibitive for dense kernel integration.

Resolution invariance, honestly stated

Because the learned weights live on a fixed set of Fourier modes, an FNO trained on a 64 by 64 grid can be evaluated on a 256 by 256 grid by feeding it the finer input and letting the FFT handle the new resolution. The FNO paper calls this zero-shot super-resolution. It works well when the solution’s energy is concentrated in the retained low modes and the underlying physics is smooth at the scale of the training grid.

It works less well when the finer grid exposes physics the training grid never resolved. A shock, a thin boundary layer or a turbulent cascade at 256 by 256 contains information that was aliased or absent at 64 by 64, and no truncated spectral layer can invent it. The correct reading is that FNO is discretisation-convergent for functions whose content is resolved, not that it recovers unresolved physics. Practitioners should test across resolutions on held-out data instead of assuming the invariance holds for their problem.

Limits of the plain FNO

The FFT assumes a regular grid and, in its basic form, periodic boundary conditions. Industrial meshes are unstructured, geometries are irregular and boundaries are rarely periodic. The community has responded with variants: geometry-aware FNO flavours that map irregular domains to latent regular grids, graph and geometry-informed neural operators for point clouds and meshes, and spherical FNO for global climate-style problems. The open-source NeuralOperator library lists several of these (FNO, TFNO, SFNO, GNO, GINO, OTNO, LocalNO, RNO and PINO) in its model zoo, and describes itself as MIT licensed and part of the PyTorch ecosystem (NeuralOperator library paper, arXiv:2412.10354). The breadth of that list is itself a signal: the plain FNO is a starting point, and real engineering geometry usually pulls you toward a variant.

Boundary handling is the most common practical failure. Non-periodic problems are often handled by padding the input before the FFT, which reduces wrap-around artefacts but does not eliminate them. If your twin deals with walls, inlets and outlets, inspect the error near the boundaries separately from the interior; global relative L2 error can hide a boundary layer that is the thing you care about.

DeepONet and PINNs: Architecture and Failure Modes

DeepONet: branch, trunk and a dot product

DeepONet represents the operator as a finite sum. For an input function a sampled at m fixed sensor points, and an output query location y, the prediction is

G(a)(y) approx sum over k of b_k(a(x_1), …, a(x_m)) * t_k(y) + bias

The branch network produces coefficients b_k from the sensor values, the trunk network produces basis functions t_k(y) from the query coordinate, and the output is their inner product. If you have met POD, the structure will look familiar: the trunk learns a nonlinear basis in place of a precomputed linear one, and the branch learns the coefficients that a Galerkin projection would otherwise compute.

Figure 3: DeepONet. The branch net encodes the input function from sensor values, the trunk net encodes the query coordinate, and their inner product gives the output value at that coordinate.

Two practical consequences follow from Figure 3. First, DeepONet can be queried at arbitrary output coordinates, because the trunk network accepts any y. That makes it comfortable with scattered, unstructured output points, which is a genuine advantage over a grid-bound FNO when your sensors or mesh nodes are irregular. Second, the input is tied to a fixed sensor layout. The branch net expects the same m locations every time. If your input sensors change between samples, you need interpolation or a different architecture.

Variants relax these limits. Physics-informed DeepONet adds a PDE residual term to the loss, MIONet handles multiple input functions, and the POD-DeepONet family replaces the trunk with a precomputed POD basis for stability. Because the lineage runs straight from the Karniadakis group’s earlier PINN work to DeepONet, the two are often combined, which we return to under physics-informed operators.

PINNs: the loss function is the model

A PINN represents the solution u_theta(x, t) with a network and trains it by minimising a weighted sum:

where N is the differential operator, evaluated at collocation points x_f using automatic differentiation, B is the boundary operator, and the last term is optional sensor data. No mesh is needed, derivatives are exact up to floating point, and inverse problems (estimating unknown coefficients from sparse data) drop out naturally by making the coefficient a trainable parameter. That last capability is why PINNs remain valuable for calibration in digital twins even when a neural operator is the forward model.

Why PINNs fail on engineering problems

The honest literature on PINNs is a catalogue of failure modes. Krishnapriyan and colleagues at Berkeley studied convection, reaction and diffusion problems and found that PINNs can fail even on moderately complex cases. They argued that the soft-constraint form of the physics loss makes the loss landscape hard to optimise, rather than the network lacking expressivity. Their remedies were curriculum regularisation, which starts with an easy version of the PDE and gradually increases difficulty, and sequence-to-sequence training over time windows, and they reported error reductions of up to one to two orders of magnitude over vanilla training on their test problems (Characterizing possible failure modes in PINNs, arXiv:2109.01050). Those figures are from their benchmarks, not a guarantee.

Several recurring mechanisms explain why:

Loss imbalance. The residual, boundary and data terms have different magnitudes and gradient scales. Without adaptive weighting, one term dominates and the network satisfies it while violating the others. Gradient-statistics and neural-tangent-kernel based weighting schemes exist, but they add hyperparameters.
Stiffness and multiscale behaviour. Equations with widely separated time or length scales, or high-wavenumber solutions, suffer from spectral bias: standard networks learn low frequencies first and struggle with high ones.
Causality violation. Time-dependent PINNs trained on the whole space-time domain at once may fit late times before early times are right. Causal training schemes weight early times first.
Poorly conditioned residuals. The PDE residual at a point is a differential expression of the network. Its minimum may not correspond to the true solution when the problem is ill-posed or the network has plenty of freedom, and a small residual is not the same as a small solution error.
Cost. Each training step needs higher-order automatic differentiation. For second-order PDEs on 3D domains, training can cost far more per step than a plain supervised network, and there is no accuracy guarantee that classical solvers give by construction.

None of this makes PINNs useless. For smooth problems, inverse problems, data-assimilation settings and cases where no simulator exists, they are valuable. The point is narrower: a PINN is a per-instance solver with a tricky optimiser, and for a repeated forward-solve workload it is usually the wrong tool.

Physics-informed operators: combining the two

PINO, the physics-informed neural operator, adds a PDE-residual term to FNO training and has been used to reduce the data needed, as the NeuralOperator library notes by including it. The appeal is clear: use a handful of simulation snapshots to anchor the model and the equation to regularise it. The caveat is that the residual term brings back the optimisation difficulties of PINNs, and computing the residual on a spectral model typically uses Fourier-based derivatives, which assume smooth, periodic fields. Treat physics-informed operators as a data-efficiency lever to test, not a default.

Neural Operators vs PINN: Side-by-Side Decision Matrix

The following matrix summarises the dimensions that matter for engineering use. It reflects the published properties of each method, not a ranking of quality; “best” depends on which row matters to you.

Dimension	PINN	DeepONet	FNO	POD-Galerkin ROM
What it learns	One solution instance	Operator, arbitrary output points	Operator on a grid	Linear basis plus reduced equations
Needs simulation data	Optional	Yes (supervised)	Yes (supervised)	Yes (snapshots)
New input at inference	Retrain	Forward pass	Forward pass	Reduced solve
Mesh handling	Meshless	Flexible outputs, fixed input sensors	Regular grid, variants for meshes	Tied to full-order mesh
Resolution transfer	N/A	Query any coordinate	Zero-shot on resolved physics	Fixed by snapshots
Training difficulty	High (loss balancing)	Moderate	Moderate	Low (SVD)
Inverse problems	Natural	Possible	Possible via gradients	Possible
Guarantees	None built in	None built in	None built in	Stability and error bounds possible
Best fit	Inverse and small smooth problems	Scattered output queries	Smooth fields on grids	Low-rank, parametric, well-understood physics

Read the last row against the first. A POD-Galerkin ROM is the only column offering structural guarantees, which matters for safety-relevant twins. Neural operators win on flexibility and on problems where the reduced-basis assumption (a low-dimensional linear manifold) breaks, such as advection-dominated flows with moving features.

Hands-On: Minimal PyTorch Implementations

The code below is deliberately small and self-contained so you can run it on a laptop. It implements a 1D spectral convolution and FNO, a DeepONet, and a PINN residual for the 1D viscous Burgers equation. These are teaching implementations; for production, use the NeuralOperator library or NVIDIA PhysicsNeMo.

A 1D Fourier neural operator

import torch
import torch.nn as nn
import torch.nn.functional as F


class SpectralConv1d(nn.Module):
    """Multiply the lowest `modes` Fourier coefficients by learned complex weights."""

    def __init__(self, in_ch: int, out_ch: int, modes: int):
        super().__init__()
        self.modes = modes
        scale = 1.0 / (in_ch * out_ch)
        self.weights = nn.Parameter(
            scale * torch.randn(in_ch, out_ch, modes, dtype=torch.cfloat)
        )

    def forward(self, x):                      # x: (batch, in_ch, n)
        n = x.shape[-1]
        x_ft = torch.fft.rfft(x, dim=-1)       # (batch, in_ch, n//2 + 1)
        out_ft = torch.zeros(
            x.shape[0], self.weights.shape[1], x_ft.shape[-1],
            dtype=torch.cfloat, device=x.device,
        )
        k = min(self.modes, x_ft.shape[-1])
        out_ft[:, :, :k] = torch.einsum(
            "bik,iok->bok", x_ft[:, :, :k], self.weights[:, :, :k]
        )
        return torch.fft.irfft(out_ft, n=n, dim=-1)


class FNO1d(nn.Module):
    def __init__(self, modes=16, width=64, layers=4):
        super().__init__()
        self.lift = nn.Conv1d(2, width, 1)     # inputs: a(x) and the grid x
        self.spec = nn.ModuleList(SpectralConv1d(width, width, modes) for _ in range(layers))
        self.skip = nn.ModuleList(nn.Conv1d(width, width, 1) for _ in range(layers))
        self.proj = nn.Sequential(nn.Conv1d(width, 128, 1), nn.GELU(), nn.Conv1d(128, 1, 1))

    def forward(self, a):                      # a: (batch, n)
        n = a.shape[-1]
        grid = torch.linspace(0, 1, n, device=a.device).expand_as(a)
        v = self.lift(torch.stack([a, grid], dim=1))
        for spec, skip in zip(self.spec, self.skip):
            v = F.gelu(spec(v) + skip(v))
        return self.proj(v).squeeze(1)


def rel_l2(pred, true):
    return (torch.linalg.vector_norm(pred - true, dim=-1)
            / torch.linalg.vector_norm(true, dim=-1)).mean()


if __name__ == "__main__":
    torch.manual_seed(0)
    model = FNO1d()
    # Toy operator: u = integral of a (cumulative sum), a smooth random input family.
    n = 128
    x = torch.linspace(0, 1, n)
    def sample(b):
        c = torch.randn(b, 4)
        a = sum(c[:, [k]] * torch.sin((k + 1) * torch.pi * x) for k in range(4))
        return a, torch.cumsum(a, dim=-1) / n
    opt = torch.optim.Adam(model.parameters(), lr=1e-3)
    for step in range(500):
        a, u = sample(32)
        loss = rel_l2(model(a), u)
        opt.zero_grad(); loss.backward(); opt.step()
        if step % 100 == 0:
            print(step, float(loss))

Two details deserve attention. The grid coordinate is concatenated to the input so the network can learn position-dependent behaviour, a standard FNO trick. And the weight tensor is complex, so the optimiser updates real and imaginary parts; PyTorch handles this natively with torch.cfloat. Because the number of retained modes is fixed at 16, you can evaluate this model on n=256 inputs without changing a parameter, which is how you would test resolution transfer on your own data.

A DeepONet

import torch
import torch.nn as nn


def mlp(sizes):
    layers = []
    for i in range(len(sizes) - 1):
        layers.append(nn.Linear(sizes[i], sizes[i + 1]))
        if i < len(sizes) - 2:
            layers.append(nn.Tanh())
    return nn.Sequential(*layers)


class DeepONet(nn.Module):
    def __init__(self, n_sensors=100, coord_dim=1, p=64, hidden=128):
        super().__init__()
        self.branch = mlp([n_sensors, hidden, hidden, p])
        self.trunk = mlp([coord_dim, hidden, hidden, p])
        self.bias = nn.Parameter(torch.zeros(1))

    def forward(self, a_sensors, y):
        # a_sensors: (batch, n_sensors); y: (n_points, coord_dim)
        b = self.branch(a_sensors)             # (batch, p)
        t = self.trunk(y)                      # (n_points, p)
        return b @ t.T + self.bias             # (batch, n_points)


if __name__ == "__main__":
    torch.manual_seed(0)
    m = 100
    xs = torch.linspace(0, 1, m)
    net = DeepONet(n_sensors=m)
    opt = torch.optim.Adam(net.parameters(), lr=1e-3)
    y = torch.linspace(0, 1, 200).unsqueeze(-1)
    for step in range(2000):
        c = torch.randn(64, 3)
        a = sum(c[:, [k]] * torch.sin((k + 1) * torch.pi * xs) for k in range(3))
        u = torch.cumsum(
            sum(c[:, [k]] * torch.sin((k + 1) * torch.pi * y.squeeze(-1)) for k in range(3)),
            dim=-1) / 200
        loss = ((net(a, y) - u) ** 2).mean()
        opt.zero_grad(); loss.backward(); opt.step()
        if step % 400 == 0:
            print(step, float(loss))

The key line is b @ t.T: the inner product between branch coefficients and trunk basis functions, evaluated for every batch item and every query point at once. Notice that y can be any set of coordinates at inference time, including points that never appeared in training, which is the structural reason DeepONet suits scattered outputs.

A PINN residual for viscous Burgers

import torch
import torch.nn as nn

nu = 0.01 / torch.pi
net = nn.Sequential(nn.Linear(2, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(),
                    nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1))


def residual(xt):
    xt = xt.requires_grad_(True)               # columns: x, t
    u = net(xt)
    g = torch.autograd.grad(u.sum(), xt, create_graph=True)[0]
    u_x, u_t = g[:, :1], g[:, 1:]
    u_xx = torch.autograd.grad(u_x.sum(), xt, create_graph=True)[0][:, :1]
    return u_t + u * u_x - nu * u_xx           # Burgers: u_t + u u_x = nu u_xx


def loss_fn(n_f=4000, n_b=200, w_pde=1.0, w_bc=10.0, w_ic=10.0):
    xt = torch.rand(n_f, 2) * torch.tensor([2.0, 1.0]) - torch.tensor([1.0, 0.0])
    l_pde = (residual(xt) ** 2).mean()
    t_b = torch.rand(n_b, 1)
    l_bc = (net(torch.cat([-torch.ones_like(t_b), t_b], 1)) ** 2).mean() \
         + (net(torch.cat([torch.ones_like(t_b), t_b], 1)) ** 2).mean()
    x_i = torch.rand(n_b, 1) * 2 - 1
    l_ic = ((net(torch.cat([x_i, torch.zeros_like(x_i)], 1)) + torch.sin(torch.pi * x_i)) ** 2).mean()
    return w_pde * l_pde + w_bc * l_bc + w_ic * l_ic


opt = torch.optim.Adam(net.parameters(), lr=1e-3)
for step in range(5000):
    loss = loss_fn()
    opt.zero_grad(); loss.backward(); opt.step()
    if step % 1000 == 0:
        print(step, float(loss))

The hand-picked weights w_bc=10 and w_ic=10 are exactly the loss-balancing knobs that cause trouble at scale. At this viscosity the problem forms a steep gradient near x = 0, which is where a plain PINN needs the most tuning. Change nu to a smaller value and watch the optimiser struggle: that is the stiffness failure mode in miniature. Notice also that this network has to be retrained for any new initial condition, while the two operator models above are trained once over a family of inputs.

Data, Cost and Deployment: Making a Neural Operator Pay Off

The amortisation arithmetic

The economics of a neural operator reduce to one inequality. Let C_sim be the cost of one high-fidelity solve, N_train the number of training simulations, C_train the training cost, C_inf the inference cost and Q the number of queries you expect to serve. The surrogate wins when N_train * C_sim + C_train + Q * C_inf is below Q * C_sim. With illustrative numbers, suppose a solve takes 30 minutes on a node, you need 2,000 training runs, and training takes a few GPU-hours. The data bill is about 1,000 node-hours. If the twin will answer 100,000 what-if queries over its life, the direct solver would need 50,000 node-hours, so the operator pays back many times over. If you only need 500 queries, it never does. These figures are illustrative, not measured; substitute your own.

This is why the hardest question in a neural operator project is not architecture but query volume. Surrogates are for loops: optimisation, calibration with many iterations, uncertainty propagation through thousands of Monte Carlo samples, and live twins that re-evaluate on every sensor update. A one-off study is better served by running the solver.

How much data is enough

There is no universal number. The FNO benchmarks used on the order of a thousand training pairs for 2D problems like Darcy flow, but your geometry, parameter dimension and target accuracy decide the real figure. Error typically falls as a power law in the training set size until it hits a floor set by model capacity, input noise or the solver’s own discretisation error. The sensible workflow is a learning curve: train on 125, 250, 500 and 1,000 samples, plot held-out error, and extrapolate before committing to a data-generation campaign. Sample the input space with a space-filling design (Latin hypercube or Sobol sequences) rather than random draws, since extrapolation outside the sampled region is where neural operators fail without warning.

Tooling in 2026

You do not need to write spectral layers yourself. Two open-source options dominate. The NeuralOperator library offers the FNO family, DeepONet-style and graph-based models, a trainer, data modules for Darcy flow, Burgers, Navier-Stokes and a 3D car-surface flow dataset, mixed-precision training and tensorised spectral weights to cut parameter counts. It is MIT licensed per its authors (NeuralOperator library, arXiv:2412.10354). NVIDIA PhysicsNeMo, the framework previously branded Modulus, documents FNOs, graph networks, diffusion models, transformer-style operator models and weather and climate models among its modules (NVIDIA PhysicsNeMo documentation). Its home page does not explicitly list DeepONet, so check the current examples catalogue before assuming it. NVIDIA positions it for building digital twins together with Omniverse, which matters if your twin already lives in that ecosystem.

Where this sits in a digital twin

A surrogate is only useful if it plugs into the twin’s data flow. Figure 4 shows a pattern that works in practice: simulation campaigns generate the training set offline, the trained operator is served behind an inference endpoint, the twin’s state estimator calls it for forward predictions, and a monitoring loop compares predictions with sensor streams to trigger fallback to the full solver or retraining.

Figure 4: Deployment pattern for a surrogate model in a digital twin. Offline simulations train the operator; at runtime an uncertainty or drift gate decides whether to trust the surrogate or fall back to the high-fidelity solver.

The gate in Figure 4 is the part teams skip, and it is the part that keeps the twin honest. Without it, the surrogate returns a confident-looking field for inputs far outside the training distribution. Add a cheap out-of-distribution check on the input parameters (distance to the training design, or an ensemble variance) and route flagged queries to the real solver. The same discipline applies to asset data: the structure and provenance of the inputs feeding the surrogate matter as much as the model, a theme explored in our piece on ISO 19650 and BIM data models for construction digital twins, where information requirements and identifiers keep model inputs traceable.

Validating and quantifying uncertainty

Report error in a way that matches the decision. Relative L2 error over the full field is the common headline, but also report error on quantities of interest (peak temperature, pressure drop, maximum stress), on boundary regions, and on a held-out set that is structurally different from training, such as a new geometry family or an extended parameter range. A model with 1% average field error can still miss a hot spot by 20%.

For uncertainty, the pragmatic options are deep ensembles (train five to ten models from different seeds and use the spread), Monte Carlo dropout, and conformal prediction on quantities of interest, which provides distribution-free coverage under exchangeability. Ensembles are expensive but robust. Bayesian variants of DeepONet exist in the literature, yet they add tuning burden. Whatever you pick, calibrate against held-out simulations and check that the claimed intervals actually cover the truth at the stated rate.

Trade-offs, Gotchas, and What Goes Wrong

Extrapolation is silent failure. A neural operator interpolates within its training distribution and degrades unpredictably outside it. Unlike a ROM built on a physical model, it has no structural reason to respect conservation laws when the inputs shift. Mitigate with the OOD gate above and with physics-based checks on outputs, such as global mass or energy balance, that flag nonsense without needing ground truth.

Conservation and constraints are not automatic. Standard FNO and DeepONet outputs do not enforce divergence-free velocity fields, positivity of density or energy balance. Hard constraints can be built in through architecture (predicting a stream function so velocity is divergence-free) or through projection post-processing. If your downstream logic assumes a conserved quantity, test it.

Long rollouts accumulate error. For time-dependent problems, autoregressive stepping compounds small errors, and a model that is accurate one step ahead can drift badly after hundreds of steps. Techniques such as noise injection during training, multi-step loss and periodic correction against sensors help, but you should measure rollout error, not single-step error.

Boundary and geometry handling. The plain FNO assumes a regular grid. Mapping an irregular CAD geometry to a grid introduces its own errors, and geometry-aware variants add complexity. If geometry varies a lot, budget for it, or consider a mesh-based graph surrogate.

High-frequency content. Truncated spectral layers and smooth MLPs both underfit sharp features. Shocks, contact discontinuities and thin layers are the hard cases, and they are often the engineering-critical ones.

Training data inherits solver error. The surrogate cannot be more accurate than the solver that made its training data. If the CFD case is mesh-dependent or the turbulence model is wrong for the regime, the surrogate faithfully reproduces that wrongness.

PINNs: over-promising. A frequent anti-pattern is choosing a PINN to avoid generating data, then spending weeks tuning loss weights and collocation sampling on a problem a coarse solver would have handled in an afternoon. If a simulator exists, use it to generate data for an operator or ROM before reaching for a PINN.

When ROM beats both. If your problem is parametric, the solution manifold is low-rank, and the physics is diffusion or elliptic-dominated, POD-Galerkin or a data-driven POD-plus-regression ROM is cheaper to build, needs less data, is easier to certify and often more accurate. Neural operators earn their keep when the linear-subspace assumption breaks (transport, moving fronts, strong nonlinearity) or when you want a single model across wide parameter ranges.

Practical Recommendations

Start from the use case, not the architecture. If you have a working simulator and need many fast evaluations, generate data and train a neural operator or a ROM. If you have sparse measurements and an unknown coefficient to estimate, a PINN, or a physics-informed fit of a trusted forward model, is the better fit. If outputs are needed at scattered coordinates and inputs come from a fixed sensor layout, DeepONet is a natural choice. If your fields live on a regular grid and are smooth, start with FNO.

Always benchmark against a baseline you already trust. A POD-based ROM, a Gaussian process on reduced coordinates, or even a well-tuned U-Net gives you a reality check on whether the neural operator adds value. In published comparisons, results vary by problem, so avoid broad claims that one method dominates.

Use this checklist before committing:

Estimate query volume and compute the amortisation break-even with your own solver cost.
Build a learning curve with at least four training set sizes before scaling data generation.
Define quantities of interest and acceptance thresholds up front, not only field-wide L2 error.
Hold out structurally different test cases to measure extrapolation.
Implement an input OOD gate and a solver fallback path.
Add physics sanity checks (conservation, positivity, boundary conditions) on every prediction.
Use an ensemble or conformal wrapper for uncertainty on decision-critical outputs.
Version the training data, solver settings and model together so the twin can be audited.
Monitor drift against sensor data in production and schedule retraining.

Frequently Asked Questions

What is the difference between neural operators and PINNs?

A PINN approximates one solution of one PDE instance by minimising the equation residual at sampled points, so a new coefficient or boundary condition means retraining. A neural operator learns the mapping from input functions to solution functions from many examples, so a new instance needs only a forward pass. PINNs can train without simulation data; neural operators usually need supervised pairs, which makes them an amortisation of solver cost rather than a replacement for it.

Is the Fourier neural operator faster than a CFD solver?

At inference, often dramatically. The FNO authors report up to three orders of magnitude over traditional solvers on their benchmarks, which were periodic-domain problems such as Burgers, Darcy and Navier-Stokes. That figure excludes the cost of generating training data and does not transfer automatically to complex industrial geometries. Compare total cost over your expected query count, and validate accuracy on quantities of interest rather than relying on a headline speed-up.

When should I use DeepONet instead of FNO?

Choose DeepONet when you need outputs at arbitrary or scattered coordinates, when your inputs come from a fixed set of sensors, or when you want a basis-function structure similar to a ROM. Choose FNO when fields sit on regular grids, are smooth and benefit from global spectral mixing. For irregular meshes, look at geometry-aware or graph-based neural operators. Benchmark both on your data, because relative performance is problem dependent.

Why do PINNs fail on some PDEs?

Research on PINN failure modes points to optimisation difficulty from soft-constraint losses, which creates ill-conditioned landscapes on convection-dominated, reaction or multiscale problems. Other causes include unbalanced loss terms, spectral bias towards low frequencies and violated temporal causality. Curriculum regularisation and sequence-to-sequence training reduced error by one to two orders of magnitude in one study, but they add design effort. Many practitioners therefore reserve PINNs for inverse and smooth problems.

Can neural operators replace reduced-order models in digital twins?

Sometimes, but not by default. POD-Galerkin ROMs are cheap to build, need little data and can offer error bounds, which suits low-rank, well-understood physics. Neural operators help when the solution manifold is not low-dimensional, as with transport-dominated flows, or when one model must span wide parameter ranges. Many teams run both and choose per subsystem; see our guide to ROMs and surrogates for real-time twins for the projection-based side.

Which libraries should I use for neural operators in 2026?

The NeuralOperator library is an MIT licensed PyTorch package with FNO variants, trainers and benchmark datasets. NVIDIA PhysicsNeMo provides FNO, graph and other physics-AI modules with GPU-oriented tooling and Omniverse integration for twins. For learning, the short implementations in this article show the mechanics. Always check each project’s current documentation for supported models and versions before standardising.