Neural solvers offer efficient surrogates for numerical simulation of partial differential equations (PDEs). For time-dependent problems, strong one-step accuracy does not necessarily translate into reliable autoregressive rollout. We observe that a solver based only on physical-state modeling can achieve lower one-step error, whereas its spectral-only counterpart can become more accurate at later rollout steps. Motivated by this observation, we present Transolver-σ, a neural PDE solver based on joint spectral--physical subspace modeling. Within each block, adaptive physical-state interactions and spectral transformations are modeled in dedicated latent subspaces, whose responses are recomposed to enable information exchange between the two representations. Within the physical subspace, we introduce Slice-Residual Physics-Attention (SRPA), which preserves an explicit slice-space identity path while retaining learnable cross-slice interaction. In parallel, an axis-factorized Fourier operator captures global spectral structure. Across five well-established PDE benchmarks spanning steady-state prediction and time-dependent dynamics, Transolver-σ achieves state-of-the-art with a benchmark-averaged relative error reduction of 33.4% over the strongest baseline for each metric, while consistently improving autoregressive rollout over single-operator counterparts. Transolver-σ further delivers strong gains on coupled multiphysics systems and real-world fluid and combustion measurements from RealPDEBench, demonstrating its effectiveness beyond standard simulation benchmarks.
Figures & tables
Figure 1: Single-step accuracy does not guarantee rollout robustness. (a) Error dynamics under teacher-forced prediction (ground-truth inputs) versus autoregressive rollout (recursive inputs) on 2D Navier–Stokes. Physolver achieves lower single-step error, whereas Specsolver achieves lower error at later rollout steps. (b) Conceptual overview of Transolver- σ , which leverages joint spectral-physical modeling via learned recomposition to achieve a balanced error trade-off.
Figure 2: Transolver- σ architecture and spectral–physical modeling. (a) Transolver- σ performs specialized physical and spectral modeling in decoupled latent subspaces, followed by learned recomposition. (b) The physical branch captures intrinsic state interactions through slice-based modeling, while the spectral branch captures global spatial variations along individual axes.
Figure 3: Comparison of slice-based operator updates. (a) Transolver updates physical states with solely attention. (b) LinearNO uses linear global-context aggregation. (c) SRPA retains slice states with scaled residual attention.
Models
Darcy
Navier– Stokes
2D Kolmogorov
3D isotropic
3D smoke
p
ω
ωavg.
ωfinal
u
p
u
d
U-Net ( 2015 )
0.80
19.82
32.81
47.24
35.76
48.49
44.10
14.95
ViT ( 2021 )
0.47
4.64
17.19
27.19
31.66
44.54
38.41
12.86
FNO ( 2021 )
1.08
15.56
29.78
45.67
33.82
46.34
42.55
13.44
F-FNO ( 2023 )
0.77
23.22
24.53
38.61
23.03
32.64
37.13
12.36
RNO ( 2025 )
0.54
8.94
18.93
29.14
22.37
30.72
28.51
11.53
Table 1: Results on Standard Benchmarks. Relative L2 errors (in %) of pressure p , vorticity ω , velocity u , and smoke density d across standard benchmarks, where “avg.” and “final” denote horizon-aggregated and final-frame errors. Our results are reported as mean ± std over three runs.
Models
Params (M)
Controlled Cylinder
FSI
Foil
Combustion
RMSE
Rel L2
fRMSE
RMSE
Rel L2
fRMSE
RMSE
Rel L2
fRMSE
RMSE
Rel L2
fRMSE
U-Net
23.08
0.80
5.55
0.10
0.85
5.83
0.07
1.00
1.59
0.11
2.16
54.87
0.26
CNO
8.00
0.81
5.83
0.09
1.05
7.41
0.09
1.36
2.53
0.18
2.48
60.30
0.32
DeepONet
3.53
3.09
23.99
0.58
3.50
25.02
0.51
2.22
3.63
0.28
2.29
57.51
0.28
FNO
109.10
0.97
7.23
0.12
1.29
8.92
0.12
1.30
2.28
0.15
2.26
56.64
0.27
WDNO
233.15
1.15
9.27
0.14
1.17
8.17
0.11
1.62
4.44
0.18
3.80
86.05
0.55
Table 2: Results on RealPDEBench. RMSE, relative L2 , and Fourier-space RMSE (fRMSE) across four real-world measurement tasks. Params denotes the mean parameter count (M) across tasks.
Figure 4: Autoregressive forecasting errors across four dynamical systems. (a–d) Relative L2 errors over prediction leads on Standard Benchmarks, with n the evaluation sample count. (e) Percentage reduction in rollout-averaged error over the stronger single-domain baseline.
Variant
Darcy
NS-2D
Smoke-3D
p
ω
u
d
w/o Physical Branch
0.65
4.96
26.29
8.97
w/o Spectral Branch
0.40
7.94
29.02
10.15
Physical → Spectral
0.42
3.04
20.75
7.78
Spectral → Physical
0.42
3.15
21.88
8.06
w/o Shared CPE
0.42
3.69
23.56
8.59
Table 3: Ablations of Transolver- σ . Relative L2 errors (%) for pressure p , vorticity ω , velocity u , and density d . Branch ablations use identity mappings; “w/o Cross-Branch FFN” removes cross-branch mixing while retaining grouped SwiGLU. SRPA → PA replaces SRPA with Physics-Attention.
Figure 5: Training efficiency and learned representations. (a,b) Epoch time and peak GPU memory on NS2D ( 642 ) and Kolmogorov flow ( 1282 ), with bubble area indicating parameter storage and error bars denoting ±1 standard deviation over three runs. (c, d) Spatial routing maps and slice-attention matrices extracted from the final layer on NS2D.
Figure 6: Case study on error maps of different models. Within each task, models share the same test case, prediction time, spatial view, and error scale. See Appendix D.2 for more details.
Model
Slice-space update
Darcy
NS2D
LinearNO
–
0.5033
6.9871
Transolver
AZWv
0.5807
9.0315
– w/o Slice Attention
Z
0.4724
7.7124
– with SRPA
Z+γℓAZWv
0.4275
6.6735
Table 4: Effects of slice-attention updates.
Appendix figures & tables17 assets
Supplementary material from the paper’s appendix.
Appendix
Benchmark
Data type
Spatial dim.
Grid
Input
Output
Canonical PDE benchmarks
Darcy
Simulation
2D
852
Coefficient a
Pressure p
Navier–Stokes
Simulation
2D
642
ω history
Vorticity ω
Kolmogorov flow
Simulation
2D
1282
ω history
Vorticity ω
Isotropic turbulence
Simulation
3D
603
(u,p) history
u,p
Smoke buoyancy
Simulation
3D
643
(u,d) history
u,d
Appendix
Table 5: Overview of the evaluated benchmarks. Dimensions are spatial; canonical forecasting configurations are summarized in Table 6 .
Task
Temporal type
Grid
Target fields
Train
Test
Tin
Tout
H
Darcy
Steady-state
852
p
1,000
200
–
–
–
Navier–Stokes
Time-dependent
642
ω
1,000
200
10
1
10
Kolmogorov flow
Time-dependent
1282
ω
100
20
10
4
16
Isotropic turbulence
Time-dependent
603
u,p
1,000
100
10
2
10
Smoke buoyancy
Time-dependent
643
u,d
2,000
200
4
4
16
Appendix
Table 6: Canonical PDE benchmarks and forecasting configurations. Train and test sizes count source sequences for time-dependent tasks and individual samples for Darcy, rather than extracted windows. Tin , Tout , and H denote the input history length, frames predicted per model call, and total forecasting horizon, respectively, all measured in frames. p , ω , u , and d denote pressure, vorticity, velocity, and smoke density.
Benchmark
Training Configuration
Model Configuration
Budget
Batch
Peak LR
Optimizer
Width
Depth
Heads
Slices
Modes
C
L
h
M
K
Canonical PDE benchmarks
Darcy
500
4
10−3
AdamW
128
8
8
32
4
NS2D
500
2
10−3
AdamW
256
8
8
32
4
KF2D
50
8
10−3
AdamW
128
8
8
32
8
Appendix
Table 7: Training and model configurations of Transolver- σ . Canonical training configurations like the AdamW optimizer ( Loshchilov and Hutter, 2019 ) are shared by the baselines trained in this study. C , L , h , and M denote latent width, depth, attention heads, and physical slices, respectively. The spectral and physical subspaces each use half of C . K denotes the retained Fourier modes.
Task
Grid
L
C
nh
M
K
Darcy
852
8
128
8
32
4
NS2D
642
8
256
8
32
4
KF2D
1282
8
128
8
32
8
IT3D
603
8
128
8
32
6
Smoke3D
643
8
128
8
32
6
Appendix
Table 8: Checkpoint configurations for layer-wise SRPA analysis. L , C , nh , M , and K denote block count, total hidden width, attention heads, slices, and Fourier modes per axis, respectively. IT3D denotes isotropic turbulence.
Figure 7: Learned residual scales and realized SRPA corrections. Columns show eight-layer profiles for Darcy, NS2D, KF2D, isotropic turbulence (IT3D), and Smoke3D. Top: solid curves show signed γℓ ; pale dashed curves show ∣γℓ∣ , with shading between the signed curve and zero. Bottom: correction ratios in slice-token space before deslicing and output projection, expressed as percentages. Thin curves show three fixed held-out cases; thick curves show their means, with shading spanning their range. Horizontal dashed lines indicate means across all cases and layers. Ratios use joint norms over heads, slices, and channels before case averaging. Layers are one-based; vertical scales vary across tasks. Neither shaded region is a confidence interval.
Figure 8: Autoregressive error maps on NS2D. Columns show leads 1 – 10 for the same test case. Rows show Physolver, Specsolver, and Transolver- σ from top to bottom, using a shared error scale. Physolver develops pronounced error bands at later leads; Transolver- σ retains weaker spatial errors across the displayed sequence.
Figure 9: Autoregressive error maps on KF2D. Columns show leads 1 – 16 for the same test case. Rows show Physolver, Specsolver, and Transolver- σ from top to bottom, using a shared error scale. Errors are more prominent in both single-domain models than in Transolver- σ .
Figure 10: Autoregressive error maps on isotropic turbulence (IT3D). Columns show leads 1 – 10 for the same test case. Rows show Physolver, Specsolver, and Transolver- σ from top to bottom, using a shared error scale. The single-domain models develop more extensive high-error regions at later leads than Transolver- σ .
Figure 11: Autoregressive error maps on Smoke3D. Columns show leads 1 – 16 for the same test case. Rows (top to bottom): Physolver, Specsolver, and Transolver- σ , all sharing one error scale. Transolver- σ has weaker localized errors at later leads; spatial patterns evolve non-monotonically.
Models
IgnitHIT
EvolveJet
Corr
train
val
test
Corr
train
val
test
CNext
96.09
2.47
3.17
2.71
92.56
1.06
4.25
2.59
CROP
94.91
0.20
5.24
4.37
84.95
1.66
3.70
3.41
DPOT
94.19
1.47
6.15
5.90
87.71
2.10
5.32
3.36
DeepONet
68.60
39.63
48.54
48.25
62.20
7.85
28.79
18.33
F-FNO
97.36
0.52
1.86
1.87
95.08
0.54
2.50
0.98
Appendix
Table 9: Results on REALM. Correlation ( ×100 , higher is better) and training, validation, and test errors (lower is better) on IgnitHIT and EvolveJet. Errors follow the original benchmark scale, and baseline results are taken from the official comparison. Our results are averaged over three training runs evaluated at their final checkpoints; standard deviations are provided in Appendix E.2 . Best and second-best results are bold and underlined, respectively.
Figure 12: Scalability on Darcy flow. Relative L2 error as a function of (a) training-set size, (b) grid resolution, and (c) hidden width. Each point is the mean over the same 200 test cases; error bars show 1.96s/200 , where s is the standard deviation of the per-case errors. The legends use S for slices and M for Fourier modes, corresponding to M and K in the manuscript.
Data scaling: 852 grid, C=128 , M=32 , K=4
Training cases
1,000
2,000
3,000
4,000
5,000
Error
3.912±1.259
3.397±1.022
3.304±0.943
3.195±0.859
3.143±0.853
Resolution scaling: 1,000 training cases, C=128 , M=32 , K=4
Grid
852
1062
1412
2112
4212
Error
3.912±1.259
3.394±1.281
2.876±1.042
2.745±1.043
2.699±0.812
Parameter scaling: 3,000 training cases, 2112 grid, M=128 , K=8
Appendix
Table 10: Darcy scaling results. Relative L2 errors ( ×10−3 ) are mean ± standard deviation across 200 test cases. Data and resolution sweeps use 200 epochs; width scaling uses 100 .
Model
Darcy
NS2D
KF2D
IT3D
Smoke3D
p
ω
Avg.
Final
u
p
u
d
Strongest baseline
0.45
4.23
12.69
21.99
12.70
19.43
23.37
9.19
Transolver- σ
0.40±0.01
2.79±0.05
2.52±0.05
4.14±0.07
10.80±0.16
15.65±0.23
17.43±0.25
7.07±0.09
Appendix
Table 11: Standard deviations on the canonical benchmarks. Relative L2 errors (%) are reported as mean ± standard deviation over three training runs. Strongest baseline values are taken from Table 1 : DRIFT-Net for Darcy, FactFormer for NS2D and Smoke3D, and EddyFormer for KF2D and IT3D.
Metric
Model
Controlled Cylinder
FSI
Foil
Combustion
RMSE
Strongest baseline
0.80
0.85
1.00
2.08
Transolver- σ
0.650±0.010
0.743±0.015
0.470±0.010
1.860±0.040
Relative L2
Strongest baseline
5.55
5.83
1.59
53.31
Transolver- σ
5.200±0.080
5.550±0.080
1.403±0.025
52.857±0.550
fRMSE
Strongest baseline
0.09
0.07
0.10
0.24
Transolver- σ
0.060±0.002
0.060±0.002
0.080±0.002
0.180±0.004
Appendix
Table 12: Standard deviations on RealPDEBench. Our entries show mean ± standard deviation across three training runs; all metrics retain the 100× display scale of the main comparison. Strongest baseline values are taken from Table 2 .
Dataset
Model
Correlation
Train error
Val. error
Test error
IgnitHIT
Strongest baseline
97.36
0.20
1.86
1.87
Transolver- σ
98.493±0.267
0.470±0.010
1.503±0.035
1.493±0.035
EvolveJet
Strongest baseline
95.08
0.54
2.15
0.98
Transolver- σ
95.930±0.255
0.030±0.001
1.310±0.030
0.630±0.020
Appendix
Table 13: Standard deviations on REALM. Our entries show mean ± standard deviation over three training runs. Correlations are multiplied by 100 ; errors retain the benchmark scale. Strongest baseline values are taken from Table 9 .
Task
n
Tout
H
NS2D
200
1
10
KF2D
400
4
16
IT3D
100
2
10
Smoke3D
200
4
16
Appendix
Table 14: Evaluation settings for the exact error decomposition. n counts trajectories (windows for KF2D); Tout denotes output frames per call and H the total horizon.
Model
NS2D (L,C,h)
Kolmogorov (L,C,h)
Operator settings
HPM
(8,256,8)
(8,128,8)
128 fixed frequencies
DRIFT-Net
Official T hierarchy
Patch 4, window 16
FactFormer
(8,256,8)
(8,128,8)
Head dimension 64
Transolver
(8,256,8)
(8,128,8)
M=32
Transolver- σ
(8,256,8)
(8,128,8)
M=32K=4 (NS2D), K=8 (KF2D)
Appendix
Table 15: Model configurations for training-efficiency measurements. L , C , and h denote backbone depth, width, and head count; M and K denote slice count and retained Fourier modes per axis. Model-specific settings preserve the respective native operator structures.
Neural surrogate models offer fast approximate mappings from PDE parameters to solutions, but they typically treat solving as a purely statistical task: once trained, they struggle to correct their own constraint violations and extrapolate beyond the training distribution. Recent hybrid methods promote physical correctness by targeting the PDE residual via gradient descent or Gauss--Newton steps, but inherit the compute cost and instability of the underlying classical optimizers. We show, theoretically and empirically, that numerically minimizing the PDE residual can be an unreliable proxy for reconstruction accuracy in ill-conditioned systems, explaining why these methods often do not make accurate predictions despite achieving low residuals. We propose error-conditioned Neural Solvers (ENS), built on a different principle: rather than an optimization target, the PDE residual field is passed as a direct input to the network at each iteration, enabling it to read the spatial structure of its own errors and learn an update policy to iteratively correct its predictions. Across four PDE families, ENS attains the highest prediction accuracy in the large majority of settings, with gains reaching 10× on turbulent Kolmogorov flow, while avoiding the expensive compute cost of hybrid methods. ENS's learned correction policy generalizes under distribution shift, including zero-shot parameter changes and cross-equation transfer, where its relative advantage is largest in the ill-conditioned regimes where residual minimization is least reliable. Project website: https://neuralsolver.github.io/.
Haina Jiang, Liam Wang, Peng-Chen Chen +4
University of Michigan · KAIST AI · Los Alamos National Laboratory
Deep learning has emerged as a transformative tool for the neural surrogate modeling of partial differential equations (PDEs), known as neural PDE solvers. However, scaling these solvers to industrial-scale geometries with over 108 cells remains a fundamental challenge due to the prohibitive memory complexity of processing high-resolution meshes. We present Transolver-3, a new member of the Transolver family as a highly scalable framework designed for high-fidelity physics simulations. To bridge the gap between limited GPU capacity and the resolution requirements of complex engineering tasks, we introduce two key architectural optimizations: faster slice and deslice by exploiting matrix multiplication associative property and geometry slice tiling to partition the computation of physical states. Combined with an amortized training strategy by learning on random subsets of original high-resolution meshes and a physical state caching technique during inference, Transolver-3 enables high-fidelity field prediction on industrial-scale meshes. Extensive experiments demonstrate that Transolver-3 can handle meshes with over 160 million cells, achieving impressive performance across three challenging simulation benchmarks, including aircraft and automotive design tasks. Code is available at https://github.com/thuml/Transolver-3.
Hang Zhou, Haixu Wu, Haonan Shangguan +4
School of Software, BNRist, Tsinghua University, China.
Neural spectral PDE solvers often learn an entire unresolved vector field even when an inexpensive approximate model can already capture most of the trajectory. Here we introduce Perturbative-NeuSA, a residual formulation that decomposes the target solution into a low-fidelity background and a high-resolution perturbation, so that only the unresolved dynamics is learned. Starting from the exact perturbation equation, the method combines a fixed spectral operator, a background-dependent correction, the background defect in the target PDE, and an optional neural closure. This construction makes the roles of physical structure and neural closure separately measurable. Across 2D Burgers, Klein-Gordon, and heterogeneous 2D wave equations, the deterministic structured solver outperforms the trained NeuSA baseline while requiring no neural-network training. The largest gains occur on Burgers, where the deterministic correction reduces training and extrapolation errors by factors of 24 and 44, respectively. In addition, a Klein-Gordon sweep over seven background resolutions shows that the effect of the closure is conditional: it improves a poor background by 3.6 times, becomes neutral at intermediate resolutions, and degrades a well-resolved background. For the wave equation, however, the closure provides an additional 18% reduction when the remaining residual is interface-localized. Multi-initial-condition diagnostics further show that the useful closure regime depends on the initial-condition spectrum and can disappear in extrapolation when structured correction already captures the dominant Burgers dynamics. Perturbative-NeuSA therefore reframes neural closure as a conditional, diagnosable correction governed by background fidelity, residual organization, and compatibility with the closure model.
Xianli Zhu, Jia Yin
School of Mathematical Sciences, Fudan University Shanghai, China