The Well, a 15-terabyte dataset spanning 16 physical simulation domains, provides an extensive benchmark for developing and evaluating machine learning surrogate models for complex physical systems. The baselines released with the benchmark were presented by their authors as illustrative examples obtained at small computational cost, trained with zero to very limited hyperparameter tuning and under a highly constraining 12-hour per-run wall-clock limit. These results have subsequently been treated as fully tuned benchmarks, and in some cases effectively interpreted as unconstrained performance baselines. We perform a short non-exhaustive hyperparameter search for the Tucker-factorized Fourier Neural Operator (TFNO) under the same 12-hour per-run wall-clock limit, while using only a single input timestep rather than the history of up to four used by the published baselines. Despite this more restrictive setup, tuning alone reduces VRMSE below both published FNO and TFNO baselines on 15 of 17 tasks, sometimes by several orders of magnitude. This sensitivity to a modest tuning effort underscores that the released neural operator baselines should not be treated as architectural benchmarks.
Figures & tables
Dataset
12h-tuned
vs. pub. TFNO
vs. pub. FNO
vs. Best Reported
acoustic_scattering_maze
0.022886
22.1×
22.1×
0.7×
active_matter
0.077859
4.6×
4.7×
1.3×
convective_envelope_rsg
0.030917
0.9×
0.9×
0.9×
euler_multi_quadrants
0.060337
6.9×
6.8×
2.5×
gray_scott_reaction_diffusion
0.022986
15.8×
5.9×
5.9×
helmholtz_staircase
0.000280
12.5×
1.8×
1.8×
Table 1: Sensitivity of the reported operator baselines to protocol and configuration. Column 2 gives the VRMSE of our 12h-tuned TFNO. Columns 3 and 4 give the ratio of the published TFNO and FNO VRMSEs from The Well [ 4 ] to ours, each computed as published value divided by tuned value ; column 5 gives the same ratio against the best (lowest) of all four published baselines on that dataset. Values above one indicate that our run scores lower than the reference; values below one indicate that it scores higher. The two baselines are reported separately rather than aggregated, since on gray_scott_reaction_diffusion and helmholtz_staircase the published FNO substantially outperforms the published TFNO, and any single aggregate would obscure this.
Configuration and outcome
Ratios
Budget
Modes
Width
Batch
LR
Rank
Epochs
VRMSE
vs. 12h
vs. operator baseline
12 h
64
32
32
2×10−3
0.20
17
0.022886
1.0×
22.1×
36 h
128
48
8
8×10−4
1.00
38
0.004875
4.7×
103.8×
Table 2: Effect of increasing the per-run wall-clock budget from 12 to 36 hours on the dataset acoustic_scattering_maze . For each budget, the configuration and cyclic learning-rate schedule are re-selected, so the comparison captures the combined effect of additional training time and the different configurations that become competitive under the larger budget; it is not a continuation of the 12-hour run. Both runs are halted by their respective wall clocks while validation error is still decreasing and should therefore be treated as budget-limited. vs. 12h is 0.022886 divided by the row’s VRMSE. vs. operator baseline uses the larger of the published FNO and TFNO VRMSEs ( 0.5062 ) divided by the row’s VRMSE. Ratios against the published operator values quantify how far those protocol-specific figures can be displaced and are not evidence of architectural superiority.
Appendix figures & tables24 assets
Supplementary material from the paper’s appendix.
Appendix
Hyperparameter
Candidate values
Layers
4, 6
Spectral modes
32, 64, 128, 256
Channel width
24, 30, 32, 40, 48, 50, 64, 96
Batch size
1–64, set by memory limits per dataset
Learning rate
1×10−4 , 5×10−4 , 1×10−3 , 2×10−3
Tucker rank
0.05, 0.10, 0.20, 0.50, 1.00
Appendix
Table 3: Candidate values per varied hyperparameter. All other settings were fixed to the values in the caption of Table 4 . The number of spectral modes was chosen separately for each dataset based on resolution.
Dataset
Layers
Modes
Width
Batch Size
LR
Rank
Patience
Epochs
ES
acoustic_scattering_maze
4
64
32
32
2×10−3
0.20
8
17
×
active_matter
6
128
48
32
1×10−3
0.20
4
171
✓
convective_envelope_rsg
6
64
24
1
5×10−4
0.05
6
40
✓
euler_multi_quadrants
6
128
30
12
5×10−4
1.00
4
4
×
gray_scott_reaction_diffusion
6
64
64
32
5×10−4
0.10
8
7
×
helmholtz_staircase
6
128
64
8
5×10−4
0.10
4
64
✓
Appendix
Table 4: Selected TFNO configurations and training outcomes for the 17 reported 12-hour runs. All runs use model=FNO , factorization=Tucker , optimizer=Adam , lr_method=cycle , normalize=1 , lifting_channels=projection_channels=64 , and max_train_hours=12 . Epochs is the number completed before termination. ES marks whether the early-stopping criterion fired within the budget (✓) or the run was halted by the 12-hour wall clock while validation error was still decreasing ( × ). A ✓ records only the behavior of the budget-shaped schedule and is not evidence of budget-independent convergence (Section 3.2 ). A rank of 1.00 denotes the dense limit. Configurations were selected by the limited search described in Section 2.5 and are not claimed to be optimal.
Figure 1: acoustic_scattering_maze , pressure p : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 2: active_matter , concentration c : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 3: convective_envelope_rsg , radial velocity vr : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 4: euler_multi_quadrants , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 5: gray_scott_reaction_diffusion , species A : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 6: helmholtz_staircase , real part of the pressure, Rep : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 7: MHD_64 , magnetic field Bx : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 8: planetswe , height h : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 9: post_neutron_star_merger , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 10: rayleigh_benard , buoyancy b : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 11: rayleigh_taylor_instability , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 12: rayleigh_taylor_instability , velocity vx : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 13: rayleigh_taylor_instability , velocity vy : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 14: rayleigh_taylor_instability , velocity vz : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 15: shear_flow , passive tracer c : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 16: supernova_explosion_64 , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 17: turbulence_gravity_cooling , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 18: turbulence_gravity_cooling , pressure p : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 19: turbulence_gravity_cooling , temperature T : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 20: turbulent_radiative_layer_2D , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 21: turbulent_radiative_layer_3D , density ρ : ground truth (top) and tuned-TFNO prediction (bottom).
Figure 22: viscoelastic_instability , conformation tensor czz : ground truth (top) and tuned-TFNO prediction (bottom).
We develop a comprehensive mathematical and computational framework for sequential surrogate modeling of three-phase black-oil reservoir dynamics using neural operators, with particular emphasis on Fourier Neural Operators (FNO) and their physics-informed variant (PINO). The application focus is the Norne benchmark reservoir, defined on a heterogeneous 46×112×22 grid (N=113,344 cells), with a production history spanning T=30 timesteps covering 3298 days. Our theoretical contributions are organized around four interlocking problems: (1) functional-analytic formulation in a product-Sobolev-space setting, including well-posedness of the implicit timestep map and sharp local Lipschitz estimates; (2) covariate shift quantification, proving that the Wasserstein-2 distance grows as W2≤ε(Ln−1)/(L−1), with exponential population-risk discrepancy for L>1; (3) physics-constrained spectral stability, showing PINO training with λR≥λR∗ reduces the learned Jacobian spectral radius to ρF+CλR−1/2, yielding uniform-in-time rollout error ∣δn∣≤ε/(1−ρ); and (4) K-step TBPTT gradient analysis, deriving geometric bias decay O(ρK), optimal window K=O(log(T/σ2)), and Adam convergence O(1/t)+O(ρK∗). Empirical validation confirms all theoretical predictions: autoregressive PINO surrogates sustain R2>0.99 (oil), R2>0.90 (gas), R2≈0.80 (pressure), and monotonically improving R2 (water) across the full 3298-day horizon, trained on eight NVIDIA B200 GPUs in under one hour. A 1000-member ensemble runs in under one minute on a single B200 GPU, giving a ∼104× wall-clock speedup over the OPM finite-volume simulator.
Neural operators (NOs) are designed to learn maps between infinite-dimensional function spaces. We propose a novel reframing of their use. By introducing an auxiliary base-space, any finite-dimensional function can be viewed as an operator acting by composition on functions of the base-space. Through a range of benchmarks on analytic functions of increasing complexity and dimensionality, we demonstrate that NOs can match or outperform standard multilayer perceptrons and Kolmogorov--Arnold Networks in accuracy while requiring significantly fewer parameters and training time. As a real-world application, we apply a two-dimensional Tensorized Fourier Neural Operator (TFNO) to the nuclear chart, learning a correction to state-of-the-art nuclear mass models as a partially observed residual field. A TFNO ensemble reaches a held-out root-mean-square error of 198.2 keV, placing it among the best recent neural-network approaches while retaining high parameter efficiency and short training times. More broadly, these results introduce NOs as a scalable framework for finite-dimensional function interpolation, from analytic benchmarks to structured scientific data.
Department of Physics, CCTP and ITCP, University of Crete, 71303, Greece · Theoretical Physics Department, CERN, Geneva, Switzerland · Physik-Institut, Universität Zürich, 8057 Zürich, Switzerland
Neural operators are widely used as surrogate solution maps for partial differential equations (PDEs), but full-size models can be costly to store, deploy, and evaluate in many-query scientific workflows. This work introduces Operator Boosting, a stagewise residual-learning framework for constructing compact neural-operator surrogates directly, rather than training a large model and compressing it afterward. Starting from the empirical mean predictor in normalized output coordinates, the method trains a sequence of tiny same-family neural operators on residual fields and incorporates each correction through validation-selected shrinkage. We instantiate the framework with Fourier neural operators (FNOs), DeepONets, and convolutional neural operators (CNOs), and compare boosted tiny stacks against full-size monolithic baselines across one-, two-, and three-dimensional PDE benchmarks from PDEBench, APEBench, and The Well. Across 30 dataset-architecture pairs, 21 show positive mean accuracy gains and 17 have positive confidence intervals, while all boosted stacks reduce trainable parameter count by approximately 72-95%. Best-model comparisons show empirical Pareto improvements on 7 of 10 completed PDE benchmarks, including two-dimensional Navier-Stokes, shallow-water dynamics, Darcy flow, one-dimensional transport and reaction systems, and three-dimensional compressible Navier-Stokes. These results show that Operator Boosting often improves the empirical accuracy-parameter Pareto frontier of neural PDE surrogates, while also exposing PDE- and architecture-dependent regimes where residual boosting fails to offset compression.
Lennon J. Shikhman
Georgia Institute of Technology College of Computing Atlanta, Georgia, USA