MultiPDENet: PDE-embedded Learning with Multi-time-stepping for Accelerated Flow Simulation
Authors: Qi Wang
Yuan Mi
Wang Haoyun
Yi Zhang
Ruizhi Chengze
Hongsheng Liu
Ji-Rong Wen
Hao Sun
Venue: ICML 2025 poster (2025)
Keywords: physics-informed learning
multiscale time stepping
spatiotemporal dynamics
PDEs
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: Introduced a PDE-embedded learning model with multi-time-stepping for accelerated simulation of flows.

Abstract: Solving partial differential equations (PDEs) by numerical methods meet computational cost challenge for getting the accurate solution since fine grids and small time steps are required. Machine learning can accelerate this process, but struggle with weak generalizability, interpretability, and data dependency, as well as suffer in long-term prediction. To this end, we propose a PDE-embedded network with multiscale time stepping (MultiPDENet), which fuses the scheme of numerical methods and machine learning, for accelerated simulation of  flows. In particular, we design a convolutional filter based on the structure of finite difference stencils with a small number of parameters to optimize, which estimates the equivalent form of spatial derivative on a coarse grid to minimize the equation's residual. A Physics Block with a 4th-order Runge-Kutta integrator at the fine time scale is established that embeds the structure of PDEs to guide the prediction. To alleviate the curse of temporal error accumulation in long-term prediction, we introduce a multiscale time integration approach, where a neural network is used to correct the prediction error at a coarse time scale. Experiments across various PDE systems, including the Navier-Stokes equations, demonstrate that MultiPDENet can accurately predict long-term spatiotemporal dynamics, even given small and incomplete training data, e.g., spatiotemporally down-sampled datasets. MultiPDENet achieves the state-of-the-art performance compared with other neural baseline models, also with clear speedup compared to classical numerical methods

MultiPDENet:
PDE-embedded Learning with Multi-iime-steppnng for Accelerated Flow Simulation
Qi Wang'
Yuan Mi!
Haoyun Wang' Yi Zhang² Ruizhi Chengze?
Hongsheng Liu?
Ji-Rong Wen Hao Sun
Abstract
1. Introduction
Solving partial differential equations (PDEs) by numerical methods meet computational cost chal- lenge for getting the accurate solution since fine grids and small time steps are required. Machine learning can accelerate this process, but strug- gle with weak generalizability, interpretability, and data dependency, as well as suffer in long- term prediction. To this end, we propose a PDE- embedded network with multiscale time stepping (MultiPDENet), which fuses the scheme of nu- merical methods and machine learning, for ac- celerated simulation of flows. In particular, we design convolutional filter based on the structure of finite difference stencils with a small number of parameters to optimize, which estimates the equiv- alent form of spatial derivative on a coarse grid to minimize the equation's residual. Physics Block with 4th-order Runge-Kutta integrator at the fine time scale is established that embeds the structure of PDEs to guide the prediction. To alleviate the curse of temporal error accumula- tion in long-term prediction, we introduce mul- tiscale time integration approach, where a neural network is used to correct the prediction error at a coarse time scale. Experiments across vari- ous PDE systems, including the Navier-Stoker equations, demonstrate that MultiPDENet accurately predict long-term spatiotemporal dy- namics, even given small and incomplete train- ing data, c.g., spatiotemporally down-sammn sampls datasets. MultiPDENet achieves the state-of-the- art performance compared with other neural base- line models, also with clear speedup compared classical numerical methods.
Complex spatiotemporal dynamical systems, e.g.. climate system (Schneider et al., 2017) and fluid dynamics (Ferziger et al., 2019), are fundamentally governed by partial differ- ential equations (PDEs). To capture the intricate behaviors of these systems, various numerical methods have been de- veloped. Direct Numerical Simulation (DNS) is a widely used method for solving PDEs. requires specifying initial conditions (ICs), boundary conditions (BCs), and PDE pa- rameters, followed by discretizing the equations using techniques like finite difference (FD), finite element (FB), finite volume (FV), or spectral methods. Despite their accuracy, traditional numerical methods face key challenges of high computational costs (Goc et al., 2021), when ad- dressing with high-dimensional problems or necessitating fine spatial and temporal resolutions.
grid
Recent advances in deep learning have introduced neural- based approaches (Lu et al., 2021; Li et al., 2021; Gupta & Brandstetter, 2023) for solving PDEs. These data-driven methods eliminate the need for explicit theoretical formula- tions, enabling networks to learn underlying patterns directly from data through end-to-end training. While promising, these approaches face notable challenges, including heavy dependence on large training datasets and limited generaliza- tion. For instance, achieving accurate predictions becomes particularly challenging when models encounter unseen ICs or scenarios beyond the training distribution.
A notable progress lies in the paradigm of Physics- informed neural networks (PINNs) (Raissi et al., 2019), which incor- porates physical prior knowledge (such as PDE residuals and 1/BCs) as constraints within the loss function. This approach allows the network to fit the data while simultane- ously maintaining a certain degree of physical consistency. Variants of PINNs (Raissi et al., 2020; Wang et al., 2020; Eshkofti & Hosseini, 2023) have shown notable success across various domains, reducing the dependency on exten- sive datasets to some degree. However, such methods still face scalability and generalizability challenges when applied to complex nonlinear dynamical systems. Additionally, op- timizing complex loss functions (Rathore et al., 2024) and ensuring model interpretability remain challenges.
'Gaoling School of Artificial Intelligence, Renmin University of China, Beijing, China HHuawei Technologics, Shenzhen, China. Correspondence to: Hao Sun <haosun ruc.edu.en>.
Proceedings of the ndnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

PDE-embedded Learning with Multi-time-stepping
A series of approaches have been proposed to integrate physics into neural networks (NNs) overcome the above challenges. For instance, PeRCNN (Rao et al., 2022; 2023), which uses feature map multiplication to construct polyno- mial combinations for approximating the underlying PDEs, can capture the latent spatiotemporal dynamics even with low-resolution, noisy, and coarse data, demonstrating strong generalizability. Nevertheless, this method suffers from er- accumulation, degrading its performance in long-term predictions. Another approach (Kochkov et al., 2021; Sun et al., 2023), combining NNs with numerical methods, aims to accelerate the simulation process on coarse grids. These hybrid methods leverage traditional solvers for stability and NNs for accuracy. However, they often rely heavily on NN capabilities and often requires large amounts of data,
where u(x,() E IR" denotes the physical state in the spa- tiotemporal domain 92x [0, T); u, the first-order time deriva- tive term; F(-) linear/nonlinear functional parameterized by PDE parameters À (e.g., the Reynolds number Re); V the Nabla operator is defined as [Ox, ay,..."; and the source term. Additionally, we define I(,, U3 x E Q2, = 0) = 0 and Bu, Vu, u..... ; E 02) = specified ICs and BCs, where 0S represents the domain boundary.
We aim to accelerate the simulation of fluid flows by using PDE-embedded network with multiscale time stepping based on a limited training data (coarse in both spatial and temporal scales). The model is capable of rapid simula- tion, achieving high solution accuracy while demonstrating strong generalizability across varying ICs, source terms, complex domains, and PDE parameters.
To overcome these limitations, we propose MultiPDENet, PDE-embedded network that incorporates multiscale time- stepping (as shown in Figure 1). to efficiently simulate spa- tiotemporal dynamics, e.g., turbulent fluid flows, on coarse spatial and temporal grids with limited data. Notably, it integrates a trainable neural solver for precise predictions at micro time scales, while employing NN correct errors at macro time steps. Additionally, by embedding PDEs, MultiPDENet offers enhanced generalizability. The primary contributions of this work are summarized as follows:
3.2. Model Architecture
In this section, we introduce MultiPDENet and show how our model efficiently captures the underlying spatiotem- poral dynamics. As illustrated in Figure ((a), predicting Us+1 from the input us involves two main components: the Physics Block and the M,NN Block.
3.2.1. MULTI-SCALE FORWARD TIME STEPPING SCIIEME
While the learnable neural solver can be used independently, its accuracy for long-term prediction is limited due to error accumulation. To address this issue, we introduce multi- scale time stepping scheme, incorporating micro-scale and macro-scale steps, to improve predictive accuracy and en- ables fast prediction of PDE solutions on coarse spatiotem- poral grids. Specifically, we define two types of time step- ping: micro-scale step and macro-scale step, to enhance the performance of spatiotemporal dynamics prediction. At the maero scale, given the coarse solution u, at time tk- MultiPDENet is expected predict the next-step solution Uk+1 at 44+1, which can be expressed as:
We developed MultiPDENet, a PDE-embedded net- work with multiscale time-stepping, for accelerated flow simulations on spatiotemporal coarse grids. By integrating neural solver with PDEs, MultiPDENet achieves great generalizability and efficiency.
Leveraging the structure of FD stencils, we introduced a symmetric convolutional filter that approximates the equivalent form of derivatives on coarse grids, aiming to reduce the residual error of the governing PDEs.
Experimental results across various datasets, covering ID and 2D equations (e.g., reaction-diffusion processes and turbulent flows), demonstrate the effecctivenes of MultiDENNet in accelerating long-term simulations.
M Σ
u&+1 Us + Su + M,NN(UE. At, dz).
(2)
2. Related Work
where At the macro-scale time interval, and de the spatial resolution of mesh grid. Here, oü is the incremental up- date by the Physics Block (see Section 3.2.2) at each micro step, as shown in Eq. (3), where M denotes the number of micro-scale time steps in one macro-scale step (e.g., M- = 4 herein). The M,NN Block (see Section 3.2.4) refines these incremental updates generated by the Physics Block on coarse grids, yielding the final update for the macro step.
Related works on numerical, machine learning, physies- inspired learning, and hybrid learning methods for simula- tion of PDE systems given in Appendix A.
3. Methodology
3.1. Problem Deseription
Let's consider a general spatiotemporal dynamical system governed by the following PDE:
3.2.2. PHYSICS BLOCK: LEARNABLE NEURAL SOLVER
To accurately predict at the micro-scale step, we developed a neural solver, referred to as the Physics Block, as illustrated
u. F(u, u2 ..... Vu, Au....
+f,
2

---

PDE-embedded Learning with Multt-time-stepping
PDE Block
60
60
60
Correction
t
Block
Physics Block
Physics Block
Physics Black
L1g
up++1
Micro Step
M,NN Block
Filter Bank
Block
b
Physics Block
81
Collect derivatives
Poisson
Correction Block
POE Block
Residue Block
B(-)
f(T)
Axt (m 136t
mm-1
Figure 1. Schematic of MultiPDENet for learning turbulent flows. (a), Model architecture. (b), Physics Block. (c), Learnable PDE block.
in Figure 1(b). This solver designed ensure the stability (Hoffman & Frankel, 2018), accuracy, and efficiency of its predictions by adhering to the Courant-Friedrichss (CFL) conditions (LeVeque, 2007). The Physics Block comprises three main components: the Poisson Block, the PDE Block, and the M,NN Block. The solution update for each micro-scale time step can be describe as t u + Su where
symmetry constraints, which calculates derivative terms based on the corrected solution produced Correction Block. These terms then combined into the governing PDEs, a learnable form of JF in Eq. (1). This process is incorporated into the RK4 integrator (see Appendix B.3) for solution update which can be expressed as
F ( approz:
k m
A
m
(4)
t = f
pla tmåt
B (
Vuk
V A).
B(ū(T), vu(T),.... TM -;))+((T)))))))))))))))))))))))) dr
Ja+(m=-166 1)5t
where B denotes the PDE block, and i the coarse solution (aka, solution on coarse grids) at micro-scale time ta + môt. Here, ú refers to the neural-corrected state of the coarse solution, which obtained through the Correction Block (see Appendix B.I for details). This corrected state û is used estimate spatial derivatives, namely, út NN(
M,NN u =
p, ú, Vp, f. Re))
(3)
Here, u represents the intermediate state at mn-th micro- scale step initialized time l (note that of = ug). We denote û(r) A u(x, r), where depicts the coordinates of coarse grid. Moreover, À can be set trainable if unknown. B represents the PDE block, used for approximating F in Eq. (I). To keep the accuracy and ensure the stability, the PDE block is designed based on the RK4 integrator (see Appendix B.3) and consists of the Correction Block and trainable filter bank. Since the considered micro-scale time interval relatively large, the M,NN Block is used as corrector refine the solution. More details can be found in Appendix B. In fact, the Physies Block can be used for prediction independently (e.g., the quantitative results for the NSE dataset predictions using purely the Physics Block are presented Table 3, labeled Model C).
Note that V and A represent trainable Nabla and Laplace operators, respectively, each consisting of symmetrically constrained convolution filter, e.g., an enhanced FD kernel to approximate spatial equivalent derivatives. By utilizing the RK4 integrator, we can project the coarse solution to the subsequent micro-scale time step. Despite the reduced res- olution causing some information loss, this learnable PDE block enables closer approximation of the equivalent form of the derivatives on coarse grids. This addition serves as fully interpretable "white box" element within the overall network structure.
Poisson Block. In solving incompressible NSE, the pres- sure term, p. is obtained by solving an associated Pois- son equation. To compute the pressure field, we imple-
PDE Block. The PDE block computes the residual of the governing PDEs. It incorporates a learnable filter bank with
3

---

PDE-embedded Learning with Multi-time-stepping
Table Overview of datasets and training configurations. Note that " denotes the downsampling process from the high resolution (simulation) to the low resolution (training and testing).
Dataset
Numerical
Spatial Grid
Time Steps
of Training Trajectories
of Testing Trajectories
Macro-step Rollout
Micro-step Rollout
Method
(Temporal Grid)
KdV
Spectral
256 + 64
10000 - 2000
3
10
10
Burgers
FD
100 25
2000 + 200
5
10
10
GS
FD
128 + 32
4000 - 200
3
10
1
4
NSE
FV
2048² + 64
153600 + 1200
5
10
I
mented pressure-solving module shown in Figure Sl(a), which solves the Poisson equation, Ap = (u)), where (b( = 2 (Ma - Un) for 2D problems (the subscripts indicate the spatial derivatives along x or directions). To compute the pressure, we employ a spectral method (Pois- son solver) based on (( to calculate Pmn As shown in Figure S1()), this approach dynamically estimates the pressure field from the velocity inputs, removing the need for labeled pressure data.
1st-order
2nd-order
derivative filter
derivative filter
as
4g
0
as
Ag
0
G
as
Q4
0
0g
04
0
Q4
03
g
a1
42
0
a2
41
g
Q1
02
02
01
ag
a4
0
4g
ag
04
0
04
83
0g
0
06
a5
0
46
as
3.2.3. ADAPTIVE FILTER WITH CONSTRAINT
Figure 2. Symmetric filter
Traditional FD methods often yield inaccurate derivatives on coarse grids. To address this, we propose learnable filter with constraints that approximates equivalent derivatives on coarse grids, minimizing the PDE residuals during training and thereby improving the model's predictive accuracy. By leveraging the symmetry of central difference stencils, our filter maintains structural integrity while enhancing network flexibility. As shown in Figure 2, we construct two 5 x symmetric matrices, each requiring only six learnable parameters due to symmetry constraints. These matrices designed to compute the first-order (g') and second- order (g") derivatives, respectively. In the matrix of g", 4 x (ag +a aa +ag) + (a, +a2). This design leverages the structural properties of central difference methods. By satisfying the Order of Sum Rules (Long et al., 2018), this filter can achieve up to 4th-order accuracy in approximating the derivatives via optimization of trainable parameters.
shown in Figure 1(b) in the upper path, u is first corrected by the Conrection Block, and f.h computed by the Poisson Block. Inputs, including solution states {a p and their derivative terms, forcing term, and Reynolds number, fed into the M,NN Block (see Figure S1(d)). The M,NN Block continuously refines the PDE block's ' outputs on the fly. For detailed information of the M,NN Block settings, please refer to Appendix Table S5.
M,NN Block. Although the Physics Block offers real-time corrections for the M,NN outputs, errors still accumulate in long-term predictions. To mitigate error accumulation in long-term predictions given training data sampled at large time steps (e.g., 128 for the NSE dataset), we introduce the M,NN Block. As depicted in Figure (aa), the M,NN Block takes the current velocity field u, as input, and up- dates the solution u which obtained by integrating the outputs from both the upper and lower paths. During the backpropagation, the M,NN Block learns to correct the coarse solution output of the Physics Block in real time, ensuring that their combined results more closely align with the ground truth. The configuration details for this block are found in Appendix Table S6.
3.2.4. NN BLOCK
To alleviate the error accumulation during long-term predic- tions on coarse grids, we introduce the M,NN and M,NN Blocks, operating at micro- and macro-scales, respectively. The MiNN Block employs a lightweight model (c.g., FNO, DenseCNN (Liu et al., 2024a)) for efficient micro-step pre- dictions, whereas the M,NN Block delivers more accurate predictions larger steps (Gupta & Brandstetter, 2023). In this study, utilized FNO the M,NN Block and UNet as the N,NN Block. The significance of these blocks is evident from the ablation studies presented in Table 3.
4. Experiment
We validate the performance of our method against baseline models on various PDE datasets. We then perform general- ization tests across different external forces (f), Reynolds numbers (Re), and domain sizes on the Kolmogorov flow (KF) dataset. Finally, we present ablation studies to demon- strate the contributions of each component in our model.
M,NN Block. The M,NN block is designed to rectify error accumulation during micro-scale time step predictions. As
4

---

PDE-embedded Learning with Multi-time-teepping
Table Results of MulGPDENet and baselines. For KdV, Burgers, and GS, we inferred upper time limits of 50 1.4 s, and 1200 8, for the test set as the system dynamics stabilized within these trajectories. These time limits were used to calculate HCT.
ing ground truth. Here, N denotes the number of batches, B the batch size, and X the trainable PDE parameters.
Model generalization. The generalization of MultiPDENet is evaluated across ICs, PDE parameters (c.g., Re), force terms, and computational domain sizes (e.g., different mesh grids). The model integrates ICs through its time-marching mechanism, ensuring robust generalization when trained
Case
Model
RMSB (Q)
MAE (1)
MNAD ())
HCT (8)
FNO
0.9541
0.4607
0,3469
10.0833
PINO
0.4120
0.3022
0.2139
13.90
UNet
19887
1,5222
1.6158
3,1250
Kdv
DeepONet
NaN
NaN
NaN
0.1500
effectively. The Reynolds number (Re) represented via (a & b), using trainable vectors and b. This embedding. applied in both the PDE and MiNN blocks, reduces error propaga- tion from the diffusion term on coarse grids and enhances generalization across Re values. The force term is incor- porated into the learnable PDE block and the MiNN block, where serves as both a PDE feature and an input feature
MultiPDENet
0.1536
0.1110
0.0833
39.8
wwo-dimensional embedding, RCcmbb
Improveeent (1)
62.7%
63.3%
61.1%
186.3% 3%
FNO
0,0980
0,0762
0.0620
0,3000
PINO
0.0832
0.0749
0.0599
0.5546
UNet
0.3316
0.2942
0.2556
0,0990
Burgers
DeepONet
0.2522
0,21066
0.1692
0.0020
PeRCNN
0.0967
0.1828
0.1875
0.4492
MultiPDNNet
0.0057
0.0037
0.0031
1.4000
Improvement (1)
93.1%
95.1%
94.8%
152.4%
FNO
8774
1303
1303
270
PINO
0.5721
0.3579
0.3520
510
map shown in Figure S1, enabling joint learning of force variations for better generalization.
UNel
NaN
NaN
NaN
20
GS
NeepONet
0.4113
0.2961
0,2898
568
PeRCNN
0.1763
0.1198
0.1198
640
Evaluation metries. We evaluate the performance of our model using four metrics: Mean Absolute Error (MAE),
MultiPDENet
0.0573
0.0294
0.0298
1400.0
Improvement (t)
67.5%
75.5%
75.1%
118.8%
Root Mcan Square Error (RMSE), Mean Normalized Abso-
FNO
L0100
0.7319
0.0887
2.5749
UNet
0.8224
0,5209
0.0627
3.9627
lute Difference (MNAD), and High Correlation Time (HCT).
LI
NaN
3,5000
NSE
TSM
NaN
NaN
NaN
3.7531
For detailed definitions, please refer to Appendix E.
DeepONet
2.1849
1.0227
0.1074
0.1126
Baseline models. To ensure a comprehensive comparison, we selected several baseline models, including FNO (Li et al., 2021), PINO (Li et al., 2024c), UNet (Gupta & Brand- stetter, 2023). TSM (Sun et al., 2023), LI (Kochkov et al., 2021), DeepONet (Lu et al., 2021), and PeRCNN (Rao et al., 2023). Details are found in Appendix F.
MultiDDNNet
0.1379
0.0648
0.0077
8.3566
Improvement (1)
83.2%
87.6%
87.7%
110.9%
4.1. Setup
Dataset. We generate the data using high-order FD/FV methods with high resolution under periodic boundary con- ditions and then downsample it spatially and temporally to a coarse grid. The low-resolution dataset is used for both training and testing. We consider four distinct dynamical systems: Korteweg-de Vries (KdV), Burgers, Gray-Scon (GS), and Navier-Stokee equations (NSE). Each dataset is divided into 90% for training and 10% for validation. We segment trajectories into data series, where each sample in- cludes multi snapshots (e.g., for the KdV dataset, the sample length is set to 10, as detailed in Table 1) separated by a time step At, the 2nd to the last snapshot serves as the training labels. During training, we use only 3-5 trajectories for each system, and evaluate them on 10 distinct trajectories. For further details, please refer to Appendix C.
4.2. Solving PDE Systems
Kdv. The primary challenge of this dataset lies in accu- rately capturing the complex interplay between nonlinearity and dispersion, leading to phenomena like soliton formation (Gardner et al., 1967). As shown in Figure 3()), each base- line model struggles to produce accurate predictions, with DeepONet exhibiting significant divergence. In contrast, MultiPDENet demonstrates superior accurate predictions for ICs outside the training range. The correlation curve in Figure 3(b) highlights the significantly higher correla- tion of MultiPDENet compared to the baselines. The error distribution in Figure 3(e) confirms its lower error levels. Ta- ble 2 shows our model's generalizability, with performance improvements ranging from 61.1% to 186.3%.
Model training. Our objective is to accelerate flow sim- ulations with all computations anchored to coarse grids. During training, the model solely predicts the solutions for subsequent time steps, employing Mean Squared Error (MSE) as the loss metric. Unlike PINNs, our MultiPDENet directly embeds PDEs into its architecture, resulting in loss function that exclusively comprises data loss, given by: N BN Ein Σ MSE (H.y. H.y), where H.j
Burgers. As shown in Figure 3()), the solution snapshots predicted by MultiPDENet are significantly more accurate than those of the bascline models. The bascline models, limited by the coarse training data, produce incorrect pre- dictions. The correlation curve in Figure 3(e) shows that MultiPDENet maintains a high correlation with the ground truth throughout the prediction, while other baseline models diverge. This is further evidenced by the error distribution in Figure 3(), demonstrating that MultiPDENet's error is
J(A) = 1 denotes the coarse solution predicted by model rollout for the j-th sample in the i-th batch, and H.y is the correspond-
5

---

PDE-embedded Learning with Multi-time-stepping
Ref.
MuleDENee
PINO
FNO
DeepONet
UNet
Simulation time
30
50
Error distibution
1.0
10"
0.8
Corcelation
10-1
0.6
Time
10
0.4
200
MulliPDENet
Ulet
NaN
0.2
10
DeepONee
FNG
PINO
Spatial (m]
0.0
200
400
Time steps
800
800
1000
10
Models
d
Simulation time
Ref.
DeepONet
0.18
1.0
0
2
6
6
10
12
10
Enror distribution
Conelation
0.00
MNAD
10
10
PUEE
10"
DeegOONet
IND
PeRCNN
0.25
0.0
60
120
140
Snapshots timestp
staps
Models
g
Simulation time
Ref.
MuAPDEENe
FND
UNit
PeRCNN
DeepONer
1.0
1400
109
Enar distrbution
0.60
0.8
10"
Condatio
0.6
MNAD
10
0.61
0.4
MultPDENct
UNet
0
0.29
0.2
PaRCNN
FND
0.0
20
40
60
80 100 120
140
10
Snepshots timestep
Time steps
Figure overview of the comparison between MultiPDENet and baselines, including predicted solutions (left), correlation curve (middle), and error distributions (right). (a)-(c) show the qualitative results on KdV. (d)-(() show the qualitative results on Burgers. (g)-(i) show the qualitative results on GS. These PDE systems are trained with grid sizes of 64, 25, 32x32, respectively.
over an order of magnitude lower than that of the baselines. Table confirms these findings with MultiDDNNet's im- provements exceeding 94. % across evaluation metries.
starts to produce incorrect patterns at t = 5 s due to error accumulation. The correlation curve Figure 4(b) supports these findings. Our model also achieves a spectrum energy distribution elosely matching the ground truth (see Figure 4(c)). Table 2 highlights performance improvement of over 83.2%. Even with 20% less training data, our model maintains strong generalizability (see Appendix Table S3).
GS. This rraction-diffusion system is nonlinear, making it challenging to capture its complex patterns (see Figure 3(g)). Only MultPDENeet accurately predicts the trajectory evolu- tion. The baseline models struggle learn the spatiotem- poral dynamics, and even PeRCNN, despite its embedded physics, produces inaccurate predictions due to the limited and coarse training data. Figure 3h) demonstrates the su- perior correlation of MultiPDENet's predictions with the ground truth. The error analysis in Figure 3(i) reveals signif- icantly lower error levels for MultPDENet. often by to orders of magnitude smaller compared to the baselines. Ta- ble 2 further validates this observation, with MultiPDENet's improvements of 67.5 37.5% 118.8% over the best baseline.
4.3. Model Generalization
We conducted generalization tests on the KF flow dataset assess model's ability capture the underlying dynam- ics. The model was initially trained on 5 sets of trajectories, where the forcing term is defined as f sin(4y)n, 0.1u with Re = 1000 and n, [1, 0j"" After training, we tested MultiPDENet on 10 different sets of trajectories, each with varying Reynolds numbers (Re) and forcing terms (f), to evaluate its performance across a range of different ICs.
NSE. We evaluate KF with Re 1000 across different ICs, governed by the NSE. Figure 4(a) shows the trajectory snapshots predicted by MultiPDENet and the baseline mod- els over 10 s. Our model outperforms DNS 512, accurately capturing both global and local correction patterns. The neural methods, particularly FNO, exhibit poor generaliza- tion, producing granular and erroneous solutions. Among the Physics ML baselines, TSM performs the best, but
Test on various Reynolds numbers. Firstly, we evaluate the generalizability MultiPDENet across four different Reynolds numbers: Re 500, 800, 1600, 2000. The vary- ing Re values result in trajectories with differing levels of complexity. Figure 5()) displays the accurate predictions made by our model at time step 300 for different Re. Fig- ure 5(b) shows the correlation curves between the predicted
6

---

PDE-embedded Learning with Multi-time-seepping
b
Simulation time
DNS 2048
2
12
14
DeepONet
TSM
Numerical Solver
DNS 64
DNS
DN5 128
Correlaton
0.6
DNS 256
DNS 66
5
DNS 1024
0.4
MuRtPOEENet
MuluPOENet
0.2
0.0
Vorticity
125
250
375
500
Time
Phy ML
1033
DeepONet
DeegONet
TSM
UNet
1012
FNO
Scaled energy sprectuum E08
DNS 128
TSM
DNS 256
DNS 512
10
DNS 1024
MultiPDENet
UNied
Data driven
1020
FNO
0.5
10a
10P
302
Wavenumber k
Figure 4. Comparison of MuliPDENet and baseline models on Kolmogorov flow with Re = 1000. (a) shows the evolution of predicted vorticity fields for reference, MultiPDENet baselines, starting from the same initial velocities. (b) shows the correlation curve across 500 time steps. (c) shows the scaled energy spectrum scaled by & averaged between time steps 100 and 500.
and ground truth trajectories, while Figure 5(c) highlights the error distributions, which remain consistently below 0.1, indicating a low error level.
training, the model was tested 10 trajectories with new ICs. Appendix Figure 6(a) illustrates the snapshots pre- dicted by MultiPDENet over 600 timesteps. demonstrating sustained accuracy even at time step 450. The correlation curve in Appendix Figure 6(b) highlights the superiority of our model compared DNS 1024. The error analysis in Appendix Figure 6(c) confirms this performance, with er- consistently below 0.01. These results demonstrate the effectiveness of MultiPDENet for higher Reynolds number, e.g., Re 4000, within domain (0, )x.
Test external forces.
Next, performed the generaliza- tion test using distinct f. As shown in Figure 5(d), Mul- tiDENNet accurately predicts the trajectories for all forces. Notably, the downward trend in the correlation curve (in Figure (e)) for fo parallels the behavior observed in Fig- ure 4(b), likely because f alters only the periodic function without changing the wave number. The error analysis in Figure 5(l) confirms that the error levels remain below 0.1.
Test on flow within larger domains. We extended the spa- tial domain from (0, 2 to 4n) to further evaluate our model's generalizability over larger mesh grids in a more complex scenario. Larger domains introduce diverse physi- cal phenomena, challenging the model capture global and local dynamics coarse grids. Using the same 64 x 64 grid, we tested our trained model on 10 unseen trajectories
Test on flow with Re 4000. Turbulence at high Re's presents significant challenges for prediction due to its non- linearity and complex vortex structures. To further demon- strate the superior capability our model, we conducted an additional experiment with a high Re 4000 (see details in Table 82) with the experimental setup in Section 4.1. After

---

PDE-embedded Learning with Multi-time-stepping
b
Simulation time
500
10
1.0
10
10c
Error distribution
Ro
MuliPDENet Rel.
0.8
10-
Correlston
0.6
MNAD
10-2
DNS
0.4
Re 500
10"2
Rc 800
0.2
10-4
Re 2000
-10
0.0
150
300
450
600
10"5
Snapshots at timestep 300
Time steps
Reynolds numbers
d
Simulation time
f
f2
fg
f4
Error distribution
10
1.0
109
f:
Ref.
0.8
10-
Correlatin
0.6
MNAD
10-2
DNS
MultiiDNNe
0.4
10"3
0.2
10-4
-10
0.0
150
300
450
600
10"
Snapshots at timestep 300
Time steps
Extemaal forcings
Figure 5. MultiPDENet can generalize to simulate different Reynolds numbers and external forcings without retraining. Vorticity snapshots predicted by MultiPDENct and ground truth timestep 300 (left), correlation curve 600 timesteps (middlc), error distribution (right). (a-c) display results for different Reynolds numbers. (d-f) show results for varying external forcings.
a
b
Simulation dime
Time stop = Time stop 150 Time soep = 300 Time stop 450 Time stop = 500
1.0
2
4
109
DNS 4096
0.8
Constatio
0.6
DNS 2048
MNAD
10
DNS 1024
MuiSPDENet
0.4
10
DN5 256
0.2
DNS 128
MutiPDENe:
256
0.0
200
400
600
10
Variiity snaps hots over 600 tinsstoss
Time steps
Models
Figure 6. MultiPDENet is applicable to high Reynolds number turbulence. (a) Trajectorics predicted by MultiPDENct at R (b-c) Correlation error distribution comparison between MulGPDENet and numerical method.
(details in Appendix C). As shown in Appendix Figure 7(a), the snapshots over 300 time steps remain accurate. The cor- relation curve in Appendix Figure 7(b) depicts our model's performance closely matches or exceeds DNS 1024. The error distribution in Appendix Figure 7(C) shows error levels comparable to DNS 1024. Notably, our model achieves speedup of 6x compared to the FV DNS method.
Poisson Block); (2) Model (no filter structure constraint);
(3) Model (only Physics Block for prediction); (4) Model
D (FD convolution instead of symmetric filter); (5) Model
E (no Correction Block); (6) Model F (no M,NN Block);
(7) Model G (no M,NN Block); (8) Model H (no Physics Block); (9) Model I (forward Euler); and (10) MultiPDENet (the full model). The results are shown in Table 3.
In summary, MultiDENett demonstrates remarkable gener- alizability, showcasing its ability to capture the underlying dynamics across multiple temporal scales. The embedded learnable PDE module within model is crucial for abling robust and accurate predictions.
Removing the Poisson Block impairs our model's perfor- mance, confirming p-u decoupling necessity the NSE, Relaxing the filter structure constraint degrades the result, validating our proposed kernel's l's efficacy. The Physics Block alone yields worse prediction compared to the full model. Replacing the Conv kernel with FD stencils results in poorer performance, indicating that fixed-value FD kernels are un- suitable for coarse grids. Omitting the Correction Block also degrades the model prediction, highlighting the necessity of field correction. While the model can still accurately predict up to 5.9 without the N,NN Block, the error increases by
4.4. Ablation Study
To quantify the contribution of each module, conducted ablation experiments the KF dataset. Specifically, we compared the following model variations: (1) Model A (no
8

---

PDE-embedded Learning with Multi-time-teepping
a
Simulation
ONS 4080
Time stop Time step
stop
15
1.0
12 14
10a
Enrr distribution
0334
10-1
DMD 312
0.6
MNAD
10-
MURPEENS
10"
MutIPDENcK
DNS 512
0.2
10-
CINS 120
1024 2048
DNS 256
:15
200
300 steps
500
10-2
Vortchty snopshats over 300 timestops
Models
Figure 7. MuliiDENet applied to a larger domain. (a) Predicted trajectories within 4 x 4m. (b-c) Correlation and error distribution comparison with the numerical method.
Table Results of the ablation study.
ity in long-term prediction of turbulent flows, effectively
capturing both global and local patterns in larger domains.
Ablated Model
RMSE (L)
MAE (L)
MNAD (J)
HCT (t)
Model A
0.1601
0.0711
0.0085
7.904
We also tested the computational efficiency of trained Multi- PDENet for accelerated flow prediction (more details shown in Appendix G.2). For a certain given accuracy (e.g., corre- lation 2 0.8), MultiPDENet achieves 2 5x speedup com- pared with GPU-accelerated DNS (Appendix Table S8), e.g., JAX-CFD, where all the tests were performed on a single Nvidia A A100 S0G GPU. However, MultiPDENet still faces some unresolved challenges. Firstly, the model currently only handles regular grids, due to the limitation of convolu- tion operation used in the model. In the future, we aim to address this issue by incorporating graph neural networks to manage irregular grids. Secondly, the model has only been currently tested on ID and 2D problems. We will extend to 3D systems in our future work.
Model B
0.2432
0.1156
0.0137
7,8633
Model C
0.2632
0.1402
0.0186
7.3146
Model D
0.2503
0,1410
0.0145
7.0783
Model E
0.2958
0,1453
0,0180
6.6856
Model F
0,4338
0.2401
0.0285
5,9768
Model G
NaN
NaN
NaN
1.4193
Model
1,2023
0.9256
0.1122
0.6227
Model
0.4357
0,2321
0.0278
6.2481
MultPDDNeet
0.1379
0.0648
0.0077
8.3566
4x, suggesting the need for the M,NN Block.
Removing the N,NN Block restricts the micro-scale accu- racy and causes macro-scale instability due to eπror accu- mulation. This underscores the importance of the M,NN Block for long-term stability. The omission of the Physics Block (e.g., solely U-Net) cripples the prediction with lim- ited coarse training data., proving its critical role. Moreover, the use of RK4 outperforms forward Euler in terms of sta- bility and accuracy. Hence, all the components are essential and contribute meaningfully to the MultiPDENet model.
A Acknowledgment
The work is supported by the National Science and Technology Major Project (No. 2022ZD0117804), the National Natural Science Foundation of China (No. 62276269 and No. 92270118), the Beijing Natural Sci- ence Foundation (No. 1232009), and the Fundamen- tal Research Funds for the Central Universities (No. 202230265). Our source codes are available in the fol- lowing GitHub repository: https: http//ttttbu.co// com/ intell-sct comput/Mull comput/MMaltiPDENNe iPDNNet.
5. Conclusion
We introduce an end-to-end physics-encodeee network (aka, MultiPDENet) with multi-scale time stepping for acceler- ated simulation of spatiotemporal dynamics such turbu- lent flows. MultiPDENet consists of multi scale temporal learning architecture, a learnable Physics Block for solution prediction at the fine time scale, where trainable symmetric filters are designed for improved derivative approximation on coarse spatial grids. Such a method is capable of long- term prediction on coarse grids given very limited training data (see the data size scaling test in Appendix D.4). Mul- tiPDENet outperforms other baselines through extensive tests on fluid dynamics and reaction-diffusinn equations. In particular, such a model excels in generalizability over ICs, Reynolds numbers, and external forces in the turbulent flow experiments. MultiPDENet also exhibits strong stabil-
Impact statement
This work aims to develop PDE-embedded network with multiscale time stepping, integrating numerical schemes and machine learning to accelerate flow simulations. This approach has the potential to benefit various fields, including weather forecasting and turbulent flow prediction. This research solely intended for scientific purposes and does not present any foresecable ethical risks.
9

---

PDE-embedded Learning with Multi-time-ttepping
References
Geneva, N. and Zabaras, N. Transformers for modeling physical systems. Neural Networks, 146:272-225 289, 2022.
Arcomano, T., Szunyogh, I., Wikner, A., Pathak, J., Hunt, B. R., and Ott, E. A hybrid approach to atmospheric modeling that combines machine learning with physics- based numerical model. Journal of Advances in Modeling Earth Systems, (403):222MMS000277 2022.
Goc, K. A., Lehmkuhl, O., Park, G. I., Bose, S. T., and Moin, P. Large eddy simulation of aircraft at affordable cost: a milestone in computational fluid dynamics. Flow, 1:E14, 2021.
Azizzadenesheli, K., Kovachki, N., Li, Z., Liu-Schiaffini, M., Kossaifi, J., and Anandkumar, A. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, pp. 1 1-9, 2024.
Gupta, J. K. and Brandstetter, J.
Towards multi- spatiotemporal-scale generalized PDE modeling. Trans- actions on Machine Learning Research, 2023.
He, K., Zhang, X., Ren, S., and Sun, J. Deep residual learn- ing for image recognition. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 770-778, 2016.
Bhatnagar, S., Afshar, Y., Pan, S., Duraisamy, K., and Kaushik, S. Prediction of aerodynamic flow fields using convolutional neural networks. Computational Mechan- ics, 64:525-545, 2019.
Hernández, Q. Badías, A., Chinesta, F., and Cueto, Thermodynamies informed neural networks for physi- cally realistic mixed reality. Computer Methods in Ap- plied Mechanics and Engineering, 407:115912, 2023. ISSN 0045-7825.
Bonev, B., Kurth, T., Hundt, C., Pathak, J., Baust, M., Kashinath. K., and Anandkumar, A. Spherical fourier neural operators: Learning stable dynamics on the sphere. In International Conference Machine Learning, pp- 2806-2823. PMLR, 2023.
Hoffman, J. D. and Frankel, S. Numerical methods for engineers and scientists. CRC press, 2018.
Brandstetter, J., Welling, M., and Worrall, D. E. Lie point symmetry data augmentation for neural pde solvers. pp. 2241-2256, 2022a.
Kochkov, D., Smith, J. A., Alieva, A., Wang, Q. Brenner, M. P., and Hoyer, S. Machine learning- accelerated com- putational fluid dynamics. Proceedings of the National Academy of Sciences, 118(21):e1001884111 2021.
Brandstetter, J., Worrall, D. E.,, and Welling, M. Message passing neural pde solvers. In International Conference on Learning Representations, 2022b.
Kossaczká, T., Ehrhardt, M., and Günther, M. Enhanced fifth order weno shock-capturing schemes with deep learning. Results in Applied Mathematics, 12:100201, 2021.
Cao, S. Choose a transformer: Fourier or galerkin. Advances in Neural Information Processing Systems, 34:24924- 24940, 2021.
Chen, H., Liu, Y. and Sun, H. PINP: Physics-informed neural predictor with latent estimation of fluid flows. In International Conference on Learning Representations, 2025.
Kumar, N. K. A review on burgers' equations and it's appli- cations. Journal of Institute of Science and Technology, 28(2):49-52, 2023.
Dresdner, G., Kochkov, D., Norgaard, P. C.,, Zepeda-Nunez, L., Smith, J., Brenner, M., and Hoyer, S. Learning to correct spectral methods for simulating turbulent flows. Transactions on Machine Learning Research, 2023.
Le Veque, R. J. Finite difference methods for ordinary and partial differential equations: steady-state and time- dependent problems. SIAM, 2007.
Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Fourier neu- ral operator for parametric partial differential equations. In International Conference Learning Representations, 2021.
Eshkofti, K. and Hosseini, S. M. A gradient- enhanced physies-informed neural network (gpinn) scheme for the coupled non-fiikiaa/noonfourrriiaa thrroolsithe analysis: A novel gpinn struc- ture. Engineering Applications of Artificial Intelligence, 126: 126:106908, 2023.
Li, Z., Han, W., Zhang, Y., Fu, Q., Li, J., Qin, L., Dong, R., Sun, H., Deng, Y., and Yang, I.. Learning spatiotemporal dynamics with pretrained generative model. Nature Machine Intelligence, (12))1566-1579 2024a.
Ferziger, J. H., Perić, M., and Street, R. L. Computational methods for fluid dynamics. springer, 2019.
Gardner, C. S., Greene, J. M., Kruskal, M. D., and Miura, R. M. Method for solving the korteweg-devries equation. Physical Review Letters, 1995,) 1967.
Li, Z., Shu, D., and Barati Farimani, A. Scalable trans- former for pde surrogate modeling. Advances in Neural Information Processing Systems, 36, 2024b.
10

---

PDE-embedded Learning with Multi-time-stepping
Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., and Anandkumar, A. Physies- informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1-27, 2024c.
partial differential equations. Journal of Compulational Physics, 378:686-707, 2019.
Raissi, M., Yazdani, A., and Karniadakis, G. E. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science, 367(6481): 2020.
Liu, X.-Y., Zhu, M., Lu, I.., Sun, H., and Wang, J.-X. Multi- resolution partial differential equations preserved learning framework for spatiotemporal dynamics. Communica- fions Physics, 7(1):31, 2024a.
Rao, C.,, Ren, P., Liu, Y., and Sun, H. Discovering nonlinear pdes from scarce data with physics-encoded learning. In International Conference Learning Representations, 2022.
Liu, Y., Hu, T., Zhang, H., Wu, H., Wang, S., Ma, L., and Long, M. itransformer: Inverted transformers are effec- tive for time series forecasting. International conference on machine learning, 2024b.
Rao, C.., Ren, P., Wang, Q., Buyukozturk, O., Sun, H., and Liu, Y. Encoding physics to learn reaction-diffusion processes. Nature Machine Intelligence, 5(7):765-779, 2023.
Liu, Z., Lin, Y., Cao, Y., Hu, H., Wei, Y., Zhang, Z., Lin, S., and Guo, B. Swin transformer: Hierarchical vision transformer using shifted windows. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 10012-10022, 2021.
Rathore, P., Lei, W., Frangella, Z., Lu, L., and Udell, M. Challenges in training pinns: A loss landscape perspec- tive. CoRR. 2024.
Long, Z., Lu, Y., Ma, X., and Dong, B. Pde-net: Learning pdes from data. In International Conference Machine Learning, pp. 3208-320 -2168322,22200 2018
Ren, P., Rao, C., Liu, Y., Ma, Z., Wang. Q., Wang, J.- X., and Sun, H. Physr: Physics-informed deep super- resolution for spatiotemporal data. Journal of Computa- tional Physics, 492:1 12438, 2023.
Long, Z., Lu, Y., and Dong. Pde-net 2.0: Learning pdes from data with a numeric-symbolle hybrid deep network. Journal of Computational Physics, 399: 108925, 2019.
Ren, P., Song, J., Rao, C., Wang, Q., Guo, Y., Sun, H., and Liu, Y. Learning spatiotemporal dynamics from sparse data via a high-order physics-encoded network. Computer Physics Communications, 11::00952,, 2025.
Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218-22) 2291, 2021.
Ronneberger, O., Fischer, P., and Brox, T. U-net: Convolu- tional networks for biomedical image segmentation. In Proceedings of the 18th International Conference Medi- cul Image Computing and Computer-Assisted Interven- tion, part III 18, pp. 234-241, 2015.
Lu, Y., Zhong, A., Li, Q.. and Dong, B. Beyond finite layer neural networks: Bridging deep architectures and numer- ical differential equations. In International Conference Machine Learning, pp. 3276-3285, 2018.
Schneider, T., Teixeira, J., Bretherton, C. S., Brient, F., Pressel, K. G., Schär, C., and Siebesma, A. Climate goals and computing the future of clouds. Nature Climate Change, 7(1):3-5, 2017.
Mi, Y., Ren, P., Xu, H., Liu, H., Wang. Z. Guo, Y., Wen, J.-R., Sun, H., and Liu, Conservation-ing innettioo----oommeddee graph learning for spatiotemporal dynamics prediction. In Pro- ceedings of the 31 st ACM SIGKDD Conference on Knowl- edge Discovery and Data Mining V.1, pp. 1056-1067, 2025.
So, C. C., Li, T. O., Wu, C., and Yung, S. P. Differen- tial spectral normalization (dsn) for PDE discovery. In Proceedings of the AAAI Conference on Artificial Intelli- gence, volume 35, pp. 9675-9684, 2021.
Moukalled, F., Mangani, L., Darwish, M., Moukalled, F., Mangani, L.,, and Darwish, M. The finite volume method. Springer, 2016.
Stachenfeld, K., Fielding, D. B., Kochkov, D., Cranmer, M., Pfaff, T., Godwin, J., Cui, C., Ho, S., Battaglia, P., and Sanchez-Gonzaled A. Learned coarse models for effi- cient turbulence simulation. In International Conference on Learning Representations, 2022.
Rahman, M. A., Ross, Z. E., and Arizzadenesheli, K. U- NO: U-shaped neural operators. Transactions on Machine Learning Research, 2023. ISSN 2835-8856.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physics- informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear
Sun, Z., Yang, Y., and Yoo, S. neural PDE solver with temporal stencil modeling. In International Conference on Machine Learning, pp. 33135-33155, 2023.

---

PDE-embedded Learning with Multi-time-teppinng
Takamoto, M., Praditia, T., Leiteritz, R., MacKinlay, D., Alesiani, F., Pflüger, D., and Niepert, M. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 35: 1596-1611, 2022.
Zeng, D., Liu, W., Chen, W. Zhou, L., Zhang, M., and Qu, H. Substructure aware graph neural networks. In Proceedings of the AAAI Conference Artificial Intelli- gence, volume 37, pp. 11129-11137, 2023.
Zhuang, Kochkov, D., Bar-Sinai, Y., Brenner, M. P., and Hoyer, S. Learned discretizations for passive scalar advection in two-dimensional turbulent flow. Physical Review Fluids, 6(6):064605, 2021.
Tang. K., Zhai, J., Wan, X., and Yang, C. Adversarial adaptive sampling: Unify pinn and optimal transport for the approximation of pdes. In International Conference on Learning Reprexentations, 2024.
Zienkiewicz, O. C.,, Taylor, R. L., and Zhu, Z. The finite element method: its basis and fundamentals. Elsevier, 2005.
Thomas, J. W. Numerical partial differential equations: finite difference methods, volume 22. Springer Science & Business Media, 2013.
Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. Atten- tion is all you need. Advances in Neural Information Processing Systems, 2017.
Wang, Q., Ren, P., Zhou. H., Liu, X.-Y., Deng, Z., Zhang, Y., Chengze, R., Liu, H., Wang, Z.., Wang, J.-X., al. PºC² Net: PDE-Preserved Coarse Correction Network for Efficient Prediction of Spatiotemporal Dynamics. Ad- vances in Neural Information Processing Systems, 37: 68897-68925, 2024.
Wang, R., Kashinath, K., Mustafa, M., Albert, A., and Yu, R. Towards physics-informed deep learning for turbulent flow prediction. Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 1457-1466, 2020.
Wang, R., Walters, R., and Yu, R. Incorporating symmetry into deep dynamics models for improved generalization. In International Conference on Learning Representations, 2021.
Wu, H., Luo, H., Ma, Y., Wang. J., and Long, M. Ropinn: Region optimized physics-infomme informed neural networks. Ad- vances in Neural Information Processing Systems, 2024.
Wu, T., Wang. Q. Zhang, Y., Ying, R., Cao, K., Sosic, R., Jalali, R., Hamam, H., Maucec, M., and Leskovee, J. Learning large-scale subsurface simulations with hybrid graph network simulator. In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 4184 4194, 2022.
Xing, L.. Wu, H., Ma, Wang, J., and Long, M. Helmfluid: Learning helmholtz dynamics for interpretable fluid pre- diction. International Conference Machine Learning, 2024.
Zeng, B., Wang, Q. Yan, M., Liu, Y., Chengze, R., Zhang, Y., Liu, H., Wang, Z., and Sun, H. PhyMPGN: Physics- encoded message passing graph network for spatiotempo- ral PDE systems. In International Conference on Learn- ing Representations, 2025.
12

---

PDE-embedded Learning with Multi-time-stepping
APPENDIX
A. Related work
Numerical Methods. Numerical methods have been extensively applied to solve PDEs. Approaches such as FD (Thomas, 2013), FE (Zienkiewicz et al., 2005), and FV methods (Moukalled et al., 2016) discretize the continuous domain into mesh grids, transforming PDEs into algebraic equations that can be solved with high accuracy. However, these methods often require fine spatiotemporal grids and substantial computational resources to achieve accurate solutions. particularly in high-dimensional spaces. This leads two main challenges: (1) the need for repeated computations when conditions change (e.g., ICs); (2) the demand for fast simulations in many industrial applications.
Machine Learning Methods. Building on the success of machine learning in fields like natural language processing (Vaswani et al., 2017) and computer vision (He et al., 2016), these techniques have also been applied to solving PDEs, With abundant labeled data, it is possible to train end-to-end models to predict solutions. Representative works include ResNet (Lu et al., 2018), CNN-based models (Bhatnagar et al., 2019; Stachenfeld et al., 2022; Gupta & Brandstetter, 2023), Transeermer-based models (Cao, 2021; Geneva Zabaras, 2022; Liu et al., 2024b; Li et al., 2024b), Graph-based models (Brandstetter al., 2022b; Wu et al., 2022; Zeng et al., 2023), and pretrained diffusion models (Li et al., 2024a). Many notable neural operators (Lu et al., 2021; Li et al., 2021; Rahman et al., 2023; Bonev al., 2023; Azizzadeneshelli et al., 2024), which learn mapping between functional spaces, enable the approximation of complex relationships in PDEs. While these methods show promise in learning complex dynamics and approximating solutions, they often require substantial amounts of labeled data.
Physics-inspired Learning Methods. Recently, physics-inspired learning methods have demonstrated impressive ca- pabilities in solving PDEs, which can be classified into two categories according to the way of adding prior knowledge: physics-informed and physics-encoded. The physics-informed methods take PDEs and I/BCs as a part of the loss functions (e.g.. the family of PINN (Raissi et al., 2019; 2020; Wang et al., 2020; Tang et al., 2024; Wu et al., 2024: Chen et al., 2025), PhySR (Ren al., 2023)). On the other hand, the physics-encoded methods employ different approach that preserves the structure of PDEs, ensuring that the model adheres to the given equations to capture the underlying dynamics, e.g., EquNN (Wang et al., 2021), TiGNN (Hermández et al., 2023), PeRCNN (Rao et al., 2022; 2023; Ren et al., 2025), PhyMPGN (Zeng et al., 2025), CiGNN (Mi et al., 2025), and HelmFluid (Xing et al., 2024). In addition, other related studies (Long et al., 2018; 2019; Kossaczká et al., 2021; So et al., 2021) have explored the use of CNN as alternative spatial derivative operators for approximating derivatives and capturing the dynamics of interest.
Hybrid Learning Methods. Hybrid learning methods combine the strengths of numerical approaches and NNs to improve prediction accuracy. For efficient modeling of spatiotemporal dynamics. these methods can be trained on coarse grids. Representative methods include FV-based neural methods (Kochkov al., 2021; Sun et al., 2023), FD-based neural methods (Zhuang et al., 2021; Liu et al., 2024a; Wang al., 2024), and spectral-based neural methods (Dresdner et al., 2023; Arcomano et al., 2022). While these approaches show efficacy in modeling spatiotemporal dynamics, their representation capacities are often limited by the fixed structure of their numerical components. As a result, most of these models still require sufficiently large amounts of training data.
B. The Details of MultiPDENet
B.1. Correction Block
The Correction Block leverages neural network refine the coarse solution, with the Fourier Neural Operator (FNO) (Li et al., 2021) as the correction mechanism within this block. FNO functions by decomposing the input field into frequency components, processing each frequency individually, and reconstructing the modified spectral information back into the physical domain via the Fourier transform. This wise update process is expressed as:
++* = o (W'v'(x) + (xX))))++))))))))))()((((() (X())))
(S1)
where v'(x) denotes the latent feature map at the I-th layer, defined on the coarse grid & The initial feature map is vxx = 3 (uk)) where P is a local mapping function that projects u' into a nigher-dimensions space. The kernel integral transformation is defined as H (z) iFFT(Rs FFF(Z)), which applies the Fourier transform, spectral
13

---

PDE-embedded Learning with Multt-time-stepping
filtering via Rg, convolution in the frequency domain, and the inverse Fourier transform to the latent feature map z. Here, represents the trainable parameters, a-)) is the GELU activation function, and W denotes the weights of the linear layer. After passing through an L-layer FNO, the refined coarse solution is computed as ü (x){ where 2 projects the latent representation of the final layer back into the original solution space.
In the Correction Block, we set I. = 2. For the Burgers equation, we configure the model with modes = 12, width = 12, and a projection from 12 channels to 50 channels. For the GS case, we use the same configuration: modes 12, width 20, with projection from 20 channels 50 channels. For the KdV equation, the setup is defined as modes 32 and width = 64, with projection from 64 channels 128 channels. The NSE case, however, requires a different setup: modes 25, width 20, with projection from channel 20 to channel 128.
B.2. Physics Block
To accurately predict the micro-scale step, developed neural solver called the Physics Block, ensuring stability, accuracy, and efficiency through adherence to the Courant-Friedrichs-Laee (CFL) conditions (LeVeque, 2007). The Physics Block comprises three key components: the Poisson Block (Figure (a)-b)), the PDE Block (Figure I(c)), and the M,NN Block (Figure SI())).
Poisson Solver.
The pressure field is computed using the spectral method, which involves solving the Poisson equation:
Ap = (u)
(S2)
Here, vp(u) = 2 (Ur"ly U2) represents the source term for the pressure.
Applying the Fast Fourier Transform (FFT) to Eq. (S2), we obtain:
- ( + "p)p" = (u),
(S3)
where fhr and fy are the wavenumbers in the and directions, respectively. Assuming 2 + 42 # 0, we can solve for the pressure in the frequency domain:
p
(S4)
( 2
Finally, the pressure field is recovered in the spatial domain using the inverse FFT (1FFT):
iFFT [p"].
(S5)
This spectral method offers efficient approach to calculating the pressure field without the need for labeled data or training.
BC encoding. To ensure that the solution obeys the given periodic boundary conditions and that the feature map shape remains unchanged after differentiation, we employ periodic BC padding (see Figure S2) in our architecture. This method of hard encoding padding not only guarantees that the boundary conditions are periodic, but also improves accuracy.
Poisson Block
b
Poisson Salver
M, Block
Filter Bank
(
FFT
IFFT
Filter
Bank
Pk
Gather feature terms u. Vú, v'a Vp, f, Re:
K
Pk
kz Wavenumbers in the x and y direction
Construct v(uu)
Poisson Solver
Neural Network
FFT Fast Fourier Transform
IFFT
Inverse Fast Fourier Transform
Figure $1, Components of Physics Block, (a), Poisson block. (b), Poisson solver. (c), M,NN block.
14

---

PDE-embedded Learning with Multi-time-sepppine
29
30
25
26
27
28
30
25
26
35
36
31
32
33
34
35
36
<1
32
6
1
2
3
4
6
1
2
11
-2
7
8
9
10
11
12
7
8
2.1
30
25
26
27
25
26
35
36
31
32
33
34
35
36
31
32
5
1
2
3
6
1
2
11
12
7
8
9
10
11
12
7
8
Symtric
Padding
Internal Nodes
Figure S2, Periodic BC padding.
B.3. RK4 Integration Scheme
RK4 is a widely used numerical integration method for solving ordinary differential equations (ODEs) and PDEs, commonly employed as time integration solver. It provides balance between computational efficiency and accuracy by calculating intermediate slopes at various points within each time step. The general numerical integration method for time marching from uz to Uty+1 can be written as:
u++1 = uy +
(K,
(S6)
Among them, U and u are solutions at time j + and j. RK4 is high-order integration scheme, which divides the time interval into multiple equally spaced small time steps to approximate the integral. The final update of the above state change can be written as:
ri B (uj, (j),
8t
5t +
F2 = B
u
x rı,
8t.
51
F3 B
U x r2, 1j +
(S7)
Г4 B (uj 8t r3, + 8t),
Uj+1 = uj + 8t(r1 + 2r + 2r ra).
where öt denotes the step size and r2, r3, F4 represent four intermediate variables (slopes). The global error is proportional to the step size to the fourth power, i.e., (
C. Data Details
KdV. The Korteweg-de Vries (KdV) equation is well-known nonlinear PDE used describe the movement of shallow water waves with small amplitude in channel. It models the dynamics of these waves and is particularly noted for its ability to represent solitary waves, or solutons, given by:
au
Ou
3u u
= 0,
(S8)
0t
dz
3r
15

---

PDE-embedded Learning with Multi-time-teepping
Table SI. Settings for generating datasets.
Parameters / Case
Kdv
Burgers
GS
NSE
DNS Method
Spectral
FD
FD
FV
Spatial Domain
[0, 64]
[0, lJ²
[0, 1J²
[0, 2x
Calculate Grid
256
1002
1282
20482
Training Grid
64
25
322
642
Simulation dt (s)
1.00
10-2
1.00 x 10-3
2.00 x 10-3
2.19 10-4
Warmup (S)
0
0.1
0
40
Training data group
3
3
5
Testing data group
10
10
10
10
Spatial downsample
4x
16x
16x
1024x
Temporal downsample
5X
10x
20x
128x x
where u u(z,t) represents the wave amplitude as a function of position and time th accounts for the nonlinear effects, while 8du captures the dispersive effects in the system. The equation illustrates how balance between these two da effects leads to the formation of solitary waves, solitons, that maintain their shape over long distances.
To generate the dataset, we employ the method of lines (MOL) using pseudospectral methods to compute the spatial derivatives (Brandstetter et al., 2022a). The dataset is initially generated on a grid of 256 points and then downsampled to grid of 64 points for numerical experiments. The simulation timestep is set to 1 x 10 seconds, with the total simulation duration set to 100 seconds. For training, we use sets of data, each comprising 1000 timesteps over At 5dt, along with ten additional testing sets with different ICs.
This equation models the behavior of a viscous fluid (Kumar, 2023), incorporating both nonlinear dynamics and diffusion effects. It finds extensive applications across various scientific disciplines, including fluid mechanics, materials science, applied mathematics and engineering. The equation is expressed as follows:
Burgers.
du u Vu, / E [0,T], 06
(S9)
where u {u, v) E 2 represents the fluid velocities, the viscosity coefficient set to 0.002, and A is the Laplacian operator.
As shown in Table S1, we generate the dataset using the finite difference method with a 4th-order Runge-Kutta time integration (Rao et al., 2023) and periodic boundary conditions over the spatial domain x E [0, 1]. The data is initially generated on a 1002 grid and subsequently downsampled to a 25 grid for use in numerical experiments. The simulation timestep is set to = 1 x 10-3 seconds, with total duration of = 1.4 seconds. During the training stage, we employ five trajectories with At= 10St, each consisting of 140 snapshots. In the testing stage, we use ten different trajectories, each containing 140 snapshots.
GS. The Gray-Scott (GS) reaction-diffusion model is system of PDEs that describes the interaction and diffusion of two reacting chemicals. It is known for its ability produce intricate and evolving patterns, making popular model for studying pattern formation. It is widely used in fields such chemistry, biology, and physics to simulate processes like chemical reactions and biological morphogenesis. The equation is expressed by:
Ou D,Au up + axl1 -),
(S10)
dv D.Av up (a k)D,
Øt
where and denote the concentrations of two distinct chemical species, with Dn, and D, indicating their respective diffusion coefficients. The first equation models the change in the concentration of u over time. The term D, represents based on
the diffusion of u, - uu? describes the reaction between u and v, and a(1 u) represents the replenishment the feed rate a. The second equation models the evolution of v. where D, accounts for diffusion, uv² represents the creation of e from the reaction with u, and -(a + 1s)e describes the decay of v, with as the decay rates.
16

---

PDE-embedded Learning with Multi-time-teppinng
Table S2. Settings for generating the NSE datasets.
Dataset
Grid
Spatial Domain
Re
Warmup time
dt
Innerstep
f f4
2048 64
(0, 2
1000
40
2.1914x10
32
Re 500
2048 + 642
(0, 2x)²
500
40
2.1914x10
32
Re = 800
2048 + 642
(0,2 (0,2
800
40
x.1914x10-+
32
Re 1000
2048 64
(0, (0,2
1000
40
2.1914x10-1
32
Re 1600
2048² + 642
(0, 2n)²
1600
2.1914x10-4
32
Re = 2000
2048
+ 64
(0,2z)²
2000
40
2.1914x10-4
32
Re 4000
4096
64
(0,2*)²
4000
40
100957x10-+
32
Re 1000
40962 + 64
(0,4x)
1000
40
1.0957x10-1
32
We also utilize the RK4 time integration method for dataset generation. In this case, we assign the values D.. 2.0 10-5 D. 5.0 x 10-6 0.04, and K 0.06. The dataset is generated using the finite difference method on 128? grid with periodic boundary conditions, spanning the spatial domain x E [o, 1J2. To generate different ICs, we first define a grid based on the spatiotemporal resolution and initialize the concentrations of two chemicals. By setting different random seeds and adding varied random noise, we create unique ICs. The simulation uses a timestep of = 0.5 s over total duration of = 1400 seconds. The data is then downsampled to a 322 grid, and the timestep increased to 10 seconds (åt 20df) for ground truth creation. We utilize three training trajectories, each with 180 snapshots, and ten additional testing sets with varying ICs.
NSE. The Navier-Stoker equations (NSE) are fundamental to the study fluid dynamics, governing the behavior fluid motion. this paper, we focus on two-dimensional, incompressible Kolmogorov flow with periodic boundary conditions, celecity-pressure expressed in velocity pressure form as:
du
+u
u
V2u Re
Vp+f, te [.
Ot
(S11)
u
where {u, v) ? denotes the fluid velocity vector, p represents the pressure. and Re the Reynolds number that characterizes the flow regime. The Reynolds number serves as scaling factor in the NSE, balancing the inertial forces, represented by the advection term (u V)u, with the viscous forces, captured by the Laplacian term Au. When Re is low, the flow remains predominantly laminar and smooth due to the dominance of the viscous forces. Conversely, at high Reynolds numbers, the inertial forces take precedence, leading to a more chaotic and turbulent flow behavior.
To create the dataset, we follow the approach outlined in JAX-CFD (Kochkov et al., 2021). We simulate data using Finite Volume Method (FVM) on a fine grid with a time step of dl (e.g., Re 1000, 2048 2048). This data is then downsampled coarse grid with At 128dt (e.g., Rt 1000, 64 64) serve as the ground truth. Different ICs are generated by introducing random noise into each component of the velocity field and subsequently filtering it to obtain divergence- frree field with the desired properties. For training, we utilize only five groups of labeled data with 4800 snapshots, while testing involves ten sets of trajectories. The model performance tests include trajectories with different Reynolds numbers Re 500, 800, 1600, 2000, 4000, different external forces fl cos(2y) n.. - 0.1u, f2 0 f3 = cos (44))m 0.lu, f, = sin(4y) 0.4u, and a larger computational domain x E (0, 2m The detailed dataset parameters are shown in Table $2.
D. Additional Experimental Results
D.1. Less Data and Added Noise
To evaluate our model's 's robustness against missing data and noise, we conducted experiments on the NSE using five sets of trajectories (5x xx000xxxxxx)) x4).. We tested two conditions: (1) randomly removing 20% of snapshots and (2) adding 0, % Gaussian noise during training. As shown in Table S3, the model's performance was only slightly impacted in Experiment 1, with low enror rates. In Experiment 2, HCT remained above s. These results highlight the model's strong generalization ability even under challenging conditions.
17

---

PDE-embedded Learning with Multi-time-ttepping
Table S3. Performance metries under different noise levels during training.
Training
RMSE
MAE
MNAD
HCT (s)
20% data
0.1935
0.0958
0.0113
8.1392
+ 1..1%% % noise
0.2083
0.1014
0.0123
8.0431
normal
0.1379
0.0648
0.0077
8,3566
Table S4. Performance metrics for different NN blocks.
Model
RMSE
MAE
MNAD
HCT
Model-a
NaN
NaN
NaN
0.8846
Model-b
NaN
NaN
NaN
5.2930
Model-c
0.2575
0.1507
0.0191
7.2930
Model-d
0.1564
0.0703
0.0083
8,0525
Model-e
0.2479
0.1242
0.0197
7.6346
MultiPDENet
0.1379
0.0648
0.0077
8.3566
D.2. Parametric Experiments on M,NN Block and M,NN Block
To investigate the role of the NN blocks in our model, we conducted additional comparative experiments with the following configurations:
Model-a: the M,NN Block was set to UNet and the M,NN Block to FNO;
Model-b: both the M,NN Block and the M,NN Block were set to FNO;
Model-c: both blocks were FNO with roll-out training applied the macro step (with unrolled step size of 8);
Model-d: the M,NN Block was set to DenseCNN and the M,NN Block to UNet;
Model-e: the M,NN Block was set to FNO while the M,NN Block was set to Swin Transformer (Liu et al., 2021).
All other experimental settings were kept consistent, and the results are presented in Table S4. Model-a and Model-b encountered NaN values, which can attributed to the M,NN Block requiring a model capable of robust predictions the macro step. Without such a model, multi-step roll-out training (as in Model-c) becomes necessary to enhance the model's stability in long-term predictions. When strong predictive module is employed at the macro step (e.g., UNet), the M,NN Block can be replaced with more purameter-efficien model, such DenseCNN (Model-d). Setting the M,NN Block to Swin Transformer resulted in a slight decrease in accuracy, which can be attributed to the relatively small size of our dataset, as the Swin Transformer typically excels on larger datasets.
D.3. Generalization Test on flow with Re = 4000.
To further demonstrate the superior capability of our model, we conducted an additional experiment with high Reynolds number Re 4000 (see details in Table S2) maintaining the experimental setup of Section 4.1. The result is shown in Figure 6.
D.4. Scaling against Data Size
As shown in Figure S3, our testing results demonstrate that the model exhibits a scaling law behavior, with the RMSE gradually decreasing as the amount of training data increases. Our model adheres to this scaling law as well. Moreover, even in scenarios with limited data regimes (c.g., only five trajectories), the model achieves a low level of error, highlighting its robustness and ability learn from limited data.
18

---

PDE-embedded Learning with Multi-time-ttepping
PeRCNN MultiPDENet
-1.0
Testing RMSE (log 10)
-1.5
-2.0
-2.5
2
2
2
2
2
Data volume
Figure Comparison PeRCNN and MultiPDENet across various training set sizes on the Burgers equation. z-axis represents data volume, defined as the product of trajectory timesteps and the number of trajectories (corresponding 5, 8, 16, 32, and 64 trajectories, respectively). Dotted lines denotes the linear interpolation.
E. Evaluation Metrics
We employ several metrics to assess the performance of the tested models, including Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Normalized Absolute Difference (MNAD), and High Correction Time (HCT) (Sun et al., 2023). RMSE quantifies the average magnitude of the errors between predicted and actual values, providing insight into the model's accuracy. MAE assesses the average absolute deviation between predicted and observed values, thereby indicating the scale of the errors. MNAD serves an important metric for evaluating the consistency of model outputs over time, calculating the average discrepancy across temporal data points and offering normalized measure of prediction error relative to the range of actual data. HCT gauges the model's capability for making reliable long-teet mergterm predictions. These metrics are defined as follows:
RMSE
H. H.
(S12)
MAE
1
Σ
H.
H.
(S13)
71 t=1
MNAD 1
H.
H.
(S14)
H
min
N
HCT Σ
Al (,C H.) > 0.8),
(S15)
where
COv(H, H.)
PCC(H, H.)
(S16)
TH, T..
Here, n represents the number of trajectories; H, denotes the ground truth for each trajectory; H, indicates the spatiotemporal sequence predicted by the model. The term "cov" refers to the covariance function, while "o" represents the standard deviation of the respective sequence. The Iverson bracket returns a value of when the condition (PCC(H, H.) > 0.8) is satisfied and 0 otherwise. The variable N signifies the total number of time steps.
19

---

PDE-embedded Learning with Multi-time-stepping
F. Baseline Models
Fourier Neural Operator (FNO). The FNO (Li et al., 2021) combines neural networks with Fourier transforms to effectively capture both global and local features of system dynamics. The architecture has two key components: First, applies Fourier transforms to the system state, performing convolution operations in the frequency domain capture global information. This is followed by an inverse Fourier transform which maps the data back the spatial domain. Secondly, convolutional layers are employed to extract local features directly from the system state. The outputs from both global and local components are then integrated through the application of activation functions, which ultimately yield the final prediction.
PINO. PINO (Li et al., 2024c) shares the same architecture as FNO but incorporates physics- constraints by embedding governing equations into the loss function. The loss function can be expressed as: [()) = LEg + L.Dat-- The first term of the loss function, LEg, represent an MSE loss computed using the analytical expressions of the dynamics, as defined in Equation 1. The second term Date is the data loss, which was previously discussed in Section 4.1.
PeRCNN. PeRCNN (Rao et al., 2023) integrates physics-based based principles directly into the learning framework by embedding governing equations into the neural network structure. The architecture features multiple parallel CNNs that model polynomial relationships via feature map multiplications. This incorporation of physical laws improves the model's generalization and extrapolation capabilities, enabling accurate predictions in dynamic systems governed by complex equations.
UNet. The UNet architecture (Ronneberger et al., 2015) adopts symmetric encoder-decoder structure originally designed for computer vision tasks. The encoder compresses the input by applying multiple downsampling layers to capture hierarchical features at various scales. Conversely, the decoder gradually restores the original spatial resolution using upsampling operations. Skip connections bridge the encoder and decoder, directly transferring feature maps to retain fine-grained details. This design allows UNet to merge high-resolution spatial information with deeper, abstract features, achieving accurate reconstructions of both local and global structures.
DeepONet. DeepONet (Lu et al., 2021) is designed to approximate operators and map inputs directly to outputs by leveraging neural networks. The architecture consists of two main components: the trunk network, which processes domain-specific information, and the branch network, which handles the input functions. This dual-structure approach enables the efficient learning of complex functional relationships and enhances the model's capability to capture detailed operator mappings across various applications.
Learned Interpolation (LI). The LI (Kochkov et al., 2021) employs a finite volume approach enhanced with neural networks as a replacement for conventional polynomial- based interpolation schemes in computing velocity tensor product. The network adapts to the local flow conditions by learning dynamic interpolation mechanism that can adjust to the characteristics of the flow. This enables LI to provide accurate fluid dynamies predictions even coarse grids, improving computational efficiency while maintaining prediction fidelity.
Temporal Stencil Modeling (TSM). TSM (Sun et al., 2023) addresses time-dependent partial differential equations (PDEs) in conservation form by integrating time-series modeling with learnable stencil techniques. effectively recovers information lost during downsampling, enabling enhanced predictive accuracy. TSM particularly advantageous for machine learning models dealing with coarse-resolution datasets.
G. Computational Details
G.1. Training Details
All experiments (both training and inference) in this study were conducted on single Nvidia A100 GPU (with 80GB memory) running on a server with an Intel(R) Xeon(R) Platinum 8380 CPU (2.30GHz, 64 cores). All model training efforts were performed on coarse grids (see Table 1).
MultiDENee.. The MultiPDENet architecture employs the Adam optimizer with a learning rate of 5 x 10-3, The model is trained over 1000 epochs with a batch size of 90. Detailed settings for the rollout timestep can be found in Table 1. Additionally, we use the StepLR scheduler to adjust the learning rate by a factor of 0.96 every 200 steps. The model hyperparameters are listed in Tables $5 and S6.
20

---

PDE-embedded Learning with Multi-time-stepping
Table $5. Overview of hyperparameters used in the M,NN Block.
Case
Hyperparameters
Value
Network
FNO (Li al., 2021)
Layers
6
Modes
30
Width
30
NSE
Blocks
1
Padding
periodical GELU
Inputs
{ Em (p, ú, ú, V f, Re)}
Network
FNO (Li et al., 2021)
Layers
4
Modes
12
Width
12
Burgers
Blocks
Padding
periodical
GELU
Inputs
{u
Network
FNO (Li et al., 2021)
Layers
6
Modes
12
Width
22
GS
Blocks
Padding
periodical GELU
Inputs
0
h m
Network
FNO (Li et al., 2021)
Layers
4
Modes
32
Width
64
KdV
Blocks
1
Padding
periodical ReLU
Inputs
{
FNO. The architecture of the FNO network closely follows that presented in the original study (Li et al., 2021), with the main adjustment being the adaptation of its training methodology to an autoregressive framework. The training utilizes the Adam optimizer with learning rate of 1 10-3 and batch size of 20. Training carried out for 1000 epochs, and the rollout timestep matches MultiPDENet.
UNet. We implement the modern UNet architecture (Gupta & Brandstetter, 2023) using its default settings, ensuring that the rollout timestep is consistent with that of the MultiPDENet. The StepLR scheduler is employed with a step size of 100 and gamma of 0.96. The optimizer is Adam, with a learning rate of 10-3 and a batch size of 10. The model is trained for 1000 epochs.
DeepONet. We utilize the default configuration of DeepONet (Lu et al., 2021) along with the Adam optimizer. The learning rate is established at 5 14-4, with a decay factor of 0.9 applied every 5000 steps. The model is trained using batch size of 16 over total of 20000 epochs.
We maintain the standard architecture of PeRCNN (Rao et al., 2023). The optimization process is executed with the Adam optimizer and employs a StepLR scheduler that reduces the learning rate by a factor of 0.96 every 100 steps. The initial learning rate set to 0.01, and the training conducted over 1000 epochs with a batch size of 32.
PeRCNN.
21

---

PDE-embedded Learning with Multi-time-ttepping
Table Só. Overview of hyperparameters used in the M,NN Block.
Case
Hyperparameters
Value
Network
U-Net (Gupta & Brandstetter, 2023)
Hidden size
[128, 128, 256, 512]
Blocks
NSE
Padding
periodical
Inputs
{u, Al, dx)
GELU
Network
U-Net (Gupta & Brandstetter, 2023)
Hidden size
[64, 64, 128, 256]
Blocks
Burgers
Padding
periodical {u, At, dx}
Inputs
GELU
Network
U-Net (Gupta & Brandstetter, 2023)
Hidden size
[64, 64, 128, 256]
Blocks
Padding
periodical {u. A, dx)
Inputs
GELU
Network
U-Net (Gupta & Brandstetter, 2023)
Hidden size
[64, 64, 128, 256]
Blocks
KdV
Padding
periodical {u. At, dx}
Inputs
ReLU
LI. We adopt the default network architecture and parameter settings for LI (Kochkov et al., 2021). The optimizer used Adam with 31 = 0.9 and B3 0.99. The batch size configured to 8, along with a global gradient norm clipping threshold 10-6
of 0.01. The learning rate is set
10 3, and weight decay is configured
TSM. We follow the default network architecture and parameter settings for TSM (Sun et al., 2023). The initial learning rate is set 1x10 with weight decay of x 10-4, The gradient clipping norm is configured to be 1 x 1x102 We use the Adam optimizer with 32 0.98, and the batch size is set to 8.
G.2. Computational Cost (Inference)
Taking NSE example, compared the inference time, RMSE, and HCT of MultiPDENet with the Direct Numerical Simulation (DNS) method across three cases. The comparison principle is based on the time required simulate the same trajectory length (T 8.4 s) under identical experimental conditions (a single A100 GPU). The inference time is measured from the moment the initial conditions (IC) are fed into the model until the trajectory of the same length is predicted.
The DNS settings follow JAX-CFD (Kochkov et al., 2021). According to the CFL condition, the simulated time step (dt) varies with the resolution of the DNS method, resulting in different numbers of timesteps required for calculation. DNS 2048, DNS 4096, and DNS 4096 are used as the ground truth for the three cases, respectively. Detailed comparison results presented in Table S7. The computational time for given accuracy (e.g., correlation 0.8) the NSE dataset is depicted in Table $8.
Nevertheless, we also would like to clarify that the DNS code used above was implemented in JAX, while our model was programmed in Py Torch. These two platforms have distinct efficiencies even for the same model. Typically, the codes under JAX environment runs much faster compared with PyTorch (up to 6x) (Takamoto al., 2022). We anticipate achieve much higher speedup of our model if also implemented and optimized in JAX, which is, however, out of the scope of the present study.
22

---

PDE-embedded Learning with Multt-time-stepping
Table S7. Performance comparison different methods on the NSE dataset across various cases.
Case
Method
Timestep
Infer Cost (s)
RMSE
HCT (s)
Re 1000
DNS 2048
38400
260
0
8,4
Re 1000
DNS 1024
19200
135
0.1267
8.4
Re 1000
DNS 512
9600
52
0.2674
6.5
Re 1000
DNS 64
1200
18
0.7818
2.7
Re 1000
MultiPDENet
300
26
0.1379
8.4
Re 4000
DNS 4096
76800
1400
0
8.4
Re = 4000
DNS 1024
19200
136
0.1463
6.8
Re = 4000
DNS 512
9600
52
0.2860
5.8
Re = 4000
DNS 128
2400
31
0.8658
3,6
Re 4000
MultiPDENet
300
26
0.1685
6.4
(0,4
DNS 4096
75750
1280
0
8.4
4x)
DNS 1024
19200
129
0.4638
6.6
4x)
DNS 512
9600
50
0.6166
5.2
E (0, 47
DNS 128
2400
30
0.8835
2.3
E (0,47)²
MultiPDENet
300
26
0.4577
6.7
Table Computational time given accuracy (e.g., correlation 2 0.8) on the NSE dataset.
Iterm
Re = 1000
Re = 4000
x <[00,4]]
DNS 1024
135
130
133
MultiPDENet
26
19 s
21 s
Speed up
5x
7x
6x
23


Reviews:
[{"review_id": "4p3bjp1uVn", "reviewer": "Reviewer_HqxA", "summary": "The authors proposed a framework for data-driven fluid flow simulation on uniform grid. Trying to soft-embed partial differential equations in data-driven flow simulations, the authors design a convolutional filter based on the constraints of the central difference discretization of the first and second order derivatives. To be able to achieve accurate large time step rollout predictions in the learned model, the authors generate each time step prediction by first iterate through a series of pseudo time steps through a \"physics block\", and then apply a corrector at the end to generate the corrected predictions at the next time step. Test results show that the proposed method outperforms existing benchmarks.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your constructive comments! We have addressed them thoroughly and added new figures/tables (see the **rebuttal.pdf** via https://anonymous.4open.science/r/Rebuttal-5D8B/rebuttal.pdf). These results will be added to the revised paper.\n\n\n### **Weaknesses**\n\n>**W1. Uniform grids.**\n\n**Re:** Thanks for your **thoughtful comments** on irregular geometries and complex boundary conditions (BCs). Our study focuses on fluid simulations on regular grids with periodic BCs, a common constraint shared with methods like FNO, LI, and TSM. This indicates that the limitation is not specific to MultiPDENet but represents a general challenge in the field. Similar to Geo-FNO [1], a geometry-aware netwok (e.g., a geometry encoder-decoder) could be established on top of MultiPDENet (latent structured geometry learner) to handle general geometries. Thank you for pointing out the important direction for our future work!\n\nTo demonstrate our model's generalization to complex BCs, we conducted experiments on Burgers equation with Dirichlet and Neumann BCs, while keeping other settings consistent with the original data generation setup. We tested our previously trained model for inference on 10 trajectories with complex BCs through BC encoding (**Table R1** in **rebuttal.pdf**). **Figure R1** in **rebuttal.pdf** shows predicted snapshots, confirming the model's generalizability over complex BCs.\n\n>**W2. Lacking 3D cases.**\n\n**Re:** Great remark! We followed your comment and tested our model on the 3D Gray-Scott (GS) equation. We generated 5 datasets (1 for training and 4 for testing).\n\nSnapshots of trajectory evolution from 0 to 600 s for MultiPDENet and baselines are shown in **Figure R2(a)** in **rebuttal.pdf**. **Figure R2(b)** in **rebuttal.pdf** shows MultiPDENet maintains a Pearson correlation coefficient $>$0.8 throughout the evolution. The error distribution in **Figure R2\\(c)** in **rebuttal.pdf** highlights our model’s superior performance. **Table R2** in **rebuttal.pdf** summarizes our model’s performance suggesting strong potential to solve 3D problems.\n\n### **Questions**\n\n>**Q1. Model inference speed.** \n\n**Re:** Thanks for raising this important point! It's crucial to **clarify** that the reported orders of magnitude speedups for popular models like FNO, LI and DeepONet (10$^3\\times$, 80$\\times$, 24$\\times$) are often **not** evaluated against numerical methods with comparable accuracy. McGreivy et al. [2] revealed this by reimplementing these models and comparing them under consistent precision, resulting in only 7$\\times$ speedup for FNO, while LI and DeepONet exhibited slower performance. Their evaluation effectively highlighted the prior discrepancies in speedup reporting.\n\nConsistent with [2], we demonstrate a 5–7$\\times$ speedup under the condition of sonsistent accuracy, althouh implementation in JAX may further yield extra speedup gains (e.g., $>5\\times$) [3]. We must admit that the rollout strategy used in our model inherently limits its speed, while reducing the network size provides only marginal gains. Nonetheless, MultiPDENet generalizes well across diverse initial conditions, varying $Re$, external forces, and larger domains. \n\n>**Q2. Does MultiPDENet still outperform other models at the same inference speed?**\n\n**Re:** Great question! As mentioned above, solely reducing network parameters does not substantially accelerate the inference due to the employed rollout strategy. However, leveraging a larger time step $\\delta t$ proves effective. MultiPDENet's multiscale time-stepping design allows it to circumvent the CFL condition, ensuring accuracy and stability even with increased $\\delta t$. Specifically, with $4\\delta t$ (MultiPDENet-L), we achieved an extra 4$\\times$ speedup (e.g., inference time of 6s) while maintaining a similar model accuracy, aligning with the speed of other models as shown in **Table A**. Note that MultiPDENet-L was trained based on a rollout strategy over 8 macro-steps.\n\nThese results demonstrate our model’s multiscale time stepping scheme consistently outperforms all baselines across all metrics, even at similar inference speeds. Hope this clarifies your concern.\n\n**Table A:** Comparison of MultiPDENet and baselines for NSE\n|Model|RMSE(↓)|MNAD(↓)|HCT(s↑)|Infer. cost(s↓)|\n|-|-|-|-|-|\n|UNet|0.82|0.06|3.96|7|\n|FNO|1.01|0.09|2.57|5|\n|LI|NaN|NaN|3.50|9|\n|TSM|NaN|NaN|3.75|9|\n|DeepONet|2.18|1.02|0.11|**1**|\n|MultiPDENet|**0.13**|**0.01**|**8.36**|26|\n|MultiPDENet-L|0.37|0.02|7.43|6|\n\n***Refs:***\n\n[1] Li et al. Fourier neural operator with learned deformations for pdes on general geometries. JLMR, 2023\n\n[2] McGreivy et al. Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations. Nature Machine Intelligence, 2024.\n\n[3] Takamoto et al. PDEbench: An extensive benchmark for scientific machine learning. NeurIPS, 2022.\n\n**Remark:** Please let us know if you have other questions. Looking forward to your feedback!", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes. The method itself leads to a series of limitations though, see \"other strength and weakness\" for details.", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "All experimental setups have been reviewed. Please see \"questions for authors\" for additional concerns and questions I have.", "relation_to_broader_scientific_literature": "The paper is related to the data-driven modeling of (fluid) flow. The authors have provided enough reviews of the related literature in the introduction.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "1. It should be noted that the proposed method not only has to be applied to structured grid (as is with most networks that utilize convolutions), but also has to be applied on uniform grid due to the limitation of the filters used. This is a significant limitation of this work which is likely not possible to resolve.\n2. The authors did not discuss the extension to 3D cases.\n3. The proposed framework will likely to suffer in cases with irregular geometries. The boundaries will be concerning when objects of non-cubic shapes are involved, since the authors only provide strategies to handle periodic boundaries of simple shape.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "1. Table S7 & S8, am I understanding it correctly that the model is only faster than Direct Numerical Simulation (which is well-known to be computationally expensive) by 5-7 times? \n2. If that is the case then I am deeply concerned about the inference speed of the proposed framework, since typical architectures in the domain are usually reported to be orders of magnitude faster than unsteady RANS or LES (both are cheaper than DNS due to less requirement on grid density). Please report the inference speed of the proposed network versus other benchmark cases.\n3. Continued from 2, please report the performance of different models when their inference speed is about the same, by adjusting the size of the models. It should be noted that such comparison should be performed by shrinking the size of the proposed network rather than increasing the size of the benchmark networks, as the benchmark networks should stay at around the recommended sizes in their respective papers.\n\nIf the authors can prove with sufficient evidence that the proposed framework still outperforms benchmark cases when the networks are adjusted to run at the same inference time per step, then I am happy to raise the score to 3.", "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": "u8iREStfen", "reviewer": "Reviewer_kUcx", "summary": "The paper introduces MultiPDENet, a PDE-embedded neural network with multiscale time stepping to accelerate flow simulations by integrating numerical methods with machine learning. It employs finite difference-based convolutional filters to approximate spatial derivatives on coarse grids, while a Physics Block with a 4th-order Runge-Kutta integrator preserves PDE structures for accurate predictions. To mitigate temporal error accumulation, a multiscale time integration approach is introduced, where a neural network corrects errors at a coarse time scale. Experiments on various PDEs, including the Navier-Stokes equations, demonstrate state-of-the-art accuracy and long-term stability with improved efficiency over traditional numerical methods and neural network baselines.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your constructive comments and suggestions! We have carefully addressed them, and the following responses have been incorporated into the revised paper.\n\n### **Questions** \n\n>**Q1. Essential References Not Discussed.**\n\n**Re:** We appreciate your comment. We have noted that our initial submission omitted several important black-box methods in ICLR 2022, such as MP-PDE [1], GMR-Transformer-GMUS [2], and SiT [3]. We will include these references in the related work section in the revised manuscript.\n\n>**Q2. Is the white-box modules and black-blox moduls applied on the coarse mesh for NSE?**\n\n**Re:** Great remark! To accelerate the prediction of spatiotemporal dynamics on coarse grids, we integrated the numerical scheme with neural networks through a multiscale time-stepping strategy. Consequently, the entire network, encompassing both white-box and black-box modules, is designed to operate on a unfied coarse mesh.\n\n>**Q3. Does the model need to satisfy the CFL condition?**\n\n**Re:** Insightful question! According to the CFL condition, the maximum allowable time step is given by: $\\delta t _{\\max}=\\mathrm{CFL} _{\\max} \\cdot \\frac{\\delta x}{\\left| u \\right| _{\\max}}$, where the standard choice for numerical simulations is $\\mathrm{CFL} _{\\max}=0.5$. Based on this, we obtained $\\delta t _{\\max} = 0.007$ for the coarse grid with $\\delta x = 2\\pi/64$, which indeed means that our model satifies the CFL condition in our original test setting. \n\nHowever, to further investigate the extent to which MultiPDENet can go beyond the CFL condition, we conducted experiments with increased time intervals, namely, $\\delta t_{\\max} = 0.028$ (four times the required smallest time step). The resulting model is called MultiPDENet-L. With this setup, we achieved an extra 4$\\times$ speedup (e.g., inference time of 6 sec vs. the original 26 sec) while maintaining a similar model accuracy, aligning with the speed of other baseline models as shown in **Table A**. Note that MultiPDENet-L was trained based on a rollout strategy over 8 macro-steps.\n\nIn summary, MultiPDENet does **not** need to adhere to the CFL condition, demonstrating its capability to operate effectively beyond conventional stability constraints. \n\n**Table A:** Comparison of MultiPDENet and baselines for NSE\n|Model|RMSE(↓)|MAE(↓)|MNAD(↓)|HCT(s↑)|Infer. cost(s↓)|\n|-|-|-|-|-|-|\n|UNet|0.82|0.52|0.06|3.96|7|\n|FNO|1.01|0.73|0.09|2.57|5|\n|LI|NaN|NaN|NaN|3.50|9|\n|TSM|NaN|NaN|NaN|3.75|9|\n|DeepONet|2.18|1.02|0.11|0.11|**1**|\n|MultiPDENet|**0.13**|**0.06**|**0.01**|**8.36**|26|\n|MultiPDENet-L |0.37  |0.18  | 0.02 | 7.42 |6|\n\n***Refs:***\n\n[1] Brandstetter et al. Message passing neural PDE solvers. ICLR, 2022.\n\n[2] Han et al. Predicting Physics in Mesh-reduced Space with Temporal Attention. ICLR, 2022.\n\n[3] Shao et al. SiT: Simulation transformer for particle-based physics simulation. ICLR, 2022.\n\n***Remark:*** Once again, we sincerely appreciate your constructive comments. Please feel free to let us know if you have any further questions. Looking forward to your feedback!", "claims_and_evidence": "Yes, at least for the 2D NSE, it is pretty impressive to me.", "methods_and_evaluation_criteria": "Yes, for fluid problem, the energy spectrum evaluation is good.", "theoretical_claims": "There is no such proofs for theoretical claims.", "experimental_designs_or_analyses": "Yes, for all the dynamical systems.", "relation_to_broader_scientific_literature": "I think the key contribution of the paper is demonstrating the black+white box method can work for NSE.", "essential_references_not_discussed": "There are missing some citations such as ICLR 2022 for black box methods.", "other_strengths_and_weaknesses": "The paper needs to add a specific algorithm for training and a specific algorithm for testing.", "other_comments_or_suggestions": "No, I do not have any.", "questions_for_authors": "For NSE, is the white box performed on the coarse mesh? \n\nDo we really need to satisfy the CFL condition?\n\nFor NSE, is the black box performed on the coarse mesh also?", "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": "8QF0wv6eBH", "reviewer": "Reviewer_XEPN", "summary": "The paper employs a multi-time stepping procedure in conjunction with a numerical integrator, physical operators, and neural networks to process downsampled data in both space and time. The proposed method applies time stepping on a finer temporal grid while utilizing the RK4 scheme for PDE integration. Supervision is provided through samples on coarser spatial-temporal grids, and the results demonstrate that this approach can yield accurate long-term predictions even when training data is limited.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your constructive comments. To enhance clarity, we have thoroughly proofread the manuscript and corrected all identified typographical errors. We believe these revisions will significantly improve the presentation. The updated version will be uploaded as soon as file submissions are enabled.\n\n### **Weaknesses**\n\n>**W1. The benchmarks include few recent pure data-driven models.**\n\n**Reply:** Great remark! Our proposed MultiPDENet integrates PDE priors with multiscale time-stepping scheme, enabling efficient spatiotemporal simulations (e.g., turbulence) on coarse grids with limited data. While recent pure data-driven models (e.g., modern-UNet [1]) perform well in data-rich scenarios (e.g., requiring **128 trajectories**), their performance deteriorates significantly in small data settings.\n\nWhile we have selectd widely recognized data-driven baselines (e.g., FNO, modern-UNet, DeepONet), we have also considered hybrid physics-learning models (e.g., LI, TSM, PhyFNO, PeRCNN) for comparison. Following your suggestion, we added the recent model ONO [2] as an additional baseline. Its inferior performance under data scarcity (**see Table A**) further highlights the robustness of our framework.\n\n**Table A:** Comparison of MultiPDENet and baselines for NSE.\n|**Model**|**RMSE**(↓)|**MAE**(↓)|**MNAD**(↓)|**HCT**(s↑)|**Infer. cost**(s↓)|\n|-|-|-|-|-|-|\n|UNet|0.8224|0.5209|0.0627|3.9627|7|\n|FNO|1.0100|0.7319|0.0887|2.5749|5|\n|LI|NaN|NaN|NaN|3.5000|9|\n|TSM|NaN|NaN|NaN|3.7531|9|\n|DeepONet|2.1849|1.0227|0.1074|0.1126|**1**|\n|ONO|0.6613|0.4441|0.0535|4.2356|5|\n|MultiPDENet|**0.1379**|**0.0648**|**0.0077** |**8.3566**|26|\n\n>**W2. The selection of baselines for varying experimental cases.**\n\n**Reply:** Insightful comment! Our baseline selection carefully balances pure data-driven models (e.g., FNO, UNet, DeepONet) and physics-inspired counterparts (e.g., PhyFNO, LI, TSM, PeRCNN) to ensure fairness and relevance. \n\nTo ensure consistent comparison, we added PhyFNO as a baseline for the Burgers and GS equations, with results shown in **Table B** below. So far, we have maintained a unified setting of baseline comparsion (FNO, UNet, DeepONet, PhyFNO, PeRCNN) for KdV, Burgers, and GS equations. In the case of NSE, we slightly altered the baseline models by considering two other well-recognized models specifically tailored for NSE (aka, LI and TSM). In all cases, we have demonstrated that MultiPDENet surpasses the baselines.\n\nWe believe the current baseline setup ensures a **fair** and **consistent** comparison. Thanks for your great suggestion!\n\n**Table B:** Comparison of MultiPDENet and PhyFNO for Burgers and GS.\n|**Case**|**Model**|**RMSE**(↓)|**MAE**(↓)|**MNAD**(↓)|**HCT**(s↑)|\n|-|-|-|-|-|-|\n|Burgers|PhyFNO|0.0832|0.0749|0.0599|0.5546|\n|Burgers|MultiPDENet|**0.0057**|**0.0037**|**0.0031**|**1.4000**|\n|GS|PhyFNO|0.5721|0.3579|0.3520|510|\n|GS|MultiPDENet|**0.0573**|**0.0294**|**0.0298**|**1400**|\n\n### **Suggestions**\n>**S1. Typos of \"Sovling PDE\".** \n\n**Reply:** We have thoroughly proofread the paper and corrected all typos in the revised version.\n\n\n### **Questions**\n\n>**Q1. Design motivation for the Correction Block and the effectiveness of state correction.** \n\n**Reply:** Great remark! The Correction Block is designed to mitigate the information loss introduced by resolution reduction before computing the equivalent derivatives, enabling the model to adapt to the coarse grid. Instead of explicitly recovering the high-resolution solution, the Correction Block implicitly corrects the coarse-grid state by adjusting the scaling of the equivalent derivative term. During training, this block focuses on minimizing the overall PDE residual rather than directly reconstructing fine-grid details. The state $\\hat{\\bar{\\mathbf{u}}}{_m^k}$, obtained via the Correction Block, represents a neural-corrected version of the coarse solution. However, this correction should not be interpreted as an explicit recovery of the fine-grid solution. Rather, it ensures that the equivalent derivative term computed via the Symmetric Filter is optimally adjusted to minimize the overall PDE residual. Our training objective is solely aimed at making the results more closely approximate the ground truth solution, as demonstrated in the ablation study (Model-E in **Table 3 on Page 8**). Therefore, there are no additional training objectives specifically designed for correction beyond this implicit adjustment mechanism.\n\n***Refs:***\n\n[1] Gupta et al. Towards Multi-spatiotemporal-scale Generalized PDE Modeling. TMLR, 2023.\n\n[2] Xiao, et al. Improved Operator Learning by Orthogonal Attention. ICML, 2024.\n\n***Remark:*** Once again, we sincerely appreciate your constructive comments. Looking forward to your feedback!", "claims_and_evidence": "The approach demonstrates efficiency in low training data regimes. Both quantitative and qualitative results underscore the effectiveness of the learning scheme, even when trained on a limited dataset (3–5 trajectories of 200–2000 timesteps), compared to the selected baseline methods.", "methods_and_evaluation_criteria": "The method leverages physical priors, a numerical integration scheme, and learnable components to facilitate the efficient training of a PDE solver. The evaluation setting challenges current data-driven and physics-informed baselines by operating in a low training data regime. Given the high computational cost required to obtain reliable simulations, this approach is particularly noteworthy.", "theoretical_claims": "Not applicable.", "experimental_designs_or_analyses": "The experiments are designed to showcase the learning efficiency of the method. While they sufficiently demonstrate the method's efficiency relative to the baselines in a low training data regime, it would be beneficial to include tests with larger datasets to evaluate whether the observed improvements remain significant across varying data quantities.", "relation_to_broader_scientific_literature": "The key contributions of the paper lie in the continuous effort of the community on crossing different modeling principles and ideas from physics, numerical scheme and deep learning.", "essential_references_not_discussed": "To the best of my knowledge the related work is sufficient for the context of the key contributions.", "other_strengths_and_weaknesses": "Strengths:\n- The method's architecture is described in sufficient detail.\n\nWeaknesses:\n- The paper appears to test relatively few recent purely data-driven models as benchmarks. Could the authors provide an explanation or clarify the rationale for this choice? The claim on the amount of data should exclude all purely data-driven methods.\n- Given the extent of physical prior information incorporated, it would be valuable to see the performance of the physics-informed counterparts of various architectures (e.g., FNO-PhyFNO, DeepONet-PhyDeepONet) across different cases. At a minimum, the authors should justify why certain baselines are included for some cases but not for others (e.g., PhyFNO is only evaluated for KdV).", "other_comments_or_suggestions": "Typos:\n- L241 (right colomn): \"Sovling PDE\".", "questions_for_authors": "It appears that the \"Correction Block\" is an FNO without additional training objectives specifically aimed at achieving the correction. Could the authors clarify this design choice or provide qualitative results that illustrate how the state $\\hat{u}_m^k$ is effectively corrected?", "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:

Meta-review: This paper studies a hybrid scheme to simulate Navier-Stokes models. The marching still follows traditional multi-step methods (RK4 is used in the paper), while the overall forcing uses a hybrid scheme (pressure is solved using an FFT Poisson solver, the convection part is taken care by NNs). All reviewers give positive feedbacks toward accept. However, one review was too short so I read the paper myself. Below are some comments:

- Equation (1), if $\mathbf{u}$ is a vector, I would not use $\mathbf{u}^2$ notation, and this confuses even more as later on $\mathbf{u}^k$ is used to denote the snapshots at the sub-timesteps.
- Sometimes, the writing uses M${}_i$ NN, while some other times, it is MiNN.
- Page 8: "confirming p-u decoupling necessity in the NSE." It should be "coupling". Also I suggest adding another ablation on the different placement of the Poisson block (e.g., after MiNN).

I also suggest the authors add the following references:
- Another post-processing scheme to correct surrogate models: Lippe et al. *PDE-Refiner: Achieving Accurate Long Rollouts with Neural PDE Solvers*.
- Parametrizing filters with symmetric (or skew-) contraint: Luan et al. *Gabor Convolutional Networks*.
- This paper studies the learning of particle-based time steppers (with a traditional integrator), similar to this submission: Li & Farimani: *Graph neural network-accelerated Lagrangian fluid simulation*.
- Using Poisson solver explicitly in substep methods, instead of solving a mixed formulation of divergence-free equation: Guermond, Minev, Shen. Comput. Methods Appl. Mech. Engrg. 195 (2006) 6011–6045.
Decision: Accept (poster)