Learning Data-Efficient and Generalizable Neural Operators via Fundamental Physics Knowledge
Authors: Siying Ma
Mehrdad Momeni Zadeh
Mauricio Soroco
Wuyang Chen
Jiguo Cao
Vijay Ganesh
Venue: ICLR 2026 Poster (2026)
Keywords: Neural Operator
PDE
Fundamental Physics Knowledge
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We propose to incorporate fundamental physics knowledge into learning neural operators to enhance its data efficiency, long-term consistency, and OOD generalization.

Abstract: Recent advances in scientific machine learning (SciML) have enabled neural operators (NOs) to serve as powerful surrogates for modeling the dynamic evolution of physical systems governed by partial differential equations (PDEs). While existing approaches focus primarily on learning simulations from the target PDE, they often overlook more fundamental physical principles underlying these equations. Inspired by how numerical solvers are compatible with simulations of different settings of PDEs, we propose a multiphysics training framework that jointly learns from both the original PDEs and their simplified basic forms. Our framework enhances data efficiency, reduces predictive errors, and improves out-of-distribution (OOD) generalization, particularly in scenarios involving shifts of physical parameters and synthetic-to-real transfer. Our method is architecture-agnostic and demonstrates consistent improvements in normalized root mean square error (nRMSE) across a wide range of 1D/2D/3D PDE problems. Through extensive experiments, we show that explicit incorporation of fundamental physics knowledge significantly strengthens the generalization ability of neural operators.
We will release models and codes at https://sites.google.com/view/sciml-fundemental-pde.

Published as a conference paper at ICLR 2026
LEARNING DATA-EFFICIENT AND GENERALIZABLE NEURAL OPERATORS VIA FUNDAMENTAL PHYSICS KNOWLEDGE
Siying Ma', Mehrdad M. Zadeh', Mauricio Soroco' !
Wuyang Chen', Jiguo Cao Vijay 2anesh
Simon Fraser University, Georgia Institute of Technology
ABSTRACT
Recent advances in scientific machine learning (SciML) have enabled neural opera- tors (NOs) serve powerful surrogates for modeling the dynamic evolution of physical systems governed by partial differential equations (PDEs). While existing approaches focus primarily on learning simulations from the target PDE, they often overlook more fundamental physical principles underlying these equations. Inspired by how numerical solvers are compatible with simulations of different settings of PDEs, we propose multiphysics training framework that jointly learns from both the original PDEs and their simplified basic forms. Our framework en- hances data efficiency, reduces predictive errors, and improves -itto---dstribution (OOD) generalization, particularly in scenarios involving shifts of physical pa- rameters and synthetic-to-real transfer. Our method is architecture-agnon demonstrates consistent improvements in normalized root mean square error (nRMSE) across wide range of 1D/2D/3D PDE problems. Through extensive experiments, we show that explicit incorporation of fundamental physics knowl- edge significantly strengthens the generalization ability of neural operators. We will release models and codes at https://sites, wttp://.ooo/eee//simmiiiiieeeeee
and
INTRODUCTION
Diffusion Navier Stakes
Recent advances in scientific machine learning (SciML.) have broadened traditional machine learning (ML) for mod- cling physical systems, using deep neural networks (DNNs) especially neural operators (NOs) (Li et al., 2021a; Pathak al., 2022; Lam et al., 2023; Bi al., 2023) fast, ac- curate surrogates for solving partial differential equations (PDEs) (Raissi al., 2019; Edwards, 2022: Kochkov al., 2021; Pathak et al., 2022). However, compared to numer- ical methods, a key disadvantage of recent data-driven SciML models is their limited integration of fundamental physical knowledge.
Scientific Machine
Siwshiisky Real World Flo
Learning
Fundamental
Diffusion Advection
Knowledge
Can SciML models (c.g., neural operators trained advanced PDEs) also understand fundamental physics knowledge (basic terms like diffusion, advection)?
Figure 1:
Numerical solvers, though tailored specific PDEs discretizations, inherently preserve physical laws (e.g., conservation, symmetry), ensuring consistent and plausible simulations across diverse conditions (physical parameters. boundaries, geometries, etc.) (Ketcheson al., 2012; Hansen et al., 2023; Mouli et al., 2024; Holl & Thuerey, 2024). In contrast, data-driven models, despite learning across PDE types (e.g., via multiphysies pretraining in SciML foundation models (McCabe al., 2023; Hao et al., 2024)), remain sensitive to training distributions, degrading under distribution shifts (Subramanian et al., 2023; Benitez et al., 2024) and requiring large, diverse datasets. This fragility is worsened by the absence of rigorous verification: Unlike classical solvers, SciML models are rarely evaluated against decomposed PDE components. This gap introduces three major challenges: High data demands: Without physics priors, neural operators require large, diverse datasets for high precision, as seen in recent SciML foundation models(Hao et al., 2024; McCabe
"Co-corresponding authors.

---

Published as a conference paper at ICLR 2026
et al., 2023) which focus on generalization without addressing training data efficiency. 2) Physical inconsistency: Lacking inductive biases, these models may violate conservation laws or produce unphysical outputs, particularly in long-term rollout predictions. 3) Poor generalization: Neural PDE solvers often struggle with unseen simulation settings and requires retraining.
Motivated by the above challenges, ask two scientific questions:
Q1: Can neural operators understand both original PDEs and fundamental physics knowledge?
Q2: Can neural operators benefit from explicit learning of fundamental physics knowledge?
In this paper, we highlight the importance of enforcing the learning of fundamental physical knowl- edge in neural operators. The key idea is idenify physically plausible busic terms that can be decomposed from original PDEs, and incorporate their simulations during training. Although often overlooked in SciML, our experiments demonstrate that these fundamental physical terms encode rich physical knowledge. Not only can they be utilized without incurring additional computational costs, but they also widely offer substantial and multifaceted benefits. This opens up new door to improve the comprehensive generalization of neural operators with data efficiency.
We summarize our contributions below:
Through comprehensive evaluations of public SciML models, we observe strong correlation between performance on original PDEs and basic PDE terms, highlighting the importance of fundamental physical knowledge in neural operators (Section 2.2).
2. We propose to explicitly incorporate fundamental physical knowledge into neural operators. Specifically, we design simple and principled multiphysics strategy to train neural operators on simulations of both the original PDE and its basic form. (Section 3).
3. Our method exhibits three major benefits: 1) data efficiency (Section 4.2), 2) long-term physical consistency (Section 4.3), 3) strong generalization in OOD (Section 4.4) and real-world (Section 4.5) scenarios. We evaluate our method on a wide range of ID/DD/3D PDEs, achieving consistent improvement in nRMSE (normalized root mean squared error, Section 4.2).
2 BACKGROUND
2.1 DEFINITION OF PDE LEARNING IN SCIML
>>>RR Rdd, defined on the temporal domain T, the spatial domain S, and the parameter space θ, with d the number of dependent variables. Numerical solvers compute (p -) from E > 1 past steps, enabling finite- difference approximations: No [0g (t (t- l- ",........ vp(t St, + vs(t,-), -), where At is the temporal resolution. SciML aims learn a surrogate operator Ne.s. parameterized by physical parameters 0 and Icarnable weights Ø, that approximates this mapping. Given N simulations, D = (o(0) ([0: fmax),-) -== 1,.... N). N.. the model is trained by minimizing a loss L, often the normalized root mean squared error: nRMSE- Mpa 222 where Upred is the model prediction.
For time-dependent PDEs, the solution a vector-valued mapping v 7 S θ ->
2.2 MOTIVATION: NEURAL OPERATORS EXIIBIT CORRELATED YET WORSE PERFORMANCE ON FUNDAMENTAL PHYSICS
Navier Stokes (nRMMSE)
0.06
Hymna
We begin with a motivating example to highlight the importance of incorpo- rating fundamental physical knowledge into neural operators. Specifically, we gather publicly released pretrained neu- ral operators, with a focus on SciMI. foundation models that are jointly pre- trained across multiple PDEs ("multi- physics") (McCabe et al., 2023; Hao et al., 2024). We evaluate these models on their ability capture fundamental physical knowledge (formally defined in Section 3),,, which none of them
0.05
DPOTT
0.04
DIOTL
DROTEM XX
0.03
HPP
0.02
0.02
MPPLL
0.00
0.10
0.25 Decomposed Convection Term of 2D Navier Stakes (nRMSE)
0:15
Figure 2: On 2D incompressible Navier Stokes, neural opera- tors and SciML foundation models (MPP (MeCabe et al., 2023). DPOT (Hao al., 2024), Hyena (Patil al., 2023)) exhibit cor- related yet worse performance on fundamental physics (x-axis).
2

---

Published as a conference paper at ICLR 2026
explicitly trained on, and compare their performance against the original PDE simulations.
From Figure 2, we observe a strong Pearson correlation (0.9625) between errors on original PDEs and their basic terms. This suggests that stronger SciML models implicitly learn basic PDE components more effectively. However, because these terms are not explicitly included their training data, their absolute errors 0.133-00.00 on the x-axis) remain much larger than those for the original PDEs (0.008-0.) on the y-axis). This gap reveals a key limitation: while the models demonstrate transfer across multiple physics, they still fail reliably capture the fundamental PDE components that underpin complex equations. This motivates us to explicitly enforce an understanding of fundamental physical knowledge within neural operators.
METHODS
Motivated by our observation in Section 2.2, incorporate fundamental physics knowledge into learning neural operators via: 1) Defining and decomposing basic forms from the original PDE (Section 3.1); 2) Jointly training neural operators simulations from both basic forms and the original PDE (Section 3.2). We provide an overview of our method in Figure 3.
Basic Forms
Cheap Simulations
Benefits (Sec. 4.2-4.5)
[Sec, 3.1)
(Table 1)
Data Efficiency
Joint Training
fmixed simulations)
Neural Operator
Long-term Consistency
Sec. 3.2
Generalization (OOD, Synthetic-to-Reall
PDE
Heavy Simulations
Figure 3: Overview of our method. Decomposed PDEs encode rich fundamental physical knowledge and introduce cheaper simulations. By jointly training on both the full PDE and its decomposed basic form, we bring multiple benefits neural operators.
3.1 FUNDAMENTAL PHYSICAL KNOWLEDGE VIA DECOMPOSED BASIC PDE FORMS
3.1.1 DEFINING FUNDAMENTAL PHYSICAL KNOWLEDGE PDEs
Let us consider the second-orde order PDE as general example:
ai u, Vu) Ōr,0z, 7
Su " Σ bila u, Vu) du
clx u, Vu) = f(x),
(1)
where u is the target solution, with æ E R the physical space (e.g., n PDEs). The coefficients aijy b,, ("physical parameters") govern the dynamics; mismatches between training and testing values cause domain shifts, leading to out-of out-of-distributtin (OOD) simulations. Finally, f denotes an external forcing function (Nandakumaran & Datti, 2020).
n= 3 for 2D time-dependent
We establish a systematic process to define the fundamental physical knowledge of PDEs, namely basic PDE terms: 1) retain terms that govern the essential and dominant physical dynamics; 2) remove terms that induce solver stiffness, increase computational cost. or contribute little to the pattern formation of interest. This procedure typically yields a simplified PDE form that can be simulated far more efficiently, while still capturing the key physical dynamics of the original system. From machine learning perspective, decomposing PDE into basic form acts as a data augmentation strategy that reduces data collection costs.
Importantly, our approach of discovering and incorporating basic PDE forms differs fundamentally from the multiphysics training in recent SciML foundation models (McCabe et al., 2023; Hao et al., 2024), which simply aggregate diverse or even weakly related PDE systems. contrast, we stress that neural operators must support basic PDE terms foundation while learning complex PDEs. Figure 4 illustrates the PDEs studied alongside their corresponding decomposed basic forms. An ablation study on different choices of terms is also provided in Appendix D.6.
3

---

Published as a conference paper at ICLR 2026
Diffusion-Reactio Reaction
Navier-Stokes (2D)
Navier-Stokes (3D)
Sitammoto-Siashinnee
0.907
0.404
63.77
PDE
(o1)
[b1]
(e1)
(d1)
0.00
BBsicc Form
%%/a
0,6
1-10, R2
64-128
149.3
106
(b2)
(d2)
0.0
-0.338
-o-teetic-t--Ralll
OOD Settings
(A3)
(b:)
30.3
(d3)
Synthetic
Pojo.
-0.3167
v= 0.05
-0.04474
(c:3)
0 =
kg
1--22
3.912
0.439
0.0
(a4)
(b4)
(d4)
ScalerFlow
Ro/o,"
5,55
0.001
-0.104
k,,
Figure 4: Visualizations of simulations of PDEs and their decomposed basic forms (Section 3.1). From left to right: Diffusion-Reactonn (activator concentration), 2D Navier-Stokes (fluid velocity). 3D Navier-Stokes (smoke density), Juramoto-Sivashing (perturbation amplitude). Basic PDE forms are used training neural operators with fundamental physics knowledge (Section 3.2), and the OOD settings are used for evaluating the generalization of neural operators (Section 4.4). De, Da: diffusion coefficients (Equation 2). y. viscosity (Equation 4), b, ln, by: magnitudes and wavenumbers (Equation 7).
3.1.2 DIFFUSION- DERESSION-REACTION
The Diffusion-RRaccioo equation models an activator-inhiiitor system, which typically happens in the dynamics of chemistry, biology, and ecology. The Inffussion-Racctin equation describes spatiotemporal dynamics where chemical species or biological agents diffuse through a medium and simultaneously undergo local reactions. These systems often used to model pattern formation, such as Turing patterns, in domains ranging from chemistry to developmental biology.
deu Duduu + Dud + Ru: 0, = D. + D,0 v R.
In Equation 2, and represent the concentrations of activator and inhibitor, respectively, and D., D. are diffusion coefficients. The nonlinear reaction terms Ra, R, model biochemical interactions. In experiments,
variant with R,, v) = u k and R, (u, v) = u v, where & = 5 x 10-3 consistent with values used in models of excitable media such neurons or cardiac tissue.
adopt the FitzHugh-Nagun
Decompased Basic Form. To isolate the fundamental transport behavior and reduce simulation cost, we consider simplified form of Equation 2 by omitting the nonlinear reaction terms R, and Rg, yielding pure diffusion equations:
d.u DuOrzu DuO O,v D,Orov + Dr0gy
Why Drop Reaction Terms? This form retains the essential dispersal dynamics but eliminates the feedback coupling between u and Nonlinear reaction terms can vary rapidly, introducing stiffness into the PDE. This stiffness necessitates smaller time steps for stable numerical integration, increasing computational cost. By omitting these nonlinear terms, the system becomes linear and more amenable to efficient numerical solutions.
Why Prioritize the Diffiusion Term? Pure diffusion, though simpler, encodes key properties such as isotropic spreading and mass conservation, providing inductive bias for learning. Unlike reaction terms, which act locally to update activator and inhibitor concentrations, the diffusion term governs spatial coupling, the primary source of pattern formation and spatial dynamics, facilitating transport and stabilization, which explains the visually similar patterns in Figures 4 al and a2.
Physical Scenarios. The emergence of spatial patterns in -ifffioonn systems governed by the ratio
D.. D which affects the relative spreading rates of inhibitor and activator (Page al.,
4

---

Published as a conference paper at ICLR 2026
2005; Asgari et al., 2011; Gambino al., 2024). Classical Turing instability arises when D, > D.. leading to diffusion- pattern formation, and the inhibitor spreads out while the activator stays localized. Following previous works (Menou et al., 2023), set DN 1 x 10 and focus learning simulations when D= 5, and possible OOD scenarios when D=
D. {1, 100).
D
3.1.3 INCOMPRESSIBLE HAVIER-STOKES
The Sewier-Stokes equations govern the dynamics of fluid flow and serve fundamentals for fluid simulations. It considers both the mass conservation and the momentum of fluid parcels.
du -u V)u u 1 Vp W V
Ōt
where u is the velocity field, v is the dynamic viscosity, p is the fluid density, p is the fluid pressure, and f is the external force field.
Decomposed Basic Form. To isolate fundamental nonlinear transport mechanisms and reduce computational complexity, we simplify Equation 4 by omitting the pressure term (incompressibility via projection) 1 Vp and diffusion term v/
du
- (u V.
(5)
Ot
This form captures inertial advection with external forcing and approximates high Reynolds number flows, where viscous effects are negligible. Such simplifications are analytically meaningful and are often used in turbulence modeling (e.g., inviscid Euler equations). Learning this reduced dynamics can help models internalize conveciiondominaae regimes.
Why Drop Pressure and Diffusion Terms? The pressure term, which enforces fluid incompressibility. requires solving large linear systems and is difficult to parallelize and computationally expensive. Omitting significantly accelerates the simulation. Similarly, the diffusion term in Navier-Stoker often uses explicit Euler integration with substeps, adding complexity. Removing simplifies the simulation further. Moreover, for many visual effects like smoke fire, viscosity is minimal, the diffusion term has little visual impact and can often be omitted without noticeable loss in realism.
Why Prioritize the Convection Term? Computationally, the convection term cheap as it describes the local transport of fluid and there is no need to iterate across the spatial domain. Meanwhile, convection the main driver of motion in most fluid flows, as it transports vorticity and mass. Without it, the fluid would just sit still or respond passively to forces. It captures nonlinear soff-interaction, which is critical for dynamic, complex-lookingg behavior.
Physical Scenarios. In Navier-Stokes, the dynamic viscosity in Equation 4 (or the Reynolds number Re puL where is the density of the fluid, u is the flow speed, L is the characteristic linear dimension) controls the fluid dynamics. It measures the balance between inertial forces (pushes the fluid particles in different directions, leading to chaotic flow patterns) and viscous forces (resists motion and smoothes out differences in velocity, promoting an orderly flow) of a fluid. Following previous works (Schlichting & Gersten, 2016; Kochkov et al., 2021; Page al., 2024), will mainly focus on learning simulations when y = 0.01, and possible OOD scenarios when v E {0.05, 0.001 Note that a smaller will lead more turbulent flows.
3D Extension. In real-world scenarios such as atmospheric or smoke dynamics, buoyancy-driven flows provide additional complexity (Eckert et al., 2019). We extend our setting to simulate 3D incompressible Navier-Stokes in a rising plume scenario (see Figure 4 cl and c3). We simulate how a plume of smoke rises and spreads in 3D box-shaped environment. Smoke is introduced from a small circular inflow region located at the bottom center of the domain, at a steady inflow rate of 0.2 units per timestep. The smoke is carried upward due to buoyancy. This setting tests the method's robustness complex spatiotemporal dynamics in three dimensions.
3.1.4 KURAMOTO- STVASSIISSY
The Kuramoto- Sivashinsky equation is a nonlinear PDE that models the interplay of instability, nonlinearity, and dissipation, making it prototype for studying spatiotemporal chaos. It is used to
Kuramoto-Sivashinsky
5

---

Published as a conference paper at ICLR 2026
simulate phenomena such wrinkled flame fronts, thin fluid film instabilities, and chaotic pattern formation. We consider the following one-dimensional Kuramoto Sivashinsky equation:
Kuramoto-Sivashinsly
tu -d- dezu Oxzzr U, on 00, L] x [0,
Ot
The solution u generally represents the perturbation amplitude of chaotic system. The spatial domain [0, L] is equipped with periodic boundary conditions. We impose the initial condition on x [0, [. as
1
K2
uo(x) cos
+ bos
ac(x)),
where uo is the superposition of two cosine modes and small mean-zero Gaussian perturbations. This initialization, common in KS studies (Papageorgiou & Smyrlis, 1991; Gudorf & Cvitanovie, 2019), combines deterministic cosine modes that inject controlled perturbations with small random noise. a, b, k1. 22 are tunable physical parameters controlling the strength and decay of the initial perturbation. Here, e(x) are i..d. standard normal samples, and we set o = 0.05. Spatial derivatives evaluated spectrally via discrete Fourier transforms, and time integration is performed with fourth-order Runge-Kiun Kutta method.
Decomposed Basic Form. isolate the linear stabilizing and destabilizing mechanisms and reduce computational complexity, we simplify the Kuramoto-Sivashinse equation by omitting the nonlinear advection term - udzu.. u. The reduced equation is
du -dzz Oyerzu, U, on [0, L] X [0.T].
(8)
Ōt
which captures the competition between the destabilizing anti-diffusion term (-0rzu) and the stabilizing diffusion term (- U-)
Why Drop the Nonlinear Advection Term? order in u, so the linearized dynamics govern the instabilities and pattern formation. Removing this term also eliminates costly Fourier transforms, yielding significant speedup in simulation.
For small amplitudes or short times, -uô__u is higher
Why Prioritize High-Order Stabilzzing/Destablliiin Terms? The balance between destabilizing - and stabilizing -dnzaz u determines how fast the chaotic model grows or decays. The fourth-order dissipation is particularly important for controlling stiffness and ensuring smoothness of solutions, making it critical for stable numerical integration. Finally, in practice, many inter- face/turbulence models reduce to precisely this "anti-diffusion + diffusion" structure, so analyzing these terms provides broadly transferable insights.
Physical Scenarios. In the initial condition (Equation 7), the amplitudes a, b control the strength of nonlinearity, while wavenumbers k1, k2 control the growth decay of chaos, yielding regimes I) For large a. b and high k,, k2, the dynamics are predominantly linear, and chaos decays over time. 2) For small b and small k ki, the system develops weakly nonlinear chaos, with instabilities growing slowly. 3) For large a, b and small k1, k2, the system exhibits strongly nonlinear chaotic growth. In our experiments, we focus on learning in the linear regime with a, b E [1, 10) and k1, K2 E (64, 128). and evaluate on two out-of -idtr--distribttion (OOD) regimes: weak nonlinearity (a, b E [0.05, 0.2), k1, k2 E [1) and strong nonlinearity (a, [1, 10), by by [1,3
3.2 JOINT LEARNING WITH FUNDAMENTAL PHYSICAL KNOWLEDGE
After defining our fundamental physics knowledge, we now explain how integrate into learning neural operators in principle from two perspectives: data composition and neural architecture.
1) Data Composition. We jointly train neural operators on simulations of both our PDE and the decomposed basic form as multiphysics training with composite dataset. We summarize our simulations in Table I. Since the simulation costs decomposed basic forms are much cheaper than the original PDE. we can "trade-in" simulations of the original PDE for more simulations of basic forms under the same simulation costs. We define the "Sample Mixture Ratio" as the rate derived from the rate of simulation costs of original PDE with its decomposed basic form while making sure

---

Published as a conference paper at ICLR 2026
Table I: Summary of simulations of PDEs and their decomposed hasic forms. GPU: NVIDIA RTX 6000 Ada.
Spatial Resolution
Temporal Steps
Simulation Costs
Sample Mixture Ratio (PDE: Basic Form)
PDF
Target Variables
(sec. per step)
Diffusion-Reaction (Eq. 2)
1.861 x 10-2 6.610 10-3
128 x 128
100
Activator u, Inhibitor
1:3
Basic Form (Eq. 3)
Navier-Stokes (2D) (Fq.. 4)
2.775
256 x 256
1000
Velocity u, Density
1:24
Basic Form (2D) (Eg. 5)
0,113
Navier-Stokes (3D) (Eq. 4)
Velocity u, Density 1.047
50 x 50 89
150
1:3
Basic Form (3D) (Eq. 5)
0.300
Prramoto-Sivashinsee (Eq. 6)
1.176 x 10-1
Basic Form (Fi.. 8)
512
200
Perturbation Amplitude
9.135 x 10-
1:12
that "exchange" one primary data with the corresponding number of basic form data will maintain comparable or reduced simulation budget when training baseline and our proposed model.
We apply a multi-task formulation, where the model learns from both the original PDE and its simplified basic form. The idea is inspired from curriculum learning (Bengio et al., 2009; Pentina al., 2014) and auxiliary task learning (Liu et al., 2019). In our method, the basic forms act as a simpler, physically motivated auxiliary task that can facilitate more efficient representation learning and accelerate the convergence on the primary task.
2) Neural Architecture. We mainly consider the Fourier Neural Operator (FNO) (Li et al., 2021a) as our neural architecture. However, we make our method agnostic to specific architectures of neural operators: We generally share the backbone of the neural operator for learning both the main PDE and its basic terms, while employing separate final prediction layers for the two tasks. We will discuss the details on model structure in Appendix B, and will include more results using transformer (Dosovitskiy et al., 2020; Tong et al.. 2022; McCabe et al., 2023; Chen et al., 2024) in Appendix D.
4
EXPERIMENTS
4.1 SETTINGS
Our baseline is learning the original PDE problem. In general, our method reallocates half of the baseline's ' simulation budget to simulate the basic PDE terms, with the sample mixture ratio defined in Table 1. To fairly compare with the baseline, adopt the same hyperparameters and optimization costs (number of gradient descent steps). Since our goal is to evaluate performance on the original PDE, we use data from the basic PDE term only during training. All testing is conducted exclusively on the original PDE data with 100 samples (Takamoto et al., 2022). We use Adam optimizer, cosine annealing learning rate scheduler, and nRMSE defined in Section 2.1. We summarize our training details in Appendix C.
4.2 DATA EFFICIENCY
2D Difluvion-Reartion (FNO)
2D Incomprressible Navier-Stnkes (PNO)
6.12
6,15
hommwboed wowwwwoddddee RIOOT
We first study our method with dif- ferent numbers of training samples, and demonstrate that neural operators trained with our method can achieve stronger performance with less train- ing data. We consider the following methods:
6.10
a.c
0.00
0.00
6.3
0.02
Simutacion Costs [Serands)
Simu silonn Costs (Secondi)
3D Incompressitle Navier-Stoios
10 Saashis--aasshiskkyy (PNO)
(PNO)
nereaaiedddee RMSS
Dus
Bormal and settoo
Cur
"Bascline": Neural operators that are only trained on simulations of the original PDE.
3.5045
3,3340
Ours: As described in Section 3.2, we can replace simulations of the original PDE with its decomposed basic form, allowing the total sim-
Simuission Costs (Scconas)
30
Sinu ation Costs (Sccands)
42 65
102
200
Figure S: Joint training neural operators on data of the original PDE and the basic form improves performance and data efficiency. Y-axis: normalized RMSE, X-axis: simulation costs (seconds).
7

---

Published as a conference paper at ICLR 2026
ulation cost of the training data to be comparable or even reduced. See Table for the sample mixture rate for fair comparison.
In Figure 5, we study prediction errors of neural operators trained with different numbers of simula- tions (measured in their simulation costs). We can see that, across PDEs and neural architectures, our method (orange square) is to the lower left of the baseline (blue circle), which means that we can achieve improved prediction errors with reduced simulation costs.
4.3 LONG-TERM PHYSICAL CONSISTENCY
20 Reaction-Uiffusion (FNO)
2D Incompressible Navier-Stokes
Next-frame prediction (Takamoto
(FND)
et al., 2022; McCabe et al., 2023;
Nomelled RMSS
0.07
Hao et al., 2024) is widely adopted evaluation, where imput frames are always ground truth. Meanwhile, autoregressive inference, where the model keeps rolling out further temporal steps with its own output as imputs, is meaningful and more challenging stress test. In autoregres- sive inference, a model forecasts fu- tures with its own (noisy) output as inputs, and thus prediction error will accumulate along rollout steps. We can see that our improvements Fig- ure 5 further persist across five toregressive steps in Figure 6, leading improved long-term consistency.
0.06
Ours
Ours
0.06
0.06
0.04
0.00
0.08
0.52
3.0
Rollout
Rollout
3D Incompressiblc NtviersSStokes (FNO)
1D Kuramoto Sivashinsky
(FND)
Baselline
Baselina
Ours
0.025
Ours
0.020
0.125
0.075
0.100
0.013
0.050
10 1.5 20
2.0
35 4.0
45
10
1.5
20
25
3.0 Rallout
4,5 4.
Rollout
Figure 6: Joint training neural operators data of the original PDE and the basic form improves performance with autoregressive inference at different unrolled steps. Models are evaluated using the best-performing checkpoints from training, shown in Figure 5.
4.4 OUT-OF-DISTE DIN-R--DDTTTIBUIOON (OOD) GENERALIZATION
We next show the benefits of method towards the generalization of neural operators in out-of- distribution (OOD) settings, where the physical parameters used during simulation are significantly shifted. We consider physical scenarios described in Section 3.1, and show results in Table 2. We can see that our method not only improves in-distribution errors, but also generalizes better on unseen physical dynamics (simulations by unseen parameters), leading to more robust neural operators.
Table 2: Comparisons OOD generalization for different training methods. Models are evaluated using the best-performing checkpoints from training, as shown in Figure 5, under comparable simulation cost settings.
Source
Target I
Target 2
PDE
Model
Setting
aRMSE
Setting
nRMSE
Setting
#RMSE
Diffusion-Reaccd (2D)
Bascline
D-
0.0289
D.
0.0413
100
0.0770
Ours
0.0231
0.0331
0.0538
Baseline
0.0487
0.0825
0.0369
Navier-Stokke Stoker-ttokss (2D)
== 0.01
# 0.05
0.0001
Ours
0.0175
0.0222
0.0125
Baseline
0.0675
0.0393
0.0836
Navier-Stokes (3D)
0.0001
Ours
0.0481
0.0329
0.0602
Kuramoto-Sivashinkky (ID)
Baseline
3. b,b=1 1
0.0037
a. b= 0.05 2
0.0021
4. bbbb1~> 10
0.0200
Ours
ly. hg = 64 128
0.0034
1y, =
0.0018
kg
1 32
0.0197
4.5 -ONTEETIC-T--RALL GENERALIZATION
Finally, we test neural operators trained on simulations of 3D Savier-Stoker in real-world scenarios. Essentially, transferring models trained on simulations real observations is synthetic-to generalization problem (Chen et al., 2020; 2021), as domain gaps between numerical simulations and real-world measurements always persist. We study the ScalarFlow dataset (Eckert et al., 2019), which is reconstruction of real-world smoke plumes and assembles the first large-scat scale dataset of realistic turbulent flows. We provide visualizations of synthetic and ScalarFlow data Figure 4.
8

---

Published as a conference paper at ICLR 2026
We show the results and visualize the ground truth as well as the model predictions on smoke plumes from ScalarFlow in Figure 7. We can see that our method outperforms the baseline model and presents qualita- tive comparison of scalar flow predictions on real data, illustrating that our joindy trained model exhibits improved synthetic-to -entetti-tto-ral generalization performance. Please read our Appendix D.2 for more results on 3D Navier- Stokes.
30.3
28.5
29.92
Scalartlow
Baseline aRMSE 0.250
Ours aRMSE 0.213
Figure 7: Visualizations of the last time step in the ScalarFlow and its predictions derived by baseline and our model.
5 RELATED WORKS
Machine Learning for Scientific Modeling Learning-basee methods have long been used model physical phenomena (Lagaris et al., 1998; 2000; Chen & Chen, 1995b;a). Physics-informed informed neural networks (PINNs) (Raissi et al., 2019; Zhu et al., 2019; Geneva & Zabaras, 2020: Gao al., 2021; Ren et al., 2022) incorporate PDEs into loss functions to enforce physical laws, but often struggle with generalization and optimization issues (Krishnapriyan al., 2021; Edwards, 2022). Operator learning methods like Fourier Neural Operators (FNO) (Li et al., 2021a; 2020; Kovachki et al., 2023) and DeepONet (Lu et al., 2019) offer greater flexibility by learning mappings between function spaces but require extensive labeled data (Raissi et al., 2019; Brandstetter et al., 2022: Zhang et al., 2023; Chen et al., 2024; Liu et al., 2026). We adopt a unique approach by evaluating and enhancing SciML through the lens of its compatibility with fundamental physical principles.
Data-Driven Neural PDE Solvers Machine learning has increasingly been used to approximate PDE solutions, with neural PDE solvers trained on diverse scenarios to mimic traditional simulators. Early models used CNNs (Guo al.. 2016; Zhu Zabaras, 2018: Bhatnagar al., 2019), while DeepONet (Lu et al., 2019) introduced neural operator (NO) framework separating input and query encodings, inspiring many extensions (Wang et al., 2021; Hadorn, 2022; Wang et al., 2022; Lin et al., 2023; Venturi & Casey, 2023; Xu et al., 2023; McCabe al., 2023; Hao al., 2024). Advances like FNO (Li et al., 10) LNO (Cao et al., 2023), CNO (Raonic et al., 2024), and KNO (Xiong et al., 2024) have expanded the field, with FNO particularly impactful across applications (Li et al., 2021b; Guibas et al., 2021; Yang et al., 2021; Rahman et al., 2022b;a; Pathak et al., 2022: Liu et al., 2022). Compared with previous neural operator works, instead of naively swapping in cheaper simulations of simplified PDEs, the core merit of our work is to emphasize the multifaceted benefits of explicit learning of fundamental physics knowledge during operator learning.
nut-of-Distribution Generalization in SciML Interest in out-of-distribune (OOD) generalization for scientific machine learning (SciML) has grown recently. Subramanian et al. (2023) showed that fine-tuning neural operators (NOs) on OOD PDEs often requires many OOD simulations, which may be impractical. Benitez et al. (2024) proposed a Heimholtz-speccfic FNO with strong OOD performance, supported by Rademacher complexity and a novel risk bound. Other work includes ensemble methods leveraging uncertainty (Hansen et al., 2023; Mouli et al., 2024), loss functions informed by numerical schemes (Kim & Kang, 2024), and meta-learning for varied geometries (Liu et al., 2023). However, varying PDE types and setups across studies hinder unified insights into OOD generalization for NOs. Our work demonstrates that neural operators explicitly trained with fundamental physics knowledge exhibit improved OOD and synthetic-to-real generalization.
6 CONCLUSION
We present a principle and architecture-agnostie approach to enhance neural operators by explicitly incorporating fundamental physical knowledge into their training. By decomposing complex PDEs into simpler, physically meaningful basic forms and using them as auxiliary training signals, our proposed method significantly improves data efficiency, long-term predictive consistency, and out-of- distribution generalization. These improvements are demonstrated across a variety of PDE systems and neural operator architectures. Our finding highlights the untapped potential of fundamental

---

Published as a conference paper at ICLR 2026
physics inductive bias in scientific machine learning, offering robust and cost-effective pathway to more reliable and generalizable surrogate models in real-world physical simulations.
THE USE OF LARGE LANGUAGE MODELS (LLMS)
LLMs did not play a significant role either the research ideation or the writing of this paper. Their use was limited to correcting minor grammatical issues and typographical errors.
ACKNOWLEDGMENT
This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-ACOO-O5CHI 12331 using NERSC award NERSC DDR-ERCAPP033662 and ASCR-ERCAPO03 1463. Cao's research partially supported by the Discovery Grants (RGPIN- 2023-04057) of the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canada Research Chair program.
REFERENCES
Yazdan Asgari, Mehrdad Ghaemi, and Mohammad Ghasem Mahjani. Pattern formation of the fitzhugh-nagumo model: cellular automata approach. IRANIAN JOURNAL OF CHEMISTRY AND CHEMICAL ENGINEERING (ICCCE) 2011.
Yoshua Bengio, Jérôme Louradour, Ronan Collobert, and Jason Weston. Currieulum learning. In Proceedings of the 26th Anmual International Conference on Machine Learning, ICML '09, pp. 41-48, New York, NY, USA, 2009. Association for Computing Machinery, ISBN 9781605585161. doi: 10. (5.114///5333741555555
URL https: /ttp::///oii. org.orrgg./..oo 100...11111111334411 1153380,
Jose Antonio Lara Benitez, Takashi Furuya, Florian Faucher, Anastasis Kratsios, Xavier Tricoche, and Maarten V de Hoop. aitrof-distribbttiona risk bounds for neural operators with applications to the helmholtz equation. Journal of Computational Physics, pp. 113168, 2024.
Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik. Predic- tion of aerodynamic flow fields using convolutional neural networks. Computational Mechanics, 64:525-545, 2019.
Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Accurate medium- range global weather forecasting with 3d neural networks. Nature, 619(7970):533-538 2023.
Johannes Brandstetter, Max Welling, and Daniel E Worrall. Lie point symmetry data augmentation for neural pde solvers. International Conference on Machine Learning, pp. 2022.
PMLR,
Qianying Cao, Somdatta Goswami, and George Em Karniadakis. Lno: Laplace neural operator for solving differential equations. arXiv preprint arXiv:2303. 10528, 2023.
Tianping Chen and Hong Chen. Approximation capability to functions of several variables, nonlinear functionals, and operators by radial basis function neural networks. IEEE Transactions on Neural Networks, 6(4):904 910, 1995a.
Tianping Chen and Hong Chen. Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems. IEEE transactions on neural networks, 6(4):911-917 1995b.
Wuyang Chen, Zhiding Yu, Zhangyang Wang, and Animashree Anandkumar. Automated synthetic- to-real generalization. In International conference on machine learning, pp. 1746-1756. PMLR, 2020.
Wuyang Chen, Zhiding Yu, Shalini De Mello, Sifei Liu, Jose M Alvarez, Zhangyang Wang, and Anima Anandkumar. Contrastive syn-to-real generalization. arXiv preprint arXiv:2104.02290 2021.
10

---

Published as a conference paper at ICLR 2026
Wuyang Chen, Jialin Song, Pu Ren, Shashank Subramanian, Dmitriy Morozov, and Michael w Mahoney. Data-efficient operator learning via unsupervised pretraining and in-context learning. arXiv preprint arXiv:2402.15734, 2024.
Michael Crawshaw, Multi-task learning with deep neural networks: A survey, 2020. URL https: 09796.
Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al. An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010. 1929, 2020.
Marie-Lena Eckert, Kiwon Um, and Nils Thuerey. Scalarflow: a large-scale volumetric data set of real-world scalar transport flows for computer animation and machine learning. ACM Transactions on Graphics (TOG), 38(6):1-1,, 2019.
C. Edwards. Neural networks learn speed up simulations. Communications the ACM, 65(5): 27-29, 2022.
Gambino, MC Lombardo, Rizzo, and Sammartino. Excitable fitzhugh-nagumo model with cross-diffusion: close and far-from- -quuiiuuuii coherent structures. Ricerche di Matematica, 73 (Suppl 1):137-156, 2024.
Han Gao, Luning Sun, and Jian-Xun Wang. Phygeonet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state pdes on irregular domain. Journal of Computational Physics, 428: 10079, 2021.
Nicholas Geneva and Nicholas Zabaras. Modeling the dynamics of pde systems with physics- constrained deep auto-regressive networks. Journal of Computational Physics, 403:109056, 2020.
Matthew Gudorf and Predrag Cvitanovic. Spatiotemporal tiling of the kuramoto-sivaahinsky equation. In APS March Meeting Abstracts, volume 2019, pp. L70-263, 2019.
John Guibas, Morteza Mardani, Zongyi Li, Andrew Tao, Anima Anandkumar, and Bryan Catan- zaro. Adaptive fourier neural operators: Efficient token mixers for transformers. arXiv preprint arXiv:2111.13587 2021.
Xiaoxiao Guo, Wei Li, and Francesco lorio. Convolutional neural networks for steady flow ap- proximation. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pp. 481-490, 2016.
Patrik Simon Hadorn. Shift-deeponet: Extending deep operator networks for discontinuous output functions. ETH Zurich, Seminar for Applied Mathematics, 2022.
Derek Hansen, Danielle C Maddix, Shima Alizadeh, Gaurav Gupta, and Michael W Mahoney, Learning physical models that can respect conservation laws. In International Conference Machine Learning, pp. 12469-12510. PMLR, 2023.
Zhongkai Hao, Chang Su, Songming Liu, Julius Bemer,, Chengyang Ying, Hang Su, Anima Anandku- mar, Jian Song, and Jun Zhu. Dpot: Auto-regressive denoising operator transformer for large-scale pde pre-training arXiv preprint arXiv:2403.03542. 2024.
Philipp Holl and Nils Thuerey. Differentiable simulations for pytorch, tensorflow and jax. In Forty-first International Conference on Machine Learning, 2024.
David I Ketcheson, Kyle Mandli, Aron J Ahmadia, Amal Alghamdi, Manuel Quezada De Luna, Matteo Parsani, Matthew G Knepley, and Matthew Emmett. Pyclaw: Accessible, extensible. scalable tools for wave propagation problems. SIAM Journal on Scientific Computing, 34(4): C210-C231, 2012.
Taeyoung Kim and Myungjoo Kang. Approximating numerical fluxes using fourier neural operators for hyperbolic conservation laws. CoRR, 2024.
11

---

Published as a conference paper at ICLR 2026
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine leaminn-acceleratter accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 11821))=2101884111 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar, Neural operator: Learning maps between function spaces with applications PDEs. Journal of Machine Learning Research, 24(89):1-97, 2023.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney, Char- acterizing possible failure modes in physics-informed neural networks. Advances in Neural Information Processing Systems, 34:26548 265660, 2021.
Isaac E Lagaris, Aristidis Likas, and Dimitrios Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions neural networks, 955)987--1000, 1998.
Isaac Lagaris, Aristidis C Likas, and Dimitris Papageorgiou. Neural-network methods for boundary value problems with irregular boundaries. IEEE Transactions on Neural Networks, (5):1041-1049, 2000.
Remi Lam, Alvaro Janchez-Gonzarnn Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, et al. Learning skillful medium-range global weather forecasting. Science, 382(6677):141--142) 2023.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33:6755-6766 6666, 2020.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021a.
Zongyi Li, Hongkai Zheng. Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. prrsies-inffomed neural operator for learning partial differential equations. arXiv preprint arXiv:2111.03799 2021b.
Guang Lin, Christian Moya, and Zecheng Zhang. B-deeponet: An enhanced bayesian deeponet for solving noisy parametric pdes using accelerated replica exchange sgld. Journal of Computational Physics, 473:11 177:3, 2023.
Burigede Liu, Nikola Kovachki, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, An- drew M Stuart, and Kaushik Bhattacharya. A learning- -baed multiscale method and its application to inelastic impact problems. Journal of the Mechanics and Physics of Solids, 158: 158:104668, 2022.
Wenzhuo Liu, Mouadh Yagoubi, and Marc Schoenauer. Meta-learning for airflow simulations with graph neural networks. arXiv preprint arXiv:2306. 10624, 2023.
Yifan Liu, Bohan Zhuang, Chunhua Shen, Hao Chen, and Wei Yin. Auxiliary learning for deep multi-task learning, 2019. URL https: / http/xxivv....//////a/ org/ httg//b////// aabs/1900 199. 02214.
Yuqiu Liu, Jingxuan Xu, Mauricio Soroco, Yunchao Wei, and Wuyang Chen. Data-efficient inference of neural fluid fields via sciml foundation model. International Conference on 3D Vision, 2026.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint
2019.
Michael McCabe, Bruno Régaldo-Sain Blancard, Liam Holden Parker, Ruben Ohana, Miles Cranmer, Alberto Bietti, Michael Eickenberg, Siavash Golkar, Geraud Krawezik, Francois Lanusse, et al. Multiple physics pretraining for physical surrogate models. arXiv preprint arXiv:2310.02994, 2023.
Lucas Menou, Chengjie Luo, and David Zwicker. Physical interactions promote turing patterns.
arXiv preprint 23Xv::230../252/ 2023.
12

---

Published as a conference paper at ICLR 2026
Grégoire Mialon, Quentin Garrido, Hannah Lawrence, Danyal Rehman, Yann LeCun, and Bobak T Kiani. Self-supervised learning with lie symmetries for partial differential equations. arXiv preprins arXiv:2307.05443. 2023.
Chandra Mouli, Danielle Maddix, Shima Alizadeh, Gaurav Gupta, Andrew Stuart, Michael W Mahoney, and Yuyang Wang. Using uncertainty quantification to characterize and improve out-of- domain learning for pdes. arXiv preprint urXiv:2403. 10642, 2024.
AK Nandakumaran and PS Datti. Partial differential equations: classical theory with modern touch. Cambridge University Press, 2020.
Jacob Page, Peter Norgaard, Michael P Brenner, and Rich R Kerswell. Recurrent flow patterns as basis for two-dimensional turbulence: Predicting statisties from structures. Proceedings of the National Academy of Sciences, 121(23))e23200077222 2024.
Karen M Page, Philip K Maini, and Nicholas AM Monk. Complex pattern formation in reaction- diffusion systems with spatially varying parameters. Physica D: Nonlinear Phenomena, 202(1-2): 95-115, 2005.
Demetrios Papageorgiou and Yiorgos S Smyrlis. The route to chaos for the kuramoto-sivashinsky equation. Theoretical and Computational Fluid Dynamics, 3(1):15-42, 1991.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. Fourcast- net: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv:2202.112/1 2022.
Saurabh Patil, Zijie Li, and Amir Barati Farimani. Hyena neural operator for partial differential equations, 2023. URL https: orgo abs / 220033063600 16544.
Anastasia Pentina, Viktoriaa Sharmanska, and Christoph H. Lampert. Curriculum learning of multiple tasks, 2014. URL https: http://awwww..aaaaaaaaaaaaaaaeae /a1..b//111111.3.
Md Ashiqur Rahman, Manuel A Florez, Anima Anandkumar, Zachary E Ross, and Kamyar Aziz- zadenesheli. Generative adversarial neural operators. arXiv preprint arXiv:2205.030/1 2022a.
Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadenesheli. U-no: U-shaped neural operators. arXiv preprint arXiv:2204.1 //27, 2022b.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci, Rima Alaifari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning of pdes. Advances in Neural Information Processing Systems, 36, 2024.
Pu Ren, Chengping Rao, Yang Liu, Xun Wang, and Hao Sun. Phycmett: Physics-informed convolutionall-ccurrenl network for solving spatiotemporal pdes. Computer Methods in Applied Mechanics and Engineering, 389: 14399, 2022.
Sebastian Ruder. An overview of tulti-task learning in deep neural networks, 2017. URL https: /176 05098.
Hermann Schlichting and Klaus Gersten. Boundary-layer theory. springer, 2016.
Shashank Subramanian, Peter Harrington, Kurt Keutzer, Wahid Bhimji, Dmitriy Morozov, Michael W Mahoney, and Amir Gholami. Towards foundation models for scientific machine learning: Charac- terizing scaling and transfer behavior. arXiv preprint arXiv:2306. 00258, 2023.
Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 35:1596-1611, 2022.
13

---

Published as a conference paper at ICLR 2026
Zhan Tong, Yibing Song, Jue Wang, and Limin Wang. Videomae: Masked autoencoders data- efficient learners for self-supervised video pre-training. Advances in neural information processing systems, 00000008-109,,, 2022.
Simone Venturi and Tiernan Casey. Svd perspectives for augmenting deeponet flexibility and interpretability. Computer Methods in Applied Mechanics and Engineering, 403: 15718, 2023.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informed deeponets. Science advances, :(44))eab88605, 2021.
Sifan Wang. Hanwen Wang. and Paris Perdikaris. Improved architectures and training algorithms for deep operator networks. Journal of Scientific Computing, 92(2):35, 2022.
Wei Xiong, Xiaomeng Huang, Ziyang Zhang, Ruixuan Deng, Pei Sun, and Yang Tian. Koopman neural operator as mesh-free solver of non-lineat partial differential equations. Journal of Computational Physics. pp. 113194, 2024.
Wuzhe Xu, Yulong Lu, and Li Wang. Transfer learning enhanced deeponet for long-time prediction of evolution equations. In Proceedings of the AAAI Conference on Artificial Intelligence, pp. 10629-10636, 2023.
Yan Yang. Angela F Gao, Jorge C Castellanos, Zachary E Ross, Kamyar Azizzadenesheli, and Robert W Clayton. Seismic wave propagation and inversion with neural operators. The Seismic Record, 12)-1344,3,, 2021.
Xuan Zhang, Limei Wang, Jacob Helwig, Youzhi Luo, Cong Fu, Yaochen Xie, Meng Liu, Yuchao Lin, Zhao Xu, Keqiang Yan, et Artificial intelligence for science quantum, atomistic, and continuum systems. arXiv preprint arXiv:2307.08423, 2023.
Yinhao Zhu and Nicholas Zabaras. Bayesian deep convolutional encoder-decoder networks for surrogate modeling and uncertainty quantification. Journal of Computational Physics, 366:415- 447, 2018.
Yinhao Zhu, Nicholas Zabaras, Phaedon-Stede Stelioo Koutsourelakis, and Paris Perdikaris. Physics- constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data. Journal of Computational Physics, 394:56-81, 2019.
14

---

Published as a conference paper at ICLR 2026
A SIMULATION SETTINGS
In this section, we detail the simulation settings for the 2D Diffusion-Reaction, 2D and 3D Incom- pressible Navier-Stokes and ID Kuramoto-Sivashinsky (see Table I for the summary). We will also explain how we prepare simulations for the "Spatiotemporal Downsampling" in Section D.I.
To ensure a fair comparison under equivalent simulation costs with the basic form of each PDE, downsample the original PDE simulations both spatially and temporally. We also introduce the settings for the corresponding reduced spatiotemporal resolution simulations. Note that the simulation cost of these downsampled settings is matched that of the basic form, which implies that their sample mixture ratios in joint training remain equivalent.
To further explore the benefits of our multiphysics joint training approach with this reduced spa- tiotemporal resolution simulation strategies (refer to Ours @ Spatiotemporal in Figure 9), we addi- tionally introduce the simulation settings for the spatiotemporally downsampled basic form of the 2D Irffusion-Reaction equation. As we will discuss in Section D.I, our framework is orthogonal to standard downsampling techniques, and combining the two can lead further reductions in simulation cost. This reduction allows for an increased proportion of basic form samples in the training mixture under a fixed computational budget (see Table 3 for details).
A.1 INFUSION-REACTIN
Our simulation setting for Diffusion- Reaction follows (Takamoto et al., 2022). Our solver is Py- Claw (Ketcheson al., 2012) that uses the finite volume method. We set the initial condition as standard normal random noise u(t = 0, xr y) ~ AO, 1.0). We use the homogeneous Neumann boundary condition. We simulate
spatial domain of Ω [-1, 1|3, with resolution 128 128. We simulate 5 seconds and save into 100 temporal steps.
Reduced Spatiotemporal Resolution. When simulating the original Diffusion- Reaction equation at low spatiotemporal grids (yellow curves in Figure 9), we reduce the spatial resolution from 128 x 128 to 96 x 96, and reduce the number of temporal steps from 100 to 50. We then upsample to 128 x 128 x 100 (steps) via bilinear interpolation to match the resolution of simulations of the original PDE. The total simulation interval is maintained at 5 seconds, preserving the underlying physical dynamics.
Similarly,
further simulate our decomposed basic form of Diffusion Reaction low spa- tiotemporal resolution (green curve in Figure 9). Table 3 shows the simulation cost of decomposed basic forms with reduced spatiotemporal resolution and the sample mixture ratio.
Table 3: Summary of 2D Diffusion-Reaction simulation and its decomposed basic forms with reduced spatiotem- poral resolution. "Sample Mixture Rate": We replace simulations of the original PDE with its decomposed basic form with reduced spatiotemporal resolution and make sure the total simulation cost of the training data can be comparable. GPU: NVIDIA RTX 6000 Ada.
Spatial Resolation
Temporal Steps
Simulation Costs (sec. slep)
Sample Mixture Ratio (PDE: Basic Forn)
PDE
Target Variables
Diffusiinn Reaction (Eq. 2) Basic Form (Eq. 3) with
128 x 128 100
1.864 x 10-2
96 x 96
50
Activator Inhibitor
2,300 x 10-3
1:8
Reduced Spatiotemporal Resolution
A.2 2D INCOMPRESSIBLE NAVIER- DAVIER-STOKES
Our simulation setting for incompressible Navier-Stokes follows (Takamoto et al., 2022). Our solver is PhiFlow (Holl & Thuerey, 2024). We simulate in spatial domain of Ω [0, 1]2 with resolution 256 x 256. We simulate seconds with a dt = steps. Our initial conditions ug and forcing term f are drawn from isotropic Gaussian random fields, where the fow-frequency components of the spectral density are scaled with scale and high- frequency components are suppressed with power-law decay by smoothness. For ug scale is 0.15 and smoothnes. is 3.0. For f. scale is 0.4 and smoothness is 1.0. Boundary conditions Dirichlet.
1 and periodically save into 1000 temporal
15

---

Published as a conference paper at ICLR 2026
Reduced Spatiotemporal Resolution. When simulating the original 2D incompressible Navier- Stokes equation at low spatiotemporal grids (yellow curves in Figure 9), we reduce the spatial resolution from 256 x 256 to 100 100. We then spatially upsample to 256 256 via bilinear interpolation to match the resolution of simulations of the original PDE. To reduce the temporal resolution while maintaining the same total simulation time and number of recorded frames, increase the time-step size and proportionally reduce the number of integration steps and output interval. Specifically, we change the time-step from dt 5 x 10 to dt 5 10-4, and reduce the total number of time steps from Rstpss 100,000 to seps 10,000. To preserve the temporal spacing between output frames, we decrease the frame interval from 100 to 10. This ensures the same total simulation duration of 5 seconds and the same number of output frames (1,000). This modification reduces computational cost by roughly 10 times.
A.3 3D INCOMPRESSIBLE MAVIER-STOKE
Our solver is PhiFlow (Holl& & Thuerey, 2024). We simulate in a spatial domain of 52 = [0, 1J with resolution 50 x 50 x 89. We simulate 150 steps with a dt 2 x 10 We set the initial uo as zero and upward buoyancy foreing term fs 5 x 10-4 Unlike the 2D Navier-Stokes, we introduce randomness of the buoyancy forcing term on horizontal directions by uniformly drawing fz fz from [-1,1] 10-4, We set Dirichlet zero boundary conditions.
A.4 ID KURAMOTO-SIVASHING
Our simulation setting for the one-dimensional Kuramoto-Sivashine equation follows (Papageor- giou & Smyrlis, 1991; Gudorf' & Cvitanovic, 2019). Our solver is a pseudospectral exponential time-differencing fourth-order Runge-Kutta scheme, implemented with 64 contour points. We use periodic boundary conditions on the spatial domain Ω [0, L] with L 64T,, discretized with equispaced 512 grid points. Spatial derivatives are evaluated spectrally via FFT (fast Fourier transform), and nonlinear terms are treated in physical space through pseudospectral squaring and differentiation. We simulate up to the final time T 30, using an internal time step di 0.1, and record stspp 200 snapshots uniformly in time.
B MODEL STRUCTURE
VideoMAE
Input Formi
Input (PDE)
We mainly consider the Fourier Neu- ral Operator and the transformer in this study. The model structures are shown in Figure 8. The basic form and the original PDE share all layers but the last prediction layer, which is agnostic to specific architectures of neural operators. This task-specific output layer is very well- well-known and widely adopted configuration in multi-task learning. has been high- lighted in survey papers that hard pa- rameter sharing (one input, shared hidden layers, multiple outputs) is the standard setup for multi-task learning due its efficiency and representa- tional benefits (Ruder, 2017; 7 Craw- shaw, 2020).
FNO
Input
Input |PDE]
Shored
Conwal
[Basic Form)
Positional Fribeddling
Shared
FFT
Tokera
4x
Fourier Layor
Linear
Liner IFFT
Encoderr Transformer Block
12X
GeLUi Inbuue
GeLLy
Liner
Decader Transformor Block
Output (Besio Formi
IFD hnsees founeer Tivastono
8X
(widtth=612 abeads =
(PDE)
Linear
Output (Besic Formi
Output (POE]
Figure 8: Our method is agnostic to specific architectures of neural operators: we always share the backbone of the model between learning the original PDE and its basic form, and decouple their predictions in the last layer.
C MORE IMPLEMENTATION DETAILS
We summarize our training details in Table 4. We conducted our experiments on NVIDIA RTX 6000 Ada GPUs, each with 48 GB of memory.
16

---

Published as a conference paper at ICLR 2026
Table 4: Training details. "DR": Diffusion-Reaction. "NS": Navier-Stokes: "KS": Kuramoto-Sivastiinkky
ZD DR (FNO)
2D TR (Transormer)
2D NS (FNO)
NS
3D NS (PNO)
3D VS (Tranformer)
ID KS (FNO)
(Trnsformer)
Input Stape Format
H W
Txx C (C--2)
H W x
CCC-
7. T x IC-
Tx
Number
Inpui Steps Section Simmlatious)
2.
32,
16,
2,4,
1024
10
10
Sample
Learning Rate
0.001
0.0003
0.001
0.001
0.001
0,00015
0.001
Baach Size for Primasy Diata
16
16
8
64
Egoohs
100
60
30
20
M1
50
Auxiliary Task Liss Weight
0.7
0.7
0.7
01.7
0.7
0.7
0.7
Training llours
0,08-1883
0.8tm-7hr
1-229
1.5:45
0.5w6.5
16-120
2 18
Gradient Descent Steps Per
23
124
70
D MORE RESULTS
Beyond the main comparison between the baseline and our proposed model using FNO showed in Section 4, we further conducted additional experiments to assess our approach. In this section, will show the results from both FNO and Transformer.
D.I COMPARISON OF DATA EFFICIENCY AND OUT-OI-DISTRIBUTION GENERALIZATION AGAINST SPATIOTEMPORAL DOWNSAMPLING
To demonstrate the effectiveness of our method, conducted an ablation study comparing it with Spatiotemporal Downsampling on both 2D Diffusion-Reaction Reaction (DR) and 2D Savier-Stoker (NS). We defined the Spatiotemporal Downsampling method as follows: Neural operators that trained with a mixture of simulations of the original PDE and simulations at low spatial and temporal resolutions (then linearly interpolated the original resolution). Similar also save simulation costs with reduced spatiotemporal resolutions. See Appendix A for rates of downsampling and spatiotemporal downsampling during simulations, our method can serve as a complementary data augmentation. On 2D -iffusion-Reaction, we can simulate our decomposed basic forms at lower spatiotemporal resolution, leading to further reduced simulation costs and improved data efficiency (green triangle), outperforming the baseline at the lower spatiotemporal resolution (yellow diamond). For 2D Navier-Stoker with the Transformer, our method still outperforms both the baseline and spatiotemporal downsampling at comparable cost.
method, we
details. Meanwhile, as our decomposed basic form orthogonal to
2D Inffusion-Reaction (FNO)
2D Incompressible Navier-Stokes
(FNO)
0.12
Normalizer lized RMSE
Baseline@Spato ocemporrl
Nozdalized RMSE
otemporal otemppp Ours
Baseline@Spat
0.10
Ours
0.08
Ours@ Burs@SSpaiioemmpaall
0.10
0.05
0.04
0.02
10t
10*
10a
Simulation Costs (Seconds)
Simulation Costs (Seconds)
2D Riffution-Reaction (Transformer)
Incompressible Navier-Stokes
(Transformer)
0.10
Normalized RMSE
0.8
Normalized RMSE
Baseline
Basel Spatiotemporal
Baseline@ Spaliottpppaaiiitrrrr
Ours
0.08
Ours
Ours@Spatictemporal
0.4
0.06
0.2
0.04
10'
10a
10c
10a
Simulation Costs (Seconds)
Simulation Costs (Seconds)
Figure 9; Joint training neural operators on of original PDE and the basic form improves performance and data efficiency. "Spatiotemporal": short for "Spatiotemporal Downsampling". Y-axis: normalized RMSE. X-axis: simulation costs (seconds).
17

---

Published as a conference paper at ICLR 2026
Table 5 reports the nut-of-distribution (OOD) generalization results across both the 2D Diffusion- Reaction and Navier-Stoker equations. Similar to the results in Table 2, here we can see that our approach not only improves in-distribution errors but also consistently enhances generalization to simulations of unseen physical parameters. This robustness holds across both FNO and Transformer architectures, leading to reliable and consistent neural operators under varying conditions.
Table 5: Comparisons of OOD generalization for different training methods with the Transformer. Models are evaluated using the best checkpoints from training in Figure 5, under comparable simulation cost settings. "Spatiotemporal" short for "Spatiotemporal Downsampling".
Souree
Target
Target 2
PDE
Model
Setting
nRMSE
Setting
nRMSE
Setting
nRMSE
Baseline
0.0289
0.0413
0.0770
Diffusion-Rendd (2D, FNO)
Baseline @ Spatiotemporal
Du
0.0234
0.0303
De
100
0.0663
Ours
0.0231
0.0331
D.
0.0538
Ourv Spatioppaiotempooral
0.0218
0.0298
0.0596
Baseline
0.1056
0.1249
0.1976
Diffusion-Reacte
Baseline Spatiotemporal
0.0542
0.0698
0.0812
100
(2D, Transformer)
Ours
0.0602
0.0782
17,
0.0853
Ours@ Spaltolettottrlralll
0.0469
0.0489
0.0671
Baseline
0.0487
0.0825
0.0369
Navier-Stokes (2D, FNO)
Buseline @ Spatiolemporal
0.01
0.0442
0.05
0.0743
0.0001
0.0269
Ours
0.0175
0,0222
0,0125
Baseline
0.0479
0.0853
0.0685
Navier-Stokes
Buscline Spatiotemporal
0.0496
0.05
0.0568
0.0001
0.0402
(2D, Transformer)
Ours
0.0265
0.0397
0.0256
D.2 DATA EFFICIENCY AND OUT-OF-DISTRIBIUS GENERALIZATION OF TRANSFORMER FOR 3D NAVIER-STOKES
Similar to what have studied in Section 4, aim also demonstrate three key benefits: data efficiency, long-term physical consistency, and strong generalization in OOD simulations using Transformer, for 3D Navier-Stokes well.
In Figure 10, we can see that joint training (orange square) both the original and basic forms of the 3D Navier-Stokes equation consis- tently reduces normalized RMSE from baseline (blue circle) across varying simulation budgets. This improvement is observed for Transformer architectures, highlighting enhanced data effi- ciency and generalization, which aligns with the results in Section 4.2.
3D Incompressible Navier-Stokes
(Transformer)
Normalized RMSE
0.4
Baseline
Ours
0.2
0.1
0.0
In Table 6, we show that our joint train- ing approach significantly improves out-of- distribution generalization on 3D Navier-Stokes across all test settings, outperforming the base- line for both FNO and Transformer models. To- gether with the results in Table 2 and 5, the consistent gains observed across all OOD setting results underscore the effectiveness and robustness of our method in generalizing to previously unseen physical regimes, particularly under significant shifts in simulation parameters.
10a
104
Simulation Costs (Seconds) Figure 10: Joint training neural operators on data of the original 3D Navier-Stoker equation and the basic form improves performance and data efficiency.
D.3 MORE LONG-TERM CONSISTENCY RESULTS OF TRANSFORMER
In our Figure 11, we show the rollout performance of the transformer the 2D Diffusion- Reaction -erfutioonRReattion and 2D incompressible checkpoints from training in Figure 9. Losses will be aggregated for five consecutive time steps.
equations. Here, we run the experiments with the best
18

---

Published as a conference paper at ICLR 2026
Table 6: Comparisons of OOD generalization on 3D NS different training methods using Transformer.
Models are evaluated using the best checkpoints from training in Figure 10.
Source
Target 1
Target 2
PDE
Model
Setting
mRMSE
Setting
nRMSE
Setting
nRMSH
Bascline
1.0114
0.0327
0.0816
3D Ntoker-Stokes
0.01
0.1
0.0001
Ours
0.0064
0.0124
0.0322
We see that our improvements in Figure 9 further persist across autoregressive steps, leading to improved long-term consistency, aligning with our results in Figure 6.
2D Diffusion-Reaction
2D Incompressible Navier-Stokes
(Transformer)
(Transformer)
Normalized aaaiiddd RMSE
0.13
Normalizee RMSE
0.14
Baseline
0.12
0.11
0.12
Ours
0.10
0.10
0.09
0.08
0.08
Baseline
0.06
0.07
Ours
0.04
0.06
1.0
1.5
2.0
2.5
3.5 4.0 4.5 5.0 Rollout
3.0
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
Rollout
Figure 11: Joint training neural operators data of the original PDE and the basic form improves performance with autoregressive inference at different unrolled steps using Transformer. Models are evaluated using best-perforiing checkpoints from training shown in Figure 9.
D.4 MORE RANDOM SEEDS
To ensure the statistical robustness of findings, now FNO using three different random seeds during initialization and training. For each configuration, we report the average performance across the three runs, and include standard deviation error bars in all plots in Figure 12. This enables a more rigorous evaluation of model performance, capturing the inherent variance and mitigating the risk of overinterpretation from single-seed outcomes. We can see that the results demonstrate that joint training of neural operators on data from both the original PDE and its decomposed basic form yields consistent improvements in predictive performance and data efficiency, highlighting the effectiveness of this multiphysics learning strategy.
2D biffusion-Reaction (FNO)
2D Incompressible Navier-Stokes (FNO)
Noraalzedd RMSE
0.12
Basei inc
Normalzzed RMSE
Basel inc
0.10
Ours
Basel e@ Spatiotemporal
0.15
Ours
Basel Ino@ spatiotemporal
0.08
Surs@SSpatictemporn
0.10
0.06
0.04
0.05
0.02
10a
10a
10a
10
Simulation Costs (Seconds)
Simulation Costs (Seconds)
Figure Model performance averaged over three random seeds. Joint training neural operators on data of the original PDE and the basic form consistently improves performance and data efficiency. Dash lines indicate the mean performance, with error hars representing standard deviation. I.egends align with the descriptions in Section 4.2. Columns: (left) 2D iiffusion-Reactionn (right) 2D Navier-Stokes. Y-axis: nRMSE. X-axis: Simulation Costs (seconds).
To further demonstrate the reliability and generalizability of our approach, take the best perforerno checkpoints from Figure 12 and evaluate their rollout accuracy as well as their out-off-istribuute behavior. From Figure 13 and Table 7, observe that the standard deviations remain exceptionally
19

---

Published as a conference paper at ICLR 2026
small, underscoring the stability of the model across different random seeds. The results show that joint training of neural operators on both the original PDE and its decomposed basic form yields consistent improvements in long-term physical consistency and generalizes better to unseen physical dynamics. Together, these findings statistically prove that the proposed multiphysics learning strategy leads to more stable, more reliable, and more physically consistent neural operators.
InffusionnReactiin (FNO)
2D Incompressible Navier-Stokes (FNO)
Baseline
Baseline Ours
0.07
Ours
0.10
Normalized maaalieed MSE
0.06
0.08
Normalized
0.05
0.06
0.04
0.04
0.03
0.02
1
2
3
4
Rollout
Rollout
Figure 13: Model performance averaged over three random seeds. Joint training neural operators on data of the original PDE and hasic form consistently improves long-term physical consistency. Dash lines indicate the mean performance, with error bars representing standard deviation. Models are evaluated using the best-perforinng checkpoints from training shown in Figure 12.
Table 7: Comparisons of OOD generalization for different training methods averaged over three random seeds. Models evaluated using the best checkpoints from training Figure 12, under comparable simulation cost settings, "Spatiotemporal": short for "Spatiotemporal Downsampling"
PDE
Model
Source
Target
Turget
Seting
MRMSE Mean (SD)
Setting
oRMSE Mean (SD)
Setting
#RMSE Mean (SD)
Baseline
0,0287 (0.0066)
0.0411 (0,008))
0.0754 (0.0014)
Dilfusion-Resn Resction (2D, FND)
Baseline @ Spatiotempora
0.0235 (0 (0.0001)
0.0305 (0.0005)
0.0664 (0.0041)
Ours
0.0232 (0.0003)
D.
0.0329 (0.00033
100
0,0532 (0.0014)
Ourste Spulitemporal
0,0216 (0,0005)
0,0296 (0,0002)
0,0548 (0.0053)
0.0522 (0.0030)
0.0855 (0.0025)
0.0367 (0.0002)
Navier-Sookes (2D) PND)
Baseling Spatiotemporal # 0.01
0.0462 (0.0017)
#= 0.05
0.0729 (0.0013)
== 0.0001
0.0272 (0.0001)
Ours
0.0174 (0.0002)
0.0222 (0.0001)
0.0131 (0.0005)
D.5 LOSs REWEIGHTING
In Section 3.2, we define our total loss for joint learning of the original PDE (Losspun) and its fundamental physical knowledge (decomposed basic form LOsSpasic) as:
Loss
LosSpall 0.7 LoSSBasic-
To address the concern regarding the fixed auxiliary loss weight, we conducted an ablation study exam- ining the effect of auxiliary loss weighting in joint training using FNO for the 2D Diffusion- Reaction system. We evaluated three auxiliary weight settings (0.5, 0.7, and 1.0) across the full range of simulation budgets used in the data-efficiency experiments (Figure 5).
As shown in Table 8, the normalized RMSE values are highly consistent across all three weight choices. Importantly, the maximum absolute deviation across weights remains below 0.0063 for all simulation costs, with most gaps around 0.002-0.004. Although performance through exhaustive hyperparameter fine-tuning, the observed variations across auxiliary weights are small relative the overall error scale and substantially smaller than the performance gains achieved by increasing the training budget, which is trivially beyond the scope of our work.
be further improved
This result demonstrates that the model is largely insensitive the specific choice of auxiliary weight and suggests that the improvements achieved through joint training are robust with respect to this hyperparameter. Based this analysis, fix the auxiliary weight to 0.7 for all the experiments.
20

---

Published as a conference paper at ICLR 2026
Table Ablation study on the auxiliary loss weight in joint training using FNO for the Diffusion- Reaction across three settings: auxiliary weights of 0.5, 0.7, and 1.0. Results indicate that our model is robust to the choice of auxiliary loss weighting. Simulation costs are aligned with the ones used in Figure 5 (top-left).
Simulation Costs (Seconds)
Weight = 0.5
Weight = 0.7
Weight =
3.73
0.1255
0.1233
0.1214
7,46
0.0674
0.0659
0.0636
14.92
0.0508
0.0478
0.0445
29.84
0.0349
0.0330
0.0314
59.68
0.0274
0.0266
0.0249
119.36
0.0245
0.0238
0.0224
238.72
0.0248
0.0236
0.0221
D.6 CHOICE OF THE FUNDAMENTAL TERM
To demonstrate the importance of the choice of dropping terms, we conduct an ablation study on 2D Diffusion-Reactionn and 2D Navier-Stokes. For Diffusion-Reaction, simulate reaction term instead of the fundamental diffusion term. The simulation cost for reaction-only term is 2.048 x 10-3 seconds per step, corresponding a 1:9 sample mixture ratio when compared to the simulation cost of the original PDE data. For Navier-Stokes, we simulate diffusion term instead of the fundamental advection term. The simulation cost for diffusion-onie onnlloyononyyooooo term is 0.234 seconds per step, corresponding to 1:12 sample mixture ratio, making sure the results comparable.
From Figure 14, we can find that 2D the fundamental diffusion term will damage the accuracy with up to 64% increase of nRMSE compared to the baseline, while applying the fundamental diffusion term keeps boosting the model performance with 11% to 24% decrease of nRMSE. Similarly, keeping diffusion term but dropping advection term in 2D Navier-Stoker Stokes can damage the accuracy up to 13% increase of nRMSE compared to baseline and up to 179% increase compared to the one using the functional advection term. This ablation study confirms that the correct fundamental basic term can improve the data-efficiency when joint training with the original data, and proved that the source of improvement in Figure is clearly from training with the fundamental term itself.
keeping the reaction term and removing
2D Diffusion-Reaction (FNO)
2D Incompressible Navier-Stokes (FNO)
0.200
Baseline
0.16
Baseline
Ours@
Ouvs@Advection
0.150
Ours Reaction
Narnalized linad RWSE
0.12
Diffuuiffusion
Normalized
0.125
0.10-
0.08
0.755
0.06-
0.050
0.025
10'
Simulation Costs (Seconds)
Simulation Costs (Secands)
Figure 14: Ablation study of joint training neural operators with two different decomposed terms: the fundamental diffusion term, and reaction term in 2D Diffusion-Reaction. and the fundamental advection term and diffusion term 2D Navier-Stokes shows the importance choice on fundamental terms from PDE equation. Y-axis: nRMSE. X- X-axis: Simulation Costs (Seconds).
To further clarify that the choice of fundamental term is not ad hoc, we also provide principled and quantitative criterion for determining which decomposed components should be retained. For each candidate fundamental operator, we compute the RMSE between the full PDE simulation and the corresponding decomposed simulation over the entire simulation trajectory, using identical initial conditions to ensure fair comparison. This RMSE here serves as a distance metric that quantifies how much each component contributes to the overall dynamics. In 2D DiffusionnRRectionn. the RMSE between the full equation and the diffusion-only simulation 0.553, whereas the RMSE with the reaction-only simulation is much larger at 5.821, which further confirms the results of the qualitative visualization of the original PDE and two decomposed term we obtained in Figure 15.
21

---

Published as a conference paper at ICLR 2026
Initial
t= 2.5
t=5
4.078
0.556
0.607
PDE
-06022
-0.636
Basic Form
0.2420
0.1687
-4.15
0.2716
-0.1992
Residual (Non- Mon-Basic Form)
4.078
0.908
0.576
-.155
-09099
-0.584
Figure 15: Qualitative visualization of the original PDE and its decomposed forms for 2D Diffusion-Reaction This results combined with the quantitative results obtained from the RMSE experiment demonstrate our choice of fundamental term is principled.
D.7 DPOT (HAO ET AL., 2024) ON 2D DIFFUSION-REACTIOO REACTION
To further demonstrate the scalability of our method to larger models and datasets, we addi- tionally evaluate our method using DPOT (Hao et al., 2024), which is a state-of-the-an founda- tion model designed specifically for large-scale PDE pretraining. We use the publicly released medium-sized DPOT checkpoint. adapt our method structure, applied the model with the -ask-specific output layers, which contains 124M model parameters. This represents sub- stantial increase scale compared to the 0.46M- parameter FNO used in our primary experiments described in Appendix B. We also increase our fine-tuning sample size from 2-128 (See Ta- ble 4) to 32-256 for DPOT.
2D Reffutoon-Reaction (DPOT)
Baseline
Normalized RMSE
0.008
Ours
0.006
0.004
0.002
100
200
300
400
Simulation Costs (Seconds)
Figure 16: Data efficiency on 2D Diffusion-Reaction using the DPOT foundation model, suggesting joint train- ing the original and basic form yields consistently lower error across simulation-cost budgets
We then fine-tune DPOT under the same simulation-cost budget for both baseline and our method. Same as the setting in Section 4, we apply our method by reallocating half of the simulation budget to generate auxiliary data from the decomposed basic PDE forms based on the sample mixture ratio defined in Table 1. All experiments use the default DPOT fine-tuning configuration, and apply the same optimization cost (number gradient descent steps) for both baseline and our method, ensuring fair and consistent comparison. As shown in Figure 16, our approach continues to improve data efficiency and achieves consistently lower nRMSE compared to fine-tuning DPOT on original regimes. Notably, these gains persist despite DPOT's
PDE data alone across all simulation-ce significantly larger capacity and extensive pretraining, demonstrating that our multiphysics auxiliary tasks provide benefits beyond what is already captured by large-scale foundation pretraining.
22

---

Published as a conference paper at ICLR 2026
D.8 LIE TRANSFORM ARGUMENT ON 2D IMCOMPRESSIBLE NAVIER- -TOOES
Lie symmetries offer a way to generate new, physically valid training examples by exploiting the analytic group transformations that map one PDE solution to another. This enables the model to learn representations that are inherently equivariant to fundamental symmetries such as translation, rotation, and scaling. To further prove the strength of our model, we leverage the implementation of Lie point symmetry augmentation from (Brandstetter et al., 2022: Mialon et al., 2023), which is orthogonal to our multiphysics joint training approach, to 2D incompressible Navier Strokes equation.
We incorporate the augmentation process into our model. We only apply Lie- transform augmentations exclusively to the velocity (u) of the original 2D incompress- ible Savier-Stokes, leaving the remaining density and all target variables from the de-
2D Incompressible Navier-Stokes (FNO, Lie Transform Augmentation)
Baseline
Normalized RMSE
0.15
Baseline@ueTranstormma
Ours
Ours@uueTranstommuu
0.10
composed basic forms unchanged. Follow-
ing (Mialon et al., 2023), the Lie transfor-
0.05
mation is implemented with second-orde Lie-Trotter splitting scheme with two steps,
10a
10P
where the five fields (x, y, L, Uz Uy) were transformed in accordance with the sampled generator strengths as follows: a maximum time shift (g1) strength of 0.1, maximum spa- tial translations (92. 93) strength of 0.1 and y respectively, a maximum scaling (g4)
Simulation Costs (Seconds)
Figure 17: Joint training neural operators data of the original PDE and the basic form, as a complementary data augmentation orthogonal to Lic-transform augumentation, can further improve performance and data efficiency. Y-axis: normalized RMSE. X-axis: simulation costs (seconds).
strength of 0.05, a maximum rotation (gs)
strength of 10", corresponding x/ 1/18 radians, a maximum x-linear boost (ge) and y-linear boost (9r) strength of 0.2 and a maximum x- and y-quadratic boosts (gs, 99) strength of 0.05.
As our decomposed basic form is orthogonal to the Lie point symmetry augmentation, method can serve as a complementary data augmentation. In Figure 17, we study prediction errors (nRMSE) of neural operators trained with different numbers of training samples (simulations). As have already seen (Figure 5), our approach (orange square) significantly outperforms the baseline (blue circle). In contrast, the Lie-transform augmentation alone (yellow diamond) only marginally improves the baseline. As a result, combining our approach with Lie transformations (green triangle) yields strong performance, but is comparable with our approach alone, underscoring the orthogonal and complementary benefits of these two techniques.
D.9 VISUALIZATION OF PREDICTIONS
To show the predicted PDE solution from our jointly training neural operators on original PDE equation and its basic form aligns with the ground truth, present qualitative visualizations of model predictions across three PDEs, 2D Diffusion-Reaction. 2D and 3D Incompressible Navier- Stokes, and 1D Siraaoot-SSvashhnskyy in Figure 18. For the first three cases, the initial state and predicted states at intermediate and final rollout times are shown. For 1D Kuramoto-Sivashinse we show the whole predicted trajectory here. The predictions are generated using the FNO model trained with our joint training framework. Across all systems and time points, the predictions closely align with the expected dynamics, accurately capturing both spatial patterns and temporal evolution. These visualizations highlight the model's 's capacity to generalize across scales and exhibit physically coherent behavior.
23

---

Published as a conference paper at ICLR 2026
Initial
2:5
t-
6,078
0.506
0,007
Diffusion-Reand -Raactnn
-4.136
-0.002
-0.336
0.509
0.612
Prediction (FNO)
-0.002
-0536
Initial 2.982
2.5
1=5
0.3330
02708
Navier- Navier-Stokes Stokes (2D)
-2.313
0.3316
0.2681
Prediction (FNO)
68
63.78
Navier Stokes (3D)
0.00
0.00
88.67
63.38
Prediction (FNO)
0.00
10
Kuramoto- Sivashinsky (1D)
235-
Ground Truth
10
Prediction (FNO)
0
250
Figure 18: Qualitative visualization of model predictions for 2D CnrfusionnRcacttionn 2D Incompressible Navier- Stokes, 3D Incompressible Navier-Stoke systems and ID Kuramoto Sivashinsky using FNO trained with joint framework. For the first three cases, the initial state and predicted states intermediate and final rollout times are shown. For Suranoto-Sivastinkkk. we show whole predicted trajectory here. results demonstrate accurate temporal evolution and spatial coherence.
24


Reviews:
[{"review_id": "YQDNMqJia5", "reviewer": "Reviewer_FhDJ", "summary": "The paper discusses basic forms of standard PDEs. It finds out that if put half of the data generation efforts into generating basic form of PDE, the performance and generalization will generally improves. The paper test on diffusion-reaction, 2d NS, 3d NS, and KS equation. The improvement on 2D NS is the most significant.", "questions": "- In the experiment, we need to be a bit careful how we measure the simulation cost (Table 1). In practice, the runtime of numerical solvers depends on the choices of gridsize, timestep, and numerical tolerance etc. Here how the simulation cost is measured? \n- If we instead add low fidelity data (with lower gridsize) with the same equation, would the performance improve?\n- How about adding smaller RE on NS, it would requires low simulation cost too?\nIn general it would be helpful to add a bit more discussion on the simulation cost.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- The paper studies an important but often overlooked probably on efficient data generation for neural operators.\n- When comparing against the baseline, the paper controls the overall run time budget.\n- The experiment show consistent improvement.", "weaknesses": "- It is not always obviously which is the best basic form to each target equation. I don't think there is a canonical choice, and the performance depends on the choise of augmented basic PDE. For example on 2D NS the treatment is more significant, but not as much on KS. It would be better to add some ablations to justify the choice.\n- Overall I like this paper but I don't like the storytelling. The main message should be \"it is helpful to generate additional data in simpler form\". It is a bit speculative to claim about \"fundamental physics knowledge\". In my opinion it can be awkward to say diffusion is the \"fundamental physics\" to diffusion reaction equation, or convection is \"fundamental physics\" to Navier Stokes equation.\n- Line 108, the authors use Spearman correlation which is defined for ordinal (rank) correlation. Instead it would be better to use Pearson correlation as the rank does not matter here.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "g7p98K39ls", "reviewer": "Reviewer_kRiv", "summary": "The paper proposes a multiphysics training scheme that jointly learns from full PDE simulations and their decomposed “basic forms” the authors argue that this injects “fundamental physics knowledge” into neural operators (NOs) and improves data efficiency and OOD generalization.  The key contribution lies in identifying and leveraging \"fundamental physics knowledge\" through decomposed basic PDE forms.  This has not  been explored extensively in the neural operator literature.  \nThe authors target two central SciML issues, 1. data hunger and 2. poor OOD transfer, for operator learning across 1D/2D/3D PDEs (Diffusion-Reaction, Navier–Stokes, Kuramoto–Sivashinsky, plus ScalarFlow). \nFormulations of PDEs and “basic forms” are standard and correctly specified\nThe paper is generally well written with helpful overview figures (Fig. 3 pipeline; Fig. 4 gallery of PDEs/basic forms) and plots tying simulation cost to nRMSE. Implementation, data splits, and training schedules are placed in appendices. \nMinor typos remain but do not impede readability.\nCentral claims are supported with experiments and results that are a bit light on content\nThe validation of physics (central theme) is light given no explicit checks on mass/energy conservations.  Another issue is the heuristic treatment of the fundamental physics term.", "questions": "See weakness section , addressing those would be good\nrecommend to:\nBetter justify the decomposition principles\nProvide theoretical analysis or at least intuition for why the approach works\nCompare with more baselines\nDiscuss limitations and failure cases more thoroughly", "limitations": "", "rating": 4, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The authors target two central SciML issues, 1. data hunger and 2. poor OOD transfer, for operator learning across 1D/2D/3D PDEs (Diffusion-Reaction, Navier–Stokes, Kuramoto–Sivashinsky, plus ScalarFlow). \nThe key contribution lies in identifying and leveraging \"fundamental physics knowledge\" through decomposed basic PDE forms.  This has not  been explored extensively in the neural operator literature.  \nProposed benefits such as:\nData efficiency, Long horizon stability , OOD generalization , are all desirable.", "weaknesses": "The term \"fundamental physics knowledge\" is somewhat vague and could be better defined\nSection 3.1 could be more systematic in explaining the decomposition principles\nSome notation inconsistencies (e.g., switching between v and u for solutions\n\nMissing error bars in main results (added later in appendix)\nthere is limited statistical analysis\nThe ScalarFlow experiment (Section 4.5) is somewhat disconnected and brief\nClaims about \"fundamental physics knowledge\" being key are not fully validated (could be just multi-task regularization)\nThere is a lack of theoretical insight, no formal explanation of why the approach works beyond intuition is provided.\nThe decomposition rules in Section 3.1.1 seem ad-hoc without principled justification\n\nEvaluation is very basic, results only compare against vanilla baseline and spatiotemporal downsampling\nThere is no comparison/discussion with other data-efficient methods or recent foundation models\nReal-world evaluation is very limited (only ScalarFlow)\nInconsistent terminology: \"Fundamental physics knowledge\" vs \"basic forms\" used interchangeably\nMissing details: How are the mixture ratios exactly determined? Training time comparisons?\nLimited discussion: When would this approach fail? What about PDEs that don't decompose nicely?\nPresentation issues: Some figures (especially in appendix) are too small to read clearly\nScalability concerns: All experiments on relatively small-scale problems\n\nGrammer-Typos\nline 483: and outha ha h-of-distribution generalization.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Rcm76I4eQ9", "reviewer": "Reviewer_dWHm", "summary": "This paper proposes a data-efficient learning framework for PDE dynamics forecasting by jointly learning from both the original PDEs and their simplified basic forms. Extensive experiments on a wide range of 1D/2D/3D PDE problems demonstrates the effectiveness of the proposed framework.", "questions": "See weaknesses.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The proposed method is well-motivated. The authors provide a critical observation by evaluating existing SciML foundation models. They find a strong correlation between a model's performance on the original PDE and its performance on the fundamental components of that PDE (e.g., pure diffusion for a reaction-diffusion system). However, the absolute error on these basic terms remains high, indicating that even powerful models lack a robust understanding of the foundational physics, which motivates the need for explicit training on these concepts.\n- Methodological Innovation:​​ The paper proposes a simple yet effective multiphysics training framework. It first derive a \"basic form\" from the original PDE by retaining terms governing essential dynamics and removing terms that cause computational stiffness or high cost. The model is trained on a composite dataset from simulations of both the original PDE and the basic form.", "weaknesses": "- Heuristic Nature of Decomposition:​​ The process for selecting terms for the \"basic form,\" while physically intuitive, remains heuristic. A more formalized principle or an ablation study discussing the impact of alternative decompositions more prominently would strengthen the methodology. \n- Inadequate Mechanistic Explanation for the Efficacy of Basic Form Data​: A significant weakness of the paper lies in its insufficient exploration of the underlying mechanisms by which the \"basic form\" data aids the learning of the original PDE. The attribution of performance gains solely to the incorporation of \"fundamental physics knowledge\" is a high-level concept that lacks granularity. A more rigorous analysis is required to dissect how the basic form data contributes. A possible explanation is that data from the basic form may provide more diverse initial conditions. Can the data from the basic form be replaced with an equivalent amount of original PDE data? Although this would incur greater simulation costs, it would help clarify the specific ways in which data from the basic form aids the model in learning the original PDE.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)