Neural operators have emerged as powerful surrogates for dynamical systems due to their grid-invariant properties and computational efficiency. However, Fourier-based variants inherently truncate high-frequency components in spectral space, resulting in the loss of small-scale structures and degraded prediction quality at high resolutions when trained on low-resolution data. While diffusion-based enhancement methods can recover multi-scale features, they introduce substantial inference overhead that undermines the efficiency advantage of neural operators. In this work, we introduce MeanFlow-Enhanced Neural Operators (MENO), a novel framework that achieves accurate all-scale predictions with minimal inference cost. By leveraging the improved MeanFlow method, MENO restores both small-scale details and large-scale dynamics with superior physical fidelity and statistical accuracy. We evaluate MENO on three challenging dynamical systems, including phase-field dynamics, 2D Kolmogorov flow, and active matter dynamics, at resolutions up to 256×256. Across all benchmarks, MENO improves the power spectrum density accuracy by up to a factor of 2 compared to baseline neural operators while achieving up to 14× faster inference than the state-of-the-art Denoising Diffusion Implicit Model (DDIM)-enhanced counterparts, effectively bridging the gap between accuracy and efficiency. The flexibility and efficiency of MENO position it as an efficient surrogate model for scientific machine learning applications where both statistical integrity and computational efficiency are paramount.
Figures & tables
Figure 1: The MeanFlow-Enhanced Neural Operators framework. Panel (a) illustrates training the MeanFlow decoder on high-resolution fields by learning denoising trajectories. Panel (b) shows training an autoregressive neural operator on low-resolution data. Panel (c) combines both components into the full MENO pipeline: given a low-resolution initial condition, the neural operator produces a low-resolution rollout, which is then decoded into a high-resolution prediction by the MeanFlow decoder.
Figure 2: Relative L2 loss over rollout time (top row) and the corresponding wavenumber energy spectra (bottom row). Shaded regions denote the standard error of the mean (SEM) computed over multiple test trajectories (PF100: 10, KF256: 8, AM256: 25). (a) PF100 20→100 : a FNO model trained on 20×20 data and its enhanced variants. (b) KF256 64→256 : performance of a UNO model trained on 64×64 data and its enhanced variants. (c) AM256 64→256 : performance of a UNO model trained on 64×64 data and its enhanced variants.
20→100
50→100
Model
RL2↓
SSIM ↑
PSDD ↓
RL2↓
SSIM ↑
PSDD ↓
Neural operators without enhancement
FNO
0.64
0.60
4.98×10−5
0.26
0.85
1.99×10−5
Diffusion-enhanced neural operators
DM-FNO
0.50
0.65
2.97×10−5
0.30
0.77
2.28×10−5
DM-UNO
0.54
0.63
3.34×10−5
0.32
0.76
2.40×10−5
Table 1: PF100 prediction metrics for direct neural operator super-resolution, diffusion-enhanced neural operators, and MENO. The best value in each column is shown in bold. The stand-alone UNO super-resolution model was not included for PF100 because it did not produce reliable autoregressive rollouts under the same training setting as FNO. Relative improvement is computed with respect to the corresponding direct neural operator baseline.
Figure 3: Free-energy trajectories for all models and resolution settings on PF100. Upper panels: UNO-based models; lower panels: FNO-based models. Left column: 20→100 ; right column: 50→100 . Only the MENO variants consistently recover the correct free-energy evolution. Shaded regions denote the SEM over 10 test trajectories.
Figure 4: Temporal autocorrelation functions on PF100 across all resolution settings and model variants. Shaded bands denote the SEM over 10 test trajectories. MENO best preserves the long-time temporal coherence of the phase-field dynamics.
Figure 5: Snapshots along an uncurated PF100 trajectory for the 50→100 setting. The top row presents the high-resolution ground truth, and the next three rows visualize absolute errors for NO-SR, MENO-NO, and DM-NO, respectively. MENO produces the lowest interface-localized error and best preserves the evolving phase morphology.
Model
Parameters, NO + decoder
Inference time per frame (s)
Diffusion-enhanced neural operators
DM-FNO
6.33M + 1.44M
0.176±0.022
DM-UNO
6.38M + 1.44M
0.178±0.021
MENO models
MENO-FNO
6.33M + 1.44M
0.028±0.022
MENO-UNO
6.38M + 1.44M
0.031±0.021
Table 2: PF100 model size and inference time for generative-enhanced neural operators. Timings are measured per predicted frame for the full pipeline, including the neural operator and the decoder.
Figure 6: Snapshots along a trajectory from the KF256 dataset at t=6,12,…,54 . The top row reports the high-resolution ground truth at 256×256 . The next three rows visualize the absolute errors of UNO-SR, MENO-UNO, and DM-UNO, respectively. The color bars indicate the vorticity field and absolute error magnitude.
Figure 7: Temporal autocorrelation functions on KF256 across all resolution settings and model variants. Shaded bands denote the SEM over 8 test trajectories. The early-time decay is emphasized because it captures the dominant temporal memory of the turbulent vorticity field.
32→256
64→256
Model
RL2↓
SSIM ↑
PSDD ↓
RL2↓
SSIM ↑
PSDD ↓
Neural operators without enhancement
FNO
0.19
0.73
1.22×10−5
0.13
0.82
1.00×10−5
UNO
0.23
0.66
1.74×10−5
0.22
0.71
1.76×10−5
Diffusion-enhanced neural operators
DM-FNO
0.34
0.64
8.79×10−6
0.27
0.72
6.70×10−6
Table 3: KF256 prediction metrics for direct neural operator super-resolution, diffusion-enhanced neural operators, and MENO. RL2 and SSIM are computed over the first 20 frames, while PSDD is computed over the full 180-frame trajectories. The best value in each column is shown in bold. Relative improvement is computed with respect to the corresponding direct neural operator baseline.
Model
Parameters, NO + decoder
Inference time per frame (s)
Diffusion-enhanced neural operators
DM-FNO
27.10M + 6.63M
0.377±0.008
DM-UNO
26.83M + 6.63M
0.376±0.008
MENO models
MENO-FNO
27.10M + 6.63M
0.030±0.004
MENO-UNO
26.83M + 6.63M
0.030±0.004
Table 4: KF256 model size and inference time for generative-enhanced neural operators. Timings are measured per predicted frame for the full pipeline, including the neural operator and the decoder.
Figure 8: Snapshots along an uncurated AM256 trajectory for the 64→256 setting. The top row presents the high-resolution ground truth, and the next three rows visualize absolute errors for UNO-SR, MENO-UNO, and DM-UNO, respectively. Frames are normalized per frame to highlight error evolution since the mean concentration drifts over time.
Figure 9: Temporal autocorrelation functions on AM256 across all resolution settings and model variants. Shaded bands denote the SEM over 25 test trajectories. The early-time decay is emphasized because it captures the dominant temporal coherence of the active concentration field.
Model
Parameters, NO + decoder
Inference time per frame (s)
Diffusion-enhanced neural operators
DM-UNO
26.83M + 6.63M
0.456±0.011
MENO models
MENO-UNO
26.83M + 6.63M
0.032±0.005
Table 5: AM256 model size and inference time for generative-enhanced neural operators. Timings are measured per predicted frame for the full pipeline, including the neural operator and the decoder.
32→256
64→256
Model
RL2↓
SSIM ↑
PSDD ↓
RL2↓
SSIM ↑
PSDD ↓
Neural operators without enhancement
UNO
2.11×10−3
0.72
1.09×10−5
2.05×10−3
0.73
1.12×10−5
Diffusion-enhanced neural operators
DM-UNO
2.17×10−3
0.56
1.24×10−5
2.14×10−3
0.56
1.30×10−5
MENO models
Table 6: AM256 prediction metrics for direct neural operator super-resolution, diffusion-enhanced neural operators, and MENO. RL2 and SSIM are computed over the first 10 frames, while PSDD is computed over the full 80-frame trajectories. The best value in each column is shown in bold. Relative improvement is computed with respect to the corresponding direct neural operator baseline.
Appendix figures & tables10 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 10: Image adapted from [ 32 ] , showing the reconstruction error trends for the 32→256 generative refinement task. The legend compares three reconstruction strategies; in this work, we focus on the Baseline method.
Figure 11: Absolute L2 loss of MENO decoder models on the KF256 dataset for the 32→256 and 64→256 generative refinement tasks.
Figure 12: Long range relative L2 loss of all UNO-based models on the KF256 dataset for the 32→256 generative refinement tasks. Shaded regions denote the SEM computed over 8 test trajectories.
Model
RL2↓
SSIM ↑
PSDD ↓
i-MF variants
i-MF:MENO-FNO
0.21
0.76
8.62×10−6
i-MF:MENO-UNO
0.21
0.76
9.06×10−6
MF variants
MF:MENO-FNO
0.23
0.75
8.96×10−6
MF:MENO-UNO
0.25
0.74
9.21×10−6
Appendix
Table 7: Ablation study on the KF256 dataset comparing i-MF- and MF-based variants of MENO. All results are reported for the 32→256 setting.
Figure 13: Training loss curves of i-MF and the original MF on the KF256 dataset, shown with moving-average smoothing. The loss of the original MeanFlow increases as training proceeds, consistent with the behavior reported in Figure 3 of [ 13 ] .
32 → 256
64 → 256
Model
RL2↓
SSIM ↑
PSDD ↓
RL2↓
SSIM ↑
PSDD ↓
DM-UNO
0.3541(4)
0.6383(2)
(9.221±0.001)×10−6
0.2581(1)
0.7178(3)
(5.764±0.001)×10−6
MENO-UNO
0.2131(1)
0.7612(1)
(9.058±0.001)×10−6
0.0813(1)
0.9321(2)
(5.161±0.001)×10−6
Appendix
Table 8: KF256 32→256 metrics for DM-enhanced NOs, and MENOs. RL2 and SSIM are computed over the first 20 frames, while PSDD is computed over the full trajectories (180 frames). Uncertainties are computed over 100 runs with different random seeds. We report values to four significant figures; uncertainties smaller than this precision are rounded and reported as ±1 in the last digit.
Parameter Name
FNO Configuration
UNO Configuration
channel_mlp_dropout
0
0
channel_mlp_expansion
0.5
0.5
channel_mlp_skip
soft-gating
linear
factorization
tucker
tucker
fno_skip
linear
linear
hidden_channels
78
128
Appendix
Table 9: FNO and UNO architectures used in all tasks for PF100 dataset.
Long-term fluid dynamics forecasting is a critically important problem in science and engineering. While neural operators have emerged as a promising paradigm for modeling systems governed by partial differential equations (PDEs), they often struggle with long-term stability and precision. We identify two fundamental failure modes in existing architectures: (1) local detail blurring, where fine-scale structures such as vortex cores and sharp gradients are progressively smoothed, and (2) global trend deviation, where the overall motion trajectory drifts from the ground truth during extended rollouts. We argue that these failures arise because existing neural operators treat local and global information processing uniformly, despite their inherently different evolution characteristics in physical systems. To bridge this gap, we propose the Dual-Scale Neural Operator (DSO), which explicitly decouples information processing into two complementary modules: depthwise separable convolutions for fine-grained local feature extraction and an MLP-Mixer for long-range global aggregation. Through numerical experiments on vortex dynamics, we demonstrate that nearby perturbations primarily affect local vortex structure while distant perturbations influence global motion trends, providing empirical validation for our design choice. Extensive experiments on turbulent flow benchmarks show that DSO achieves state-of-the-art accuracy while maintaining robust long-term stability, reducing prediction error by over 88% compared to existing neural operators.
Huanshuo Dong, Hao Wu, Hong Wang +2
University of Science and Technology of China, Hefei, China · Tsinghua University, Beijing, China · Institute of Automation, Chinese Academy of Sciences, Beijing, China +1
Fast surrogate modeling for high-dimensional physical dynamics requires more than low short-term error: useful models must roll out efficiently while preserving the statistical structure of long trajectories. Neural operators provide inexpensive autoregressive forecasts but can drift in turbulent regimes, whereas rolling diffusion and latent generative surrogates can represent stochastic transitions at the cost of multi-step denoising, noise-schedule design, or auxiliary compression models. We propose MeanFlow Long-term Invariant Spatiotemporal Consistency Autoregressive Models (MeLISA), a latent-free autoregressive generative surrogate built on pixel-space MeanFlow. MeLISA defines a blockwise stochastic transition kernel that generates each forecast block with a single model evaluation, avoiding latent encoders and iterative diffusion solvers at inference time. To stabilize long-horizon rollouts, MeLISA combines a Window-Consistency MeanFlow objective that learns conditional spatiotemporal generation from partially observed temporal windows with a Time Increment Consistency loss that constrains multi-lag finite increments and targets temporal-correlation structure. We evaluate MeLISA with compact UNet and scalable DiT backbones on two high-resolution benchmarks, extended 2D Kolmogorov flow at 256×256 and turbulent channel-flow slice at 192×192. MeLISA outperforms neural-operator baselines on short-term forecasting accuracy and long-horizon statistical metrics, including energy spectra, turbulent kinetic energy, and mixing-rate-related dynamics, while achieving inference speeds comparable to, and in some cases faster than, neural operators. Compact 3.7-5.7M-parameter variants already deliver strong parameter efficiency, and DiT variants provide a scalable path up to 150M parameters. Overall, MeLISA benefits both rollout efficiency and long-horizon statistical accuracy.
Tianyue Yang, Xiao Xue
The Center for Computational Science University College London
Neural operators enable fast PDE forecasting, but repeated predictions accumulate errors and fine-scale structures remain difficult to resolve. We introduce a frequency-decomposed finite-time flow-map neural operator (F3NO) that leverages updated low-frequency features to guide nonlinear refinement of high-frequency information. Within each layer, this cross-scale conditioning connects global spectral processing with local detail refinement. The model directly predicts states at specified future times and adjusts the contributions of the two branches according to the prediction interval. For longer trajectories, it combines parallel predictions within short temporal segments with recursive propagation between segments. Experiments on five PDE benchmarks demonstrate improved forecasting accuracy over autoregressive and direct-prediction baselines. Ablations show that frequency-decomposed refinement can improve accuracy with fewer parameters, while the benefits of segmentation depend on spatial resolution and dynamical regime.