Physics-informed Reduced Order Modeling of Time-dependent PDEs via Differentiable Solvers
Authors: Nima Hosseini Dashtbayaz
Hesam Salehipour
Adrian Butscher
Nigel J. W. Morris
Venue: NeurIPS 2025 poster ()
Keywords: Reduced Order Modeling
Physics-informed Machine Learning
Differentiable Physics
Primary area: machine_learning_for_sciences
TL;DR: A physics-informed reduced-order model that integrates differentiable PDE solvers into training to couple latent dynamics with governing physics, enabling generalization to unseen dynamics, long forecasting, and recovery from sparse observations.

Abstract: Reduced-order modeling (ROM) of time-dependent and parameterized differential equations aims to accelerate the simulation of complex high-dimensional systems by learning a compact latent manifold representation that captures the characteristics of the solution fields and their time-dependent dynamics. Although high-fidelity numerical solvers generate the training datasets, they have thus far been excluded from the training process, causing the learned latent dynamics to drift away from the discretized governing physics. This mismatch often limits generalization and forecasting capabilities. In this work, we propose **Ph**ysics-**i**nformed **ROM** ($\Phi$-ROM) by incorporating differentiable PDE solvers into the training procedure. Specifically, the latent space dynamics and its dependence on PDE parameters are shaped directly by the governing physics encoded in the solver, ensuring a strong correspondence between the full and reduced systems. Our model outperforms state-of-the-art data-driven ROMs and other physics-informed strategies by accurately generalizing to new dynamics arising from unseen parameters, enabling long-term forecasting beyond the training horizon, maintaining continuity in both time and space, and reducing the data cost. Furthermore, $\Phi$-ROM learns to recover and forecast the solution fields even when trained or evaluated with sparse and irregular observations of the fields, providing a flexible framework for field reconstruction and data assimilation. We demonstrate the framework’s robustness across various PDE solvers and highlight its broad applicability by providing an open-source JAX implementation that is readily extensible to other PDE systems and differentiable solvers, available at https://phi-rom.github.io.

Physics-informed Reduced Order Modeling of Time-dependent PDEs via Differentiable Solvers
Nima Hosseini Dashtbayaz' Department of Computer Science University of Western Ontario London, Ontario, Canada nhosse5Buwo.co ca
Hesam Salehipour Autodesk Research
Toronto, Ontario, Canada
hesam.. h www.salhiiiprrrattttttt.coom
Adrian Butscher Autodesk Research Toronto, Ontario, Canada
Nigel Morris Autodesk Research
Toronto, Ontario, Canada nigel higgll.oorrss@@uodess..oom
adrian. www.scchrrraattteese.coom
Abstract
Reduced-order modeling (ROM) time-dependent and parameterized differen- tial equations aims to accelerate the simulation of complex nigh-dimensional systems by learning a compact latent manifold representation that captures the characteristics of the solution fields and their time-dependent dynamics. Although high-fidelity numerical solvers generate the training datasets, they have thus far been excluded from the training process, causing the learned latent dynamics drift away from the discretized governing physics. This mismatch often limits generalization and forecasting capabilities. In this work, propose Physics- informed ROM (4-ROM) by incorporating differentiable PDE solvers into the training procedure. Specifically, the latent space dynamics and its dependence PDE parameters are shaped directly by the governing physics encoded in the solver, ensuring strong correspondence between the full and reduced systems. Our model strate-
outperforms state-of-the data-driven ROMs and other gies by accurately generalizing to new dynamics arising from unseen parameters, enabling long-term forecasting beyond the training horizon, maintaining continuity in both time and space, and reducing the data cost. Furthermore, 4-ROM learns recover and forecast the solution fields even when trained or evaluated with sparse and irregular observations of the fields, providing a flexible framework for field reconstruction and data assimilation. We demonstrate the framework's 's robustness across various PDE solvers and highlight its broad applicability by providing an open-source JAX implementation that is readily extensible to other PDE systems and differentiable solvers, available at (ttttts:: htopps://hhh //phi rom. github. io
Introduction
Many-query problems in engincering, such as design optimization, optimal control, and inverse problems, rely on exploring the solution manifold of a governing PDE in a large parameter space. Reduced order modeling (ROM) provides an appealing alternative to computationally expensive high-fidelity numerical simulations 2 Consider an autonomous dime-dependent PDE of the form
ú Nu; u(t,x): T x Ω 1
(1)
"This study was conducted during the author's internship at Autodesk. Corresponding author
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

where is an m-dimensional vector field of interest, à is its time derivative, .N is a nonlinear differential operator in space (possibly parameterized by some parameters 3), and T and & R are the temporal and spatial domains, respectively. When accompanied by appropriate boundary conditions, the goal is to evolve an initial condition u(0, ar) to any time 1 >> 0. While full-order models discretize Ω2 to grid of N coordinates and solve the resulting nigh-dimensional system of ODEs, ROMs aim to find a manifold in a low-dimensional space IR* with k < N that captures the characteristics of the solution fields and their temporal dynamics. While the computational cost of full-order models grows rapidly with the spatial dimension d and the complexity of the PDE, ROMs provide computationally efficient alternative by working in space with prescribed low-dimensional structure, allowing for real-time simulations.
Traditionally, projection-basse ROMs find linear subspace via linear dimensionality reduction techniques (e.g. Principal Component Analysis or PCA) [3. while more recently, nonlinear manifold ROMs use neural networks (e.g. auto-encoders) to learn compact latent representation a E R to better capture the nonlinear dynamies of the system 4 Among the latter, most works utilize the temporal structure of the latent space for temporal evolution of the solution fields by either finding the governing system(s) of ODEs the latent space through sparse symbolic regression [ or learning separate neural network to realize the latent dynamics through recurrent neural network 10, a transformer network M or a Neural ODE I2
In all the aforementioned works, the latent space is constructed in data-driven manner. i.c. a numerical solver first generates a large dataset of solution trajectories u(t, x) with different initial conditions and/or parameters for some time window up to time T, and then, models are trained on the generated data to learn the latent space and its evolution up to T. Nevertheless, the nonlinear manifold obtained by the data-driven training is not guaranteed to be consistent with the true physical dynamics of the system and often fails to generalize to new parameters or unseen dynamics. In particular, it is common to observe growing inaccuracies over time due to error accumulation [10 [6 where the model fails to forecast beyond the time horizon of the training dataset.
In order to enhance data-driven ROMs to be more faithful to the governing conservation laws, (i.e. to render them dthysics-informed" two categorically different methods have been proposed in the literature thus far; (i) augmenting the "PINN" loss [7 to the loss function during the offline training needed for learning the latent space 18 20
and (ii) evaluating the solutions directly in the full physical space during inference based on the exact PDEs as evaluated on a sub-sampled spatial grid. This second approach was recently proposed by Chen et al. (2 in which reduction is achieved only during inference as result of spatial sub-sampling and not directly due to the reduced space of the learned manifold. Both methods leverage Implicit Neural Representations (INRs) and automatic differentiation (AD) to formulate the PDEs used in either offline training or online inference. In all these previously published works, however, the numerical solver, which contains the true governing physics of the system, discretized, is discarded after the data generation and is left out of the training process. See Appendix for an extended discussion of the related works.
In this paper, we propose novel category of physics-informed ROMs, namely one in which the conservation laws are embedded within the training procedure through differentiable solver. We propose a Physics-informan of the system on the latent space by learning the latent dynamics directly from the numerical solver. We build the recent DINo framework introduced by Yin et al. [3] that consists of conditional INR decoder D, which maps latent coordinates a to reconstructed fields, and a latent dynamics network V (a)) = d, which models the temporal dynamics of the latent space, making the framework time and space continuous (i.e. mesh-free). DINo was shown to be effective compared to other neural PDE surrogates. Taking inspiration from hyper-reduction techniques in computational physics (22). we replace the numerical integration of & during training in DINo with direct training of á with given by the PDE solver. The feedback from the differentiable solver results in a well-structured latent space and a regularized decoder that is consistent with the true physical dynamics of the system and generalizes better unseen dynamics induced by new parameters and initial conditions. Taking advantage of the mesh-free decoder, 4-ROM also learns from irregular and sparse observations of the solution fields, and recovers the full dynamics in the sparse training setting.
Reduced- Ordder Model (6-ROM) that imposes the true physical dynamics
In summary, our contributions are as follows. (1) We introduce the novel 4-ROM framework that efficiently incorporates differentiable numerical solvers into the training process of nonlinear ROMs. (2) We show that Q-ROM generalizes better to new parameters and initial conditions, and extrapolates
2

---

beyond the training temporal horizon compared to purely data-driven training. (3) We demonstrate the effectiveness of 4-ROM compared to other alysies-informed training methods for nonlinear ROMs. (4) We further show how 4-ROM can be trained with sparse observations of the solution fields while maintaining we show that 6-ROM is robust to the underlying numerical method of the solver and its discretization scheme, and can be easily extended to other PDEs and differentiable solvers using our open-source codebase in JAX.
accuracy, providing a robust framework for data assimilation. (5) Finally,
2 Background
2.1 Problem setting
Consider the PDEs described in Eq. ① Given initial condition u(0, and parameters B, we interested in evolving the PDE in time to find the state u(t, x) for any time E T. To this end, u is discretized grid X' of N spatial locations in 12 some time steps E [0, 77, creating snapshots u' that form trajectory u {u", u!.....u u²}. A dataset U then consists of M trajectories u!:T, where each trajectory is associated with different parameter 3 and/or initial condition. In our experiments, we train the models on [0, Ter] sub-interval of each trajectory in the training dataset Utr. and evaluate the models on [0, T.r] (interpolation) and [Trr, Tiel (extrapolation) intervals of the trajectories in a test dataset U.e consisting of unseen initial conditions or parameters. Note that we do not make any assumptions regarding the spatiotemporal discretization method its consistency across the samples during training and inference.
2.2 Auto-decoder for continuous spatial reduction
To have mesh-free spatial reduction pipeline, adopt the auto-decoding scheme used in [K] where the encoder network is replaced with an inversion step (23), and the decoder is a conditional INR (2). For decoder network Da parameterized by 0 and any snapshot grid X, the goal is to learn a compact latent manifold with coordinates E IR* such that
Do(a,,
where ú is approximate reconstruction of u. We drop X' from Dg for brevity if there is ambiguity. Note that Do reconstructs a field at each coordinate x E X' separately and stacks them form the full reconstrueted field ú (see Fig. 5 in Appendix). During training, for each training snapshot of E U.,, a latent vector of is initialized with zeros and is optimized along with the decoder parameters 0 by minimizing the reconstruction loss Treai
0,1 = argmin
Σ (u).a))
0,,
0.r
E(U...))
Lyeel (ul: 0, a = El ((Do(a)). u'),
where T is the set of all training latent vectors (i.e. r= (:: hegtthee..o (e...M..)) and is measure of error.
We use normalized error for as defined in Appendix B.3
When trained, the latent coordinates
corresponding to new state are computed by inversion:
= (1, X') = argmin (u: 0,
(4)
where DJ is the pscudo-inverse of Do (replacing the traditional encoders in auto-encoders). The minimization problem in Eq. (4) is solved iteratively using a pradient-based optimizer. Note that inversion is only performed once during inference when forecasting an initial condition. See Appendix Clfor a detailed description of the inversion step during inference.
2.3
Reduced-ordered temporal evolution
The spatial reduction described in Section 2.2 forms low-dimensional latent manifold that captures the spatial information of the high-dimensional physical space. Similar to previous works, we realize the temporal structure of the latent manifold by training another neural network. This dynamics
3

---

Inversion
u X
A'
X
f
X
Do
D
To
To
V
Integrate-
a +
T (aa))dtt
Figure 1: Inference with R-ROM. The initial condition and the forecast solution at target time T. can be observed and reconstructed on arbitrary (and possibly different) grids X' and X"
network learns the temporal evolution in the latent space, which in turn can be used for time-stepping
in R2 12 13
The dynamics network, V is therefore defined such that for any latent coordinate a,
da dt
T (a) à
(5)
When trained, the inference involves evolving an initial state U (discretized on some grid X') to desired time T., For that purpose, first the inversion in Eq. 4 is solved to find the corresponding a and then the dynamics network is integrated in time by numerical ODE solver to find a-::
S/ ((
T.
t = á = D)(u'r, X)
a = ā +
(6)
With a. the predicted state of the PDE may be reconstructed ú Dp(ā "x') (possibly different) grid X". This formulation of ROM has been used before [12 19]. and as will be discussed in Section 3 build on the same framework in this paper.
Remark 2.1. Parameterized dynamics network Lee and Parish (25 introduced parameterized variant of V, where the parameter 3 of the PDE (e.g. Reynolds number, diffusivity, etc.) concatenated with before being fed into the dynamics network, allowing the dynamics network to learn the temporal evolution the latent space conditioned on the PDE parameter, i.c. (a, 3). We found that applying trainable linear transformation to 3 before concatenation with significantly improves its performance and generalization across the parameters. As such, we adopt this approach in this work whenever the PDE is parameterized (see Appendix B for the detailed architecture). Henceforth, we drop / from the notation for brevity.
2.4 Differentiable solvers
We define differentiable solver as discrete and differentiable nonlinear operator S. For time- dependent PDEs with parameters 3. evolves the state u on grid Xs such that u ju', B, Xs] We may use S to obtain dt/d which represents the temporal behaviour of the discretized PDE, analogous to its continuous counterpart /V introduced in Eq. T See Appendix B.4 for further details on deriving du/dt from S.
When is implemented within a differentiable framework such JAX or PyTorch, can be seamlessly integrated with other neural network components part of a unified architecture. This integration allows gradients computed during backpropagation to flow directly through the solver. Specifically, during training, we require the derivative of S with respect to the network parameters denoted as 0000. This is computed using AD and the chain rule as 0S/00 = (ass/OO)(a/0000 /00) where ú is the reconstructed field as defined in Eq. 2 Importantly, the term as captures key sensitivity information of the governing PDEs with respect to the state variables-informmtioon that often used in constructing adjoint equations for PDE-constrained optimization.

---

E U fr
#2, 2ah E drr
Xx
Do
Lrec
Xss
S
ú
du/ dt
Do
Difforentiale Solver
L dym
r
B
V
do dt
Figure 2: (Left) Training 6-ROM with the reconstruction loss and dynamics loss. Red and purple dashed lines show the gradient path taken by Lree and Ldyn. The reconstruction loss trains Dg and with spatial information, while the dynamics loss trains the V and regularizes the decoder and solution manifolds with physics information from S. (Right) 4-ROM can be trained with partial and irregular observations of the solution fields and recover the full fields and dynamics.
Physics-informed training with a differentiable solver
3.1 Training objective
6-ROM trained to preserve the spatial information of the physical space by minimizing the reconstruction loss Lree in Eq. 3 and to learn the temporal dynamics within the latent space by minimizing a physics-informed dynamics loss Laya.
Dynamics loss. By taking the time derivative of the decoder in Eq. 2). we get
dü
Jo(ajá di"
where (a) is the Jacobian of the decoder w.r.t. a. The time derivative of the reconstructed state ú is given by the numerical solver as explained in Appendix B.4 Ideally, the reconstructed dynamics represented by JD (a:) should follow the dynamics given by the solver. However, rather than forming the loss function (a))á dů/ct in the full physical space IR N=m we project Eq. 7 into the latent space R& with l < N by taking the pseudo-inverse to get
à Jb() dů
(8)
dt
where (a) is the pseudo- seeudo-inverse inverse of the Jacobian matrix. In other words, we essentially solve the resulting squares problem for á in Eq. 77) find â as defined in Eq. (8). We will elaborate on the computation of à via hyper-reduction in Section [33] Finally, since the latent dynamics á is realized by the dynamics network V we define the dynamics loss as:
Lddn (a, 0 = (VV
Training $-ROM. With the reconstruction loss from Eq. 3 and the new dynamics loss in Eq. we define the training objective for 4-ROM as
9
LI-ROM (u; 0, a, Φ) ALreo (a:reo(u:0 (u; + (1 A) (0, a, Φ)
(10)
minimized over all snapshots ul E Utr:
0,1, Φ = argmin
Lo-Rom(u;;0,
0.1.
(u) c(U...
(11)
5

---

In Eq. (00. À controls the regularization effect of the dynamics loss. Choice of À depends on the sensitivity of S to the noise and errors in the reconstructed fields. We chose 0.5 0.8 depending
on the solvers (see Appendix F.3 and G for the choice of λ). See Appendix B for the detailed training procedure.
Remark 3.1. Note that 6" and ú in Eq. (8) are both functions of and a. Since S is differentiable, gradients of the dynamics loss w.r.t. and @ are backpropagated through the solver discussed earlier in Section 2.4 (see also Fig. 2 As a result, Ldyn has regularizing effect on the decoder and the latent manifold, aligning them with the true dynamics of the solver. This is in contrast with [13, where the decoder and its resulting manifold are trained solely by Lyee and are not influenced by the dynamics network.
Remark 3.2. One can alternatively obtain the derivatives du/dt dt by finding the right-hand side of the exact PDEs using AD similar to PINNs SPINN-ROM, (see 19 for a closely related technique).
3.2 Training with data on irregular grids
As pointed out in Section and shown schematically in Fig. 2 we do not make any assumptions about the spatial sampling grid of the data and its consistency across snapshots. However, the solver may require its input fields to be on a specific grid, Xs. while the training data may be on a different (possibly irregular) grid, Xyx. Since 4-ROM uses INR decoder, the fields on Xs to compute the dynamics loss, while the reconstruction loss is computed on X,r- In Section [4.4, we will demonstrate that 4-ROM can be trained reliably with such irregular data sampled at a limited number of sparse locations and recovers the full fields successfully. This ability is particularly helpful for field reconstruction and data assimilation tasks.
still reconstruct
3.3 Hyper-reduction for computing the dynamics loss
QR decomposition to solve the overdetermined system of equations given in Eq. (7) in the least-squaree sense. Note that for an m-dimensional state reconstruction ú I N xm on grid Xs with N coordinates, the Jacobian Jp(a) of size N m k. whereas á R* where again < N is the reduced latent dimension. While QR decomposition enables solving such an overdetermined system of equations, the computational cost of forming the Jacobian and solving the system grows with N (and the spatial dimension d). To mitigate the computational cost, we adopt the hyper-reduction technique (see e.g. (22) in which a subset X3S of Xss with size >N is sub-sampled for each training snapshot in order to solve the least-squares problem only on the reduced X3- In this work, stochastic hyper-reduction with 7 0.1 (as opposed to other techniques such as greedy algorithms) proved to be sufficient for this purpose in all our experiments. Furthermore, we use forward-mode AD to compute the Jacobian matrix. See Appendix B.2 for detailed description of the dynamics loss computation.
In order to solve for dt defined in Eq. (8)
4 Results
4.1
Experimental settings
We evaluate @-ROM on five problems: 1D viscous Burgers' (Burgers') and 2D diffusion (Diffusion) equations, both solved by finite-difference solvers, 2D Korteweg-De Vries (KdV) equation from the spectral solver Exponax (26), 2D Navier-Stokes decaying turbulence problem (N-S) with the finite volume-based solver JAX-CFD (27], and the 2D flow over a cylinder (LBM) using the Lattice Boltzmann solver XLB (28). This diverse set of PDEs and solvers, spanning fundamentally different numerical methods, highlights the robustness and versatility of R-ROM. See Appendix D for detailed descriptions of each problem.
For each problem, models are trained on Mrr trajectories of Ter time steps and evaluated on Mi trajectories to assess forecasting and generalization performance. Evaluation is performed by fore- casting solutions for unseen initial conditions or parameters in U.e for Ter (interpolation) and Tee (extrapolation) time steps. Among the problems, Burgers' and LBM are parameterized by their source term and Reynolds number, respectively, while the others vary in their initial conditions. Please refer to Appendix [F for a detailed description of our experimental settings and Appendix H for extended results.

---

Table Diffusion and Burgers' forecasting errors (in RNMSEs) 4-ROM, DINo, and other physics-
informed ROMs.
Diffusion
Burgers'
Time
[0, Tr
[Tr Tte]
[0, [.r]
[Ter Ttel
Ø-ROM
0.080
0.034
0.021
0.028
DINo
0.089
0.051
0.021
0.060
PINN-ROM
0.081
0.042
0.088
0.348
FD-CROM
0.131
0.351
0.001
0.044
AD-CROM
0.093
0.106
0.121
0.196
JADCROM
0.456
0.856
0.090
0.212
==
T
Tte
0.0
19
0.
0.0
n0
0.2
5.
Ground Truth
50 < Ter
4.5
9-ROM PINN ROM VI-CROM
4.0
LAD CROM
3.5-
3.0
4-HOM
2.5
2.0
1.
DINo
5.0
T <
150 Te
4.5
n ROM
4.0
1 DINo i
2001
155
1.0
50
1
(a) Comparing the predicted solution to 1D Burgers' equation for an unseen initial state during (top) interpo- lation and (bottom) extrapolation periods.
(b) Forecasted fields and their corresponding absolute errors for an unseen initial condition of N-S. Plotted velocity magnitudes.
Figure 3: Forecasting Burgers' (left) and N-S (right) test samples.
4.2 Physics-informed strategies
We first compare 4-ROM with two other alternative techniques for constructing physics-informed ROMs, namely CROM [2]] and PINN-ROM (see Remark 3.2 For reference, we also provide comparisons with the data-driven method of DINo [13] In this section, we only discuss relatively simple PDEs, namely 2D Diffusion (linear) and ID Burgers' (nonlinear) equations. All models have the same decoder architecture, while RORM PINN-ROM, and DINo also have the same dynamies network architecture. We evaluate CROM in three settings: (FD-CROM) where the temporal evolution is based on the full grid using finite difference discretization, (AD-CROM) where spatial derivatives of the PDE are calculated using AD on the full grid, and (LAD-CROM) where the spatial derivatives are calculated using AD but on sub- sub-sampled grid to achieve meaningful reduction. all three settings, the training identical, and the same trained decoder used for evaluation.

---

Table 2: Forecasting enrors (RNMSE) of $-ROM and DINo for N-S and KdV problems trained with U. on grid &er. Errors are computed for initial conditions (ICs) sampled from U,, or Uze and observed on test grid Xe, and are averaged over either the interpolation window [0, Tid or the extrapolation window (Ttr, Tea]-
N-S
KdV
Time IC Set
[0, Terl
(Tr: Tiel
10, Tir
Ter, Tie Utc
Urr
Utc
Utr
Un
Utr
Ute
Utr
Nerr
X's
Nee
4-ROM
0.064
0.170
0.136
0.373
0.219
0.233
0.475
0.486
DINo
0.036
0.580
0.692
1.543
0.098
0.459
0.399
0.728
Full training
Xs,
= Ks|
4-ROM
0.066
0.170
0.139
0.378
0.219
0.233
0.475
0.487
DINo
0.041
0.580
0.692
1.580
0.099
0.459
0.399
0.728
Xrr = Xs.
Le
= 2%
Ksl
6-ROM
0.078
0.177
0.160
0.393
0.680
0.764
0.901
0.967
DINo
0.087
0.587
0.693
1.558
0.523
0.639
0.677
0.796
|
5%|
X's
|, Xie = Xss
Sparse training
4-ROM
0.077
0.192
0.148
0.397
0.238
0.248
0.485
0.499
DINo
0.046
0.584
0.717
1.450
0.165
0.543
0.498
0.851
|/tr1|
2%| X
|, Xxe
X's
6-ROM
0.092
0.189
0.169
0.394
0.268
0.280
0.546
0.567
DINo
0.183
0.594
0.726
1.517
0.759
0.902
1.222
1.396
Table shows the forecasting errors of all models on the Diffusion and Burgers' problems. For the linear Diffusion problem with smooth Gaussian blobs as the initial conditions, both PINN-ROM and @-ROM closely outperform DINo, and they all perform better than CROM, even without reducing the spatial grid. We note that AD-CROM slightly improves FD-CROM, showing its effectiveness in simple linear problems. However, for the non-linear Burgers', both AD-CROM and PINN-ROM fail grossly to capture the dynamics and exhibit higher errors than FD-CROM and DINo (see Fig. 3a). Comparing FD-CROM with o-ROM in Fig. 3a both methods accurately reproduce the dynamics of the non-linear PDE (especially the travelling shock wave) during the training window, but unlike O-ROM, FD-CROM fails to forecast the solution evolution for & >> Ter. Failure of AD-CROM and PINN-ROM in the non-linear case highlights the challenges with optimizing PINN-type losses as studied in the literature (e.g. spectral bias) 21-31
In contrast R-ROM, CROM (21 evolves the solution directly in physical space using the exact form of the governing PDE. Nevertheless, its computational savings do not stem from a reduced representation of a nonlinear manifold but rather from using limited number of sub-sampled grid points relative to the full-order discretization (i.e. JAD-CROM). Moreover, constructing PDE residuals not always feasible when sub- sampling an INR that does not output all required physical equations requires both pressure and velocity fields,
fields. For instance, solving the Navier-Stoke and thus cannot be done with an INR that only predicts velocity. Without sub-sampling, CROM offers no computational savings and still incurs significant cost during online inference, as it requires solving a least-squares optimization problem at every time step in the full space to recover the corresponding reduced latent coordinates. For these reasons and considering the PINN-ROM inaccuracies highlighted above, we only compare d-ROM with DINo in the remainder of this paper.
4.3 Forecasting performance on full grid
(N-S) We trained 4-ROM and DINo with identical networks on N-S with M,, = 256 trajectories and Tir = 50 time steps on 64 x 64 regular grid Xx Xs and evaluated them with Mr. 64 trajectories of unseen initial conditions with T.. = 200. Table 2 (top rows) reports the average forecasting errors for both test and train initial conditions (see also Fig, 3b for visual comparison).
8

---

While DINo achieves low average error of 0.02 on the training initial conditions and time horizon, 4-ROM consistently outperforms DINo in test samples and temporal extrapolation. Notably, 4-ROM performs more than four times better than DINo in the test extrapolation, demonstrating its superior generalization and temporal stability. Both 4-ROM and DINo are capable of forecasting initial conditions new (and irregular) spatial grids Xx F X,,- As shown in Table 4-ROM is able to recover the full solution field from a partial observation of the initial conditions with only 5% and 2% of the grid points.
We further trained both 4-ROM and DINo with 128 and 512 trajectories to evaluate their data efficiency. Fig. 4 shows how their forecasting errors (with initial conditions sampled from ei- ther test or training sets) change as the train- ing dataset size increases. We observe that the physics-informed training of R--RM signifi- cantly improves the generalization and prevents overfitting compared to DINo, especially with fewer training trajectories. With Mr. 512.
1.50
4-ROM
1.25
DINO
[0.7ie]
RNMSE
1.00
[Te. Ti.
0.75
Train
0.50
0.25
e-ROM closes the gap between the train and test errors.
0.00
256
(KdV) Table [ shows test and train forecast errors for the KdV equation. Both models are trained with Mr = 512 trajectories of Ter = 40 time steps and evaluated with Mie = 64 test trajectories for The 80 time steps. f-ROM again outper- forms DINo on test interpolation and extrapolation errors, while DINo fails to generalize to unseen initial conditions. Notably, however, both models perform poorly in sub-sampled inference with only 2% grid observation of the initial conditions.
Mer
Figure 4: RNMSE for N-S based on e-ROM and DINo for [0, Terl and [Tir: Tia] with increasing size of the training dataset, showing data efficiency.
64 x 64 regular grid
(LBM) For the problem of flow over a cylinder, we keep initial conditions the same and evaluate the performance of &-ROM for different PDE parameters (here 3 represents the Reynolds num- ber). We train the LBM models, with 3,, rang- ing from 100 to 200 for training and 3re ranging from 200 to 300 for parameter extrapolation. We use 40 training trajectories with Ter 100 and evaluate with Tie 125 time steps. Table 3 reeports forecasting errors for LBM with un- seen parameters 3, and 3te. DINo captures the dynamics of the parameters within the training range and achieves smaller errors than $-ROM. However, for out-of distribbutii parameters Sre. DINo fails to forecast the PDE with the unseen dynamics, while R-ROM maintains a small error. Both models are able to recover the solutions when only 2% of the training grid is used for inference.
Table 3: LBM forecast errors (RNMSE) for interpolation and extrapolation parameter (i.e.. Reynolds number) and time intervals.
Time B
[0, Tirl
Ter: Tiel BeC
Per
BiC
Ber
Xe. = Xs
-Xie
$-ROM
0.049
0.115
0.116
0.180
DINo
0.011
0.457
0.108
0.566
Xx. Xs, [/
2% X'sl
$-ROM
0.049
0.182
0.116
0.302
DINo
0.011
0.400
0.108
0.507
|
2%|
|.
Xxe
Xss
$-ROM DINo
0.065
0.188
0.150
0.303
0.369
0.412
0.420
0.439
4.4 Forecasting performance with sparse training
In the previous experiments, we trained all models with data on a regular grid, i.e. Xr, = -s. When trained with irregular and sparse observations at small number of fixed locations, 4-ROM still able to learn the dynamics and recover the full solution fields. Table 2 (bottom rows) reports the forecasting errors for sparse training of N-S and KdV with training grids that are 5% and 2% of their full grids in size. In all cases, f-ROM maintains its accuracy close to the full training scenario, while errors for DINo grow significantly, especially in the 2% case. Similarly, we evaluate sparse training for LBM with 2% grid size (see bottom rows of Table 3) Φ Q-ROM again recovers the full solution

---

fields and maintains its generalization to out-of-domain parameters, while DINo fails to capture the dynamics from the partial observations.
Conclusion
We introduced @-ROM, physics-informed ROM that embeds the governing physies in the reduced latent space by learning the dynamics directly from a differentiable PDE solver. We showed that compared other physics-informed strategies and data-driven methods for constructing ROMs, Φ- ROM generalizes better to unseen dynamics induced by new initial conditions or parameters, forecasts beyond the training temporal horizon, and recovers the full solution fields when trained tested with partial observations. The model provides accelerated forecasting within a reduced space and benefits from mesh-free construction that continuous in space and time. Furthermore, 4-ROM adheres to the convergent discretization schemes of numerical PDE solvers when learning the temporal dynamics, making it robust to various physical phenomena compared to other physics-infone -nnfore- approaches.
Despite the promising results, acknowledge that requiring access to differentiable PDE solvers lim- its &-ROM's immediate applicability, although ongoing developments in differentiable programming offer promising avenues. The inclusion of the PDE solver in the training process also increases the training cost compared to purely data-driven methods. While we believe the superior data efficiency of e-ROM justifies the increased training costs (as demonstrated in Fig. 4 tractable physics-informed training for large-scale 3D problems would require further study and optimization of the dynam- ics loss and its components. Additionally, although Ø-ROM was successfully trained with sparse data, extending the framework to data assimilation settings with limited sensor measurements is compelling area for further exploration. Finally, this study focused on PDEs with first-order time derivatives, and extending to higher-order time-dependent (e.g. wave equations) or -meeindependent PDEs left for future rescarch.
References
[1] Clarence W Rowley and Scott TM Dawson. Model reduction for flow analysis and control. Annual Review of Fluid Mechanics, 49(1)387-4477, 2017.
[2] Steven L. Brunton and J. Nathan kutz.. Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control. Cambridge University Press, 2022.
[3] Peter Benner, Serkan Gugercin, and Karen Willcox. A survey of projection-based model reduction methods for parametric dynamical systems. SIAM review, 57(4)3483-531 2015.
[4] Kookjin Lee and Kevin T Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics, 404: 108973, 2020.
[5] Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932-333) 3937, 2016.
[6] Kathleen Champion, Bethany Lusch, J Nathan Kutz, and Steven L Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 1 116((5))22445-222) 22451,, 2019.
[7] William D Fries, Xiaolong He, and Youngsoo Choi. Lasdi: Parametric latent space dynamics identification. Computer Methods in Applied Mechanics and Engineering, 399:115436 2022.
[8] Taichi Nakamura, Kai Fukami, Kazuto Hasegawa, Yusuke Nabae, and Koji Fukagata. Convolu- tional neural network and long short-term memory based reduced order surrogate for minimal turbulent channel flow. Physics of Fluids, 33(2), 2021.
[9] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash. Reduced-order modeling of advection donennten-dominated systems with recurrent neural networks and convolutional autoencoders. Physics of Fluids, 33(3), 2021.
10

---

[10] Rakesh Halder, Mohammadmehci Ataci, Hesam Salehipour, Krzysztof Fidkowski, and Kevin Maki. Reduced-order modeling of unsteady fluid flow using neural network ensembles. Physies of Fluids, 36(7), 2024.
[11] Alberto Solera- Ricoo, Carlos Sanmiguel Vila, Miguel Gómez-López, Yuning Wang, Abdulrah- man Almashjary, Scott TM Dawson, and Ricardo Vinuesa. 3-variational autoencoders and transformers for reduced-order modelling of fluid flows, Nature Communications, 15(1):1361, 2024.
[12] Dmitrii Kochkov, Alvaro Sanchez-Gonzalez, Jamie Alexander Smith, Tobias Joachim Pfaff, Peter Battaglia, and Michael Brenner. Learning latent field dynamics of pdes. In Third Workshop on Machine Learning and the Physical Sciences (NeurIPS 2020), 2020.
[13] Yuan Yin, Matthieu Kirchmeyer, Jean- Yves Franceschi, Alain Rakotomamonjy, and patrick gallinari. Continuous PDE dynamics forecasting with implicit neural representations. In The Eleventh International Conference on Learning Representations, 2023.
[14] Louis Serrano, Lise Le Boudec, Armand Kassai Koupai, Thomas Wang, Yuan Yin, Jean-Noël Vittaut, and Patrick Gallinari. Operator learning with neural fields: Tackling pdes on general geometries. Advances in Neural Information Processing Systems, 36:70581-70611 2023.
[15] Ashish S Nair, Shivam Barwey, Pinaki Pal, Jonathan MacArt, Troy Arcomano, and Romit Maulik. Understanding latent timescales in neural ordinary differential equation models for advection-dominated dynamical systems. Physica D; Nonlinear Phenomena, page 134650, 2025.
[16] Ricardo Vinuesa and Steven L. Brunton. Enhancing computational fluid dynamics with machine learning. Nature Computational Science, 2(6):358-366, 2022.
[17] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informat neural networks:
A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[18] Kookjin Lee and Kevin T Carlberg. Deep conservation: A latent-dynamics model for exact satisfaction of physical conservation laws. Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 277-285, 2021.
[19] Minji Kim, Tianshu Wen, Kookjin Lee, and Youngsoo Choi. Physics-informed reduced order model with conditional neural fields. arXiv preprint a3X1v:2412.05233, 2024.
[20] Wenqian Chen, Qian Wang, Jan S Hesthaven, and Chuhua Zhang. Physics-informed informed machine learning for reduced-order modeling of nonlinear problems. Journal of computational physics, 446:1 10666, 2021.
[21] Peter Yichen Chen, Jinxu Xiang. Dong Heon Cho, Yue Chang, G. A. Pershing, Henrique Teles Maia, Maurizio M. Chiaramonte, Kevin T. Carlberg, and Eitan Grinspun. CROM: continuous reduced-order modeling of pdes using implicit neural representations. In The Eleventh Interna- tional Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net, 2023. URL
net http://wwww......ooo.ooooooaann.
[22] Youngkyu Kim, Youngsoo Choi, David Widemann, and Tarek Zohdi. fast and accurate physies-informed neural network reduced order model with shallow masked autoencoder. Journal of Computational Physics, 451: 2022.
[23] Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove. Deepsdf: Learning continuous signed distance functions for shape representation. In Proceed- ings of the IEEF/CVF conference on computer vision and pattern recognition, pages 165-174, 2019.
[24] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Advances neural information processing systems, 33:7462-7473, 2020.

---

[25] Kookjin Lee and E. Parish. Parameterized neural ordinary differential equations: Applications to computational physics problems. Proceedings of the Royal Society A, 2020. doi: 10.1098/ rspa.2021.0162.
/oab///0000.. 14685v1
[26] Felix Koehler, Simon Niedermayr, Rudiger Westermann, and Nils Thuerey. APEBench: A benchmark for autoregressive neural emulators of PDEs. Advances in Neural Information Processing Systems (NeurIPS), 38, 2024.
[27] Dmitrii Kochkov, Jamie A. Smith, Ayya Alieva, Qing Wang, Michael P. Brenner, and Stephan Hoyer. Machine learning-accelerate accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 118(21), 2021. ISSN 0027-8424. doi: 10.1073/pmas.2101784..
BtlLhtts:://wwe pna.swwwwwwpppppp............... pnas httpa///oooeeee/////////00001010
[28)
Mohammadmehdi Ataei and Hesam Salehipour. XLB: A differentiable massively parallel lattice Boltzmann library in Python. Computer Physics Communications, 300:109187, 2024.
[29] Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. Advances in Neural Information Processing Systems,
2021.
[30] Sifan Wang. Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient patholo- gies in physics-informed neural networks. arXiv preprint arXiv:200/.04530 2020.
[31] Ghazal Farhani, Nima Hosseini Dashthayaz, Alexander Kazachek, and Boyu Wang. A simple remedy for failure modes in physics informed neural networks. Neural Networks, 183: 883:106963, 2025.
[32] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learn- ing (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint 17X1v:171.100661 2017.
[33] Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv: 1910.03193, 2019.
[34] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoor- thi, and Ren Ng. Nerf: Representing scenes as neural radiance fields for view synthesis. Communications of the ACM, 65(1):99-106, 2021.
[35] Federico Pichi, Beatriz Moya, and Jan S Hesthaven. A graph convolutional autoencoder approach to model order reduction for parametrized pdes. Journal of Computational Physics, 101:1 2024.
[36] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[37] Zong-Yi Yi Li, Nikola B. Kovachki, K. Azizzadenesheli, Burigede Liu, K. Bhattacharya, Andrew M. Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. International Conference on Learning Representations, 2020.
[38] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhat- tacharya, Andrew M. Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.
[39] Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural PDE solvers. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022.
[40] Zongyi Li, Miguel I Liu-Schiaffini, Nikola Kovachki, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning dissipative dynamics in chaotic systems. arXiv preprint arXiv:2106.0668 06898, 2021.
12

---

[41] Nima Hosseini Dashtbayaz, Ghazal Farhani, Boyu Wang. and Charles X Ling. Physics-informe neural networks: Minimizing residual loss with wide networks and effective activations. In IJCAI, pages 5853-5861, 2024.
[42] Zong-Yi Li, Hongkai Zheng, Nikola B. Kovachki, David Jin, Haoxuan Chen, Burigede Liu, K. Azizzadenesheli, and Anima Anandkumar. Physics- inforsed neural operator for learning partial differential equations. ACM/ / IMS Journal of Data Science, 2021. doi: 10. 10.1145/3648506.
URL httpp://axxii........ org/ abs/ 01111.03794v4
[43] Deniz A Bezgin, Aaron Buhendwa, and Nikolaus A Adams. Jax-fluids 2.0: Towards hpc for differentiable cfd of compressible two-phase flows. Computer Physics Communications, 308: 109433, 2025.
[44] B. Bonev, T. Kurth, Christian Hundt, Jaideep Pathak, Maximilian Baust, K. Kashinath, and Anima Anandkumar. Spherical fourier neural operators: Learning stable dynamics on the sphere. International Conference on Machine Learning, 2023. doi: 100448850/aXXX12236000883 URL
https://arxiv.. http://wwww..rraa//////aa/////// abs/ 03838v1
[45] Antonio Stanziola, Simon R. Arridge, Ben T. Cox, and Bradley E. Treeby. j-wave: An open- source differentiable wave simulator. 2022.
[46] Tianju Xue, Shuheng Liao, Zhengtao Gan, Chanwook Park, Xiaoyu Xie, Wing Kam Liu, and Jian Cao. Jax-fem: differentiable gpu-accelerated 3d finite element solver for automatic inverse design and mechanistic data science. Computer Physics Communications, page 108802, 2023.
[47] Kiwon Um, Robert Brand, Yun Raymond Fei, Philipp Holl, and Nils Thuerey. Solver-in-the- loop:
Learning from differentiable physics to interact with iterative pde-solvers. Advances in neural information processing systems, 33:61 1-6122, 2020.
[48] lya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv preprint 1.005000 2017.
[49] Przemyslaw Bogacki and Lawrence Shampine. A (2) pair of runge-kutte formulas. Applied Mathematics Letters, 2(4):321-325 1989.
[50] Patrick Kidger and Cristian Garcia. Equinox: neural networks in JAX via callable PyTrees and filtered transformations. Disserentiable Programming workshop at Neural Information Processing Systems 2021, 2021.
[51] Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
13

---

Appendix
A
Related works
Implicit Neural Representations Implicit Neural Representations (INRs) (also known as Neural Fields) are neural networks that map coordinates (such spatiotemporal coordinates) to function or signal values and provide a continuous and differentiable representation of signals, complex geometries, and physical quantities (32/44]. Conditional INRs extend INRs by conditioning the represented function on additional parameters, such as a latent code, to model parametric functions [33 or their evolution in time [3] 21. Using conditional INRs, Park et al. [2] proposed auto- encoder network is replaced with an inversion step, which requires solving a least-squares problem to find the latent code, resulting in continuous and discretization-fre (or mesh-free) encoding and decoding. In this work, we similarly use conditional INRs for autodecoding to represent the PDE solution fields.
decoder as an alternative to outo-encoders, where the petrreiizationdeppenedee
Reduced order modeling Traditionally, linear reduction methods (e.g. PCA) have been employed to identify reduced basis functions for constructing ROMs using a collection of snapshots 2 In projection-based reduced basis methods, the governing equations are then intrusively modified to arrive at truncated set of parameterized equations in the reduced space, while non-intrusive ROMs rely on interpolation or regression between reduced basis coefficients for a new parameter or time value (see (3] for more details). For highly nonlinear systems where the Kolmogorov n-width decays slowly, linear subspace methods struggle to capture the underlying patterns of the system efficiently and accurately 4 22 A prominent method for constructing nonlinear ROMs, first proposed by leverages auto-encoders to learn a nonlinear, compact, and low-dimensional
Lee and Carlberg manifold of the solution space. Almost always, this reduced space employed for inferring unseen dynamics at new parameter and time values (except for CROM 21 to be discussed below). In these cases, when the governing PDEs are parameterized or deneddeenn,n the realized latent space also needs to become parameterized or time-dependent. This has been achieved either by (i) discovering equations that govern the latent space through sparse symbolic regression [6] 7 for example,e using SINDY 5 or (ii) training a secondary network to incorporate the additional dependencies of the latent space. The latter includes a whole host of methods that vary in their choice of this secondary network to learn the dynamic evolution 1115 parametric dependency [25] B5). Notably, Yin et al.
[3] proposed DINo as a time and space continuous ROM that uses a conditional INR to reduce the spatial dimension and a Neural ODE (36] to evolve the latent coordinates.
In contrast to the aforementioned works, CROM 21 evolves the solution directly in the physical space using the exact form of the continuous governing PDE (similar to PINNs) at the inference time. Although CROM enforces the underlying equations, making it more interpretable than purely data-driven approaches, its computational savings do not stem from a reduced latent representation of a nonlinear manifold. Instead, it relies on the use of limited number of sub-sampled grid points that are reconstructed and used to calculate the temporal derivatives of the PDE using automatic differentiation. The time derivatives are then used to find the solution at the next time step which is subsequently encoded to its corresponding latent code. However, due to the errors in the gradients when taking the subsampled time derivatives, CROM is not always able to perform subsampling and instead uses the discretized (e.g. fiite-differences) PDE on the full grid to compute the time derivatives. As such, CROM essentially solves the full-order PDE in such scenarios.
Our proposed framework of d-ROM adopts similar setup to DINo 13 but benefits from a physics- informed training strategy. 4-ROM directly trains the dynamics network with time derivatives of the reconstructed fields, obtained from differentiable solver, in contrast to DINo, which relies on time- integgttrrttoon during training using neural ODEs to evolve the latent coordinates. We also leverage hyperreduction [22] for fast and efficient training of the dynamics network in R-ROM, unlike other similar approaches that also rely on the time derivatives of the reconstructed fields [9]. Moreover, o 9-ROM differs substantially from CROM in the sense that it enforces the governing physics in the latent space and evolves the solution directly in the latent space during the training as opposed to the full physical space.
and Physies-Informed Neural Networks (PINNs) 322 are the other prominent machine learning methods for solving PDEs. Neural Operators learn mappings between function spaces and can solve parameterized
Machine learning for physics Apart from ROMs, Neural Operators
PDEs at either fixed times B7 38 or autoregressively [39 40). Fourier Neural Operators perform
14

---

kernel convolution with the Fourier basis to solve PDEs on regular grids [38, while graph neural operators extend to irregular grids and free-form geometries [9 37,, although they are still limited in generalizing to unseen grids. While Neural Operators are often completely data-driven methods, given a PDE and its initial/boundary conditions, PINNs use AD and INRs for solving PDEs without requiring labelled data. While effective for solving simple PDEs, the original formulation of PINNS introduces optimization challenges with stiff equations due to its complex loss landscape [41 29 31). Furthermore, while continuous in time and space, PINNs are limited to single parameterization of PDE and are unable to extrapolate beyond their training time horizon. To alleviate these limitations, similar mhysics-informan training methods have been proposed for training neural operators 42
Differentiable solvers Although numerical solvers have long been studied and developed in computa- tional physics and engineering, the growth of machine learning and the development of differentiable programming have prompted the development of differentiable solvers that can be intertwined with machine learning pipelines. JAX-CFD [27]. XLBB [28] and JAX-Fluids 43 for computational fluid dynamics, Torch Harmonics 44 for spherical signal processing. j- j-We 45 for acoustic simulation,
and JAX-FEM 446 for finite element methods are some recent examples of such solvers implemented in JAX and Py PyTorch. Pioneered by Um et al. [47], neural networks have been used adjacent with numerical solvers to improve their accuracy and convergence by learning corrections to the numerical errors caused by discretization [27]. More recently, Apebench (26) proposed a training algorithm for auto-regressive neural PDE surrogates where a differentiable solver makes time-steps on the network predictions to improve the forecasting stability. In both methods, the interaction between the neural network and the solver is crucial for the neural network to learn the underlying physies of the system. h-ROM is more closely related to the latter, where the neural network models the physics, and the solver is only used for training. However, 4-ROM creates the correspondence between the solver and the neural network in the latent space of a ROM rather than the full physical space. To demonstrate e-ROM in this paper, we will report results based on JAX-CFD, XLB, and Apebench libraries, where each library is used for different PDE.
B Training details
B.1 Training procedure
We train 4-ROM by minimizing the reconstruction and dynamics losses together in two phases. In the first warm-up phase, I' and Dg are trained with only the reconstruction loss, and V. is trained with the dynamies loss. The warm-up epochs ensure that the reconstruction loss decreases to a small value before the full training phase begins. In the full training phase, the aggregate loss in Eq. (1O) is minimized w.r.t. all model parameters. See Algorithm [or the detailed training procedure. We choose AdamW W [48) optimizer to train 4-ROM and apply exponential learning rate decay with decay rate of 0.985 every 50 epochs. (See Appendix [F3 for complete hyperparameters.) Note that, unlike Yin et al. 13,, both loss terms are minimized together using the same optimizer in both training phases.
B.2 Dynamics loss and hyper-reduction
We used stochastic hyper-reduction to randomly sample a subset of the points the time derivatives given by the solver and solve the least-squares problem for the dynamics loss with only the selected points to reduce the training computational cost. See Algorithm 2 for a step-by-step description of the dynamics loss. Note that the random reduced coordinates X are independently sampled for each training snapshot every training epoch.
We solve the least squares problem on line in Algorithm 2 by first computing the QR decomposition of the reduced Jacobian matrix J E RRNN*xk xk.
J J+= =QQR
where Q E Nmx and R E Nmxk are orthogonal and upper triangular matrices, respectively. Replacing J with QR in the least squares problem and multiplying both sides by Q" we get:
J+a it JLQQR QR
it xQ
Rd = Q Tut..
15

---

Algorithm Training procedure of 6-ROM
1: Input: er u ir Br, der
E {0, ..., Mrr]
2: I' - -
(0. ..... Merli 1 {0,... }
3: Randomly initialize 0 and &
4: for each warm-up epoch do
5: for each (uf. a.)
)) E (Urr, T) do
6:
8..1~ + AVerLrec ( a + A) &yy ddnnm at ④)
7:
8:
0,1, optimizer
9: end for
end for
11:
for each training epoch do
12:
for each (u af) E (Utr, do
13:
0r, D , t- Vor.,boo.oo O.T,d +- optimizer(O", r ,²) 15: end for
14:
16:
end for
Since R is an trianglr matrix, we can solve the problem for á". In practice, we used NumPy's QR decomposition and SciPy's "solve_triangull"
implementations JAX for differentiability.
Moreover, we use forward-mode AD to compute the Jacobian matrix required in solving Eq. 7 This because Dg R b Xm and k < 7.m. Note that forward mode AD achieves (kk memory complexity and O(KTD) time complexity (where I'D is the computational cost of the forward pass), which is significantly cheaper than the (zyNm memory and O(yNm'I') time complexity of the standard reverse-mode AD.
Algorithm 2 Dynamics loss of ROM
1: function Layn (8, o, a)
2: ú Døla,
Reconstruct the field on the solver grid (ú R .N = (Xsl) Compute time derivative, Appendix B.4
3:
4:
X Xss
Randomly sample N coordinates (where N (X's) b Sample time derivatives at X and flatten (û+ E R2
it + ulxg
5:
6:
J v.Do(a, X3)
Jacobian via forward mode AD (J+ N Solve least-squares (à* R$) b Return dynamics loss
7: à' argmin& át- ståll
8: return ((V (a),
9:
end function
B.3 Loss measure
We use RNMSE to measure the training error in both the reconstruction and the dynamics losses, ensuring that both loss terms are of the same scale. For reconstructed field ù and label on the same grid X, the error is calculated as
Lá
((ú ().u)
(12)
where || - || is the Frobenius over the spatial axes, and the average is taken over the field channel axis (e.g. velocity directions) indexed by i.
The RNMSE for two latent time derivative du, d2 E R in the dynamics loss is calculated slightly
differently as
||án 2zll2
E(ân,
(13)
[ST(az)]1
16

---

where || 1/2 is the 2-norm and ST stops the gradients from backpropagating through the denominator. Since d2 is function of trainable model parameters 0 and ag T gradient stopping prevents the denominator and the gradients from growing large and stabilizes the training.
B.4 Computing the time derivatives
We used numerical PDE solvers to compute dú/ dt or the time derivatives of a reconstructed state ú. If du/dt dt is not output of differentiable solver after time step, we may resort to a first-order approximation of the time derivative to compute du/dt given by:
S ú
Ats
where Ais is the step size of the solver. This approach proved to be sufficient in our N-S, KdV, and LBM experiments.
Some numerical PDE solvers, on the other hand, may be sensitive to noise and errors in the recon- structed fields and return corrupted time derivatives that are suitable for training. In such cases, we let the solver apply multiple time steps on the reconstructed field to correct the noise and calculate an averaged time derivative. We employed this method in our LBM experiments by applying five consecutive solver steps.
C Inference details
Inference in @-ROM consists of three steps (similar DINo): (1) inversion, (2) integration, and (3) decoding. Algorithm provides the pseudo-code for these three steps. (1) Given an initial condition n discretized problem defined in Eq. (4) to find the corresponding latent coordinates â". The inversion is performed by minimizing the same RNMSE loss defined in Eq. (12) using AdamW optimizer (same as the training optimizer) with a zero initialization for â We perform 1000 optimization steps across all our experiments with a learning rate of 0.1. (2) For desired time T. to, the dynamics network is integrated in time to find the final latent coordinate a"' In this work, we used third-order explicit Runge-Kutee Kutta ODE solver with adaptive step size, Bogacki Bogacki-Shampiin Shampine method 49) for numerical integration. We did not observe any improvements by using other higher-order ODE solvers. (3) The predicted solution at time T. obtained by simply decoding T on any grid X', where X' and X"
(possibly irregular) grid X, 4-ROM first solves the least squares
can be different. Note that the model may not have seen the integration time T, during training, and model is continuous in time.
Algorithm 3 Inference with f-ROM
1: Input: (u X)), X", T., Do, V
2: 6 <-
3: for each inversion step do
Inversion (D)(n'o ,2))
4: a= f- Val'Dpāā", ,2), )
5: á + optimizer(a")
6: end for
7: T. + ODEESolue(a', to,T. T.)
Integrating the latent coordinates Deeodding the final state
8: T. Do(ā"-,
D Detailed description of the datasets
Table 4 summarizes the datasets used in the experiments. The datasets are generated using different PDEs and solvers, and the details of each dataset are described below. all experiments, the same solver used for data generation is also used for training 4-ROM.
17

---

Table 4: Summary of the datasets used the experiments. My, Mie. Meat are the number of training. test, and validation trajectories, Trr and Tıe are the number of training and testing time steps, and m the number of output scalar fields.
Mrr
Mie
Meal
T
Tte
A's
Diffusion
100
32
16
25
200
42 x
Burgers'
8
9
100
200
256
1
N-S
256
64
50
200
x
KdV
512
64
16
40
80
64 x 64
1
LBM
40
40
+ 40
100
125
180 x
D.1
Diffusion equation
We consider the 2D diffusion equation with constant diffusivity x: = defined on a bounded spatial domain with zero Dirichlet boundary conditions.
Du 0t
= K ( 1Pu a2u
K 2,
2x.
[r, y, 0 for (x, y) E O0.
We let Ω = [-20, 20 and sample Gaussian blobs with means Ma, E [-12, 12]. standard deviations ranging from 10, and amplitudes ranging from 0.5 to 2 sampled uniformly for the initial conditions. We use a uniform grid of 42 x 42 points for the spatial domain and a time step of 0.1 for the simulation and solve the equation using a finte-ddifference solver. We generate 100 training and 32 test trajectories with Ttr 25 and Tee 200.
D.2 Burgers' equation
We consider a ID inviscid Burgers' equation,
dw
2w 0.02e
+ 0.5
ØL
Or
parameterized by in the source term. We use the same initial profile as in Chen et al. [21] and Lee and Carlberg [4]. Training parameters (μ)) are sampled from (0.015, 0.03] to generate M,, = training trajectories with
(0.015, 0.0171, 0.199, 0.0214, 0.0236, 0.0257, 0.0279, 0.03)
We similarly generate M,.- = 9 test trajectories with f E {0.0129, 0.0161, 0.0182, 0.0204, 0.0225, 0.0246, 0.0268, 0. 0.0289, 0.0321};
which includes two extrapolating parameters. Trajectories are simulated for 200 time steps with finite-difference solver on 256 ID grid, where the first 100 time steps are used for training and the rest for extrapolation.
D.3 Navier-Stokes decaying turbulence
We model a decaying turbulence problem defined by the incompressible Navier-Stokes equations on periodic domain in 2D as
- (u + Ou Øt
1 Vp. P
2u
0,
Re
where u is the velocity field, p is the pressure, Re = UL/V L is the Reynolds number and is the density. The spatial domain a periodic box of size L, = и discretized on a uniform grid of 64 64 points. The initial conditions are generated randomly and filtered to be divergence-free and have maximum velocity magnitude of
and peak wavenumber of Density set to 1 and the kinematic viscosity is 0.01, giving Re 310.
18

---

We used JAX-CFD, a JAX-basel library that provides finite volume solvers for differentiable compu- tational fluid dynamics simulations [27]. to simulate the Wavier-Stoker equations. We chose a time step size of Als = 0.0245, and saved the velocity field every 5 steps for a total of T = 200 (1000 solver steps). The dataset contains a total of 1024 trajectories from random initial conditions, where Mr = 182,, 256, 512]. Mze 64, and 32 trajectories are used for validation. Training is performed on the first Tor = 50 time steps of the trajectories, and testing is performed on Tie 200 time steps.
D.3.1 Dataset for training with sparse irregular grid
We prepare the training data for experiments in Section 4.4 by first upscaling the 64 x 64 solver grid 128 x 128 grid by linear interpolation, and then randomly sampling 5% and 2% of the original 204 and [0.02 x 64 81 probing locations are sampled from the upscaled 128 x 128 grid for the 5% and 2% cases, respectively. We perform the linear interpolation and upscaling of the grid to ensure that the training grid is not a subset of the solver grid, which is used for inference. The selected sparse locations are fixed and the same across all training snapshots and trajectories.
grid size, i.e. [0.05 x 64*]
D.4 Korteweg- Vries equation
26 We use Exponax to simulate the Exponax library provides spectral solvers for semi-linear PDEs KdV equation defined as:
du
3V (u +1 (V @ (Vu)) = 0.03 ((V V) (V e V))
Ōt
in 2D) periodic domain. Initial conditions are randomly sampled from a random Fourier series with maximum wavenumber of 2 and normalized to have a maximum amplitude of 1. The periodic domain uniform grid of 64 x 64 points. The training dataset contains
length is 15 and is discretized Mer = 512 trajectories with T. 40 with Ats = 0.04, and the test dataset contains Mrc = 64 trajectories with Tie 80, and 16 trajectories are used for validation.
D.4.1 Dataset for training with sparse irregular grid
We sampled the sparse data for KdV similar to N-S as explained in Appendix D.3.1
D.5 Lattice Boltzmann equation
XLB is a differentiable library for solving computational fluid problems based on the Lattice Boltz- mann Method (LBM) (28) and implemented in JAX and Warp. The continuous Boltzmann equation for a distribution function f(r, E, t) is defined as
af Ōt
0f E Ox;
F, 0f ((f), p 0c.
+
(14)
where E; = dx; /tt is the particle velocity in direction i, Fi/p dE./dt is the specific body force, and the source S2(f) is the collision operator that represents the local redistribution of f due to particle collisions. XLB solves the discretized version of this equation.
We used XLB in 2D to simulate flow around circular cylinder with the D2Q9 lattice BGK collision operator, and suitable boundary conditions representing the flow inlet and outlet, and the no-slip boundaries for set of Reynolds numbers, Re UD/v where U is the mean speed of the incoming flow, D is the diameter of the cylinder, and v is the kinematic viscosity. For sufficiently large Reynolds numbers, this configuration leads to the formation of the Kármán vortex street caused by vortex shedding. We generate the data on a 180 x 82 grid with Reynolds numbers ranging from 100 to 200 and simulate for 4000 solver steps. We discard the data from the first 3000 steps of the solver when the flow is still forming. Furthermore, since LBM makes very small steps, we save the remaining simulated data only every steps. As such, each trajectory has 200 time steps, where T₁, 100 are used for training and the rest are used for testing. However, we report the time extrapolation errors in the paper for only Tie 125 as both 4-ROM and DINo fail to capture the long-horizon dynamics of LBM. (See Appendix H.3) ). The training dataset consists of No, = 40 trajectories with linearly spaced Reynolds numbers Re E 3er = [100, 200). and the test dataset contains N = 80 trajectories where 40 trajectories have Re E 3tr (parameter interpolation) and another 40 trajectories
19

---

have Re E Ste [200, 300] (parameter extrapolation). In the paper. we report the test errors for parameter interpolation and parameter extrapolation settings separately.
Finally, the data generated by XLB is represented at mesoscopic scales on D2Q9 lattice with discrete lattice velocities e; for i {1, 9). As interested in simulating the velocity, convert the simulated distributions f to velocity given by
u(x, 1) p(x,t) t)
c ((,,(),
(15)
p(r, 1) Mi,,,) 1=1
(16)
where indexes the discrete lattice directions, and is the density, See Ataci and Salehipour [28] for detailed description of the LBM equations. We use the converted velocity fields for training DINo
and use the distribution function f for training 4-ROM as explained in Appendix D.5.2
D.5.1 Dataset for training with sparse irregular grid
We sample the sparse data for LBM similar to previous cases, except that the coordinates lying inside the cylinder are added to the sampled data to ensure that the cylinder (as a boundary condition) is fully represented in the data.
D.5.2 Training and inference for LBM
Training. As noted earlier Appendix D.5 XLB models the flow using the LBM equations and hence it represents the macroscopic variables (i.e. velocity and pressure) mesoscopic scales. As such, the state variable f(z, t) is a distribution function discretized on the D2Q9 lattice, i.c. f(x, t) R To train 4-ROM, need to reconstruct the distribution functions instead of the velocity u, in order to pass the state to the solver and get the time derivatives. Thus, we train 4-ROM using the simulated distribution functions, where each snapshot has nine channels. Furthermore, modify the reconstruction loss to further constrain the reconstructions using the velocity formula in Eq. 15 For training snapshot f E 18 and the corresponding latent coordinate a, the reconstruction loss for LBM is
(17) where f = Da(a, X) and rho[-] and vel[-] return the velocity and density defined in Eq. 15 This reconstruction loss ensures that the reconstructed states meet the constraints required by the LBM solver, which would otherwise be sensitive to the errors in the reconstruction.
Lrec (f: 0, a) {(f. f) ((vel[f), vel(f)) + (rhooff) rholf)).
Inference. While we train with the distribution function, we still evaluate the model using the velocity fields (similar to DINo). As such, for an initial velocity field u, we modify the reconstruction loss used for inversion Eq. to be
[ a, 0) [(vel[Do(a)),
(18)
Note that the above reconstruction loss for inference is different from Eq. 17 which was used during training. Using the new reconstruction loss, 4-ROM still operates on velocities. Furthermore, during evaluation, we calculate R--OM'ss errors by first converting its reconstructions to velocities and then measuring the errors, so that it can be compared to the results from DINo.
E
Architecture details
o-ROM consists of two neural network components: the INR decoder De and the dynamics network L. In this paper, we used the hyper decoder of [13] and conditional SIREN (24] as decoders, depending on the problem, and an MLP with a parameterized Swish activation function (similar to [13] as the dynamics network. We briefly describe the architectures below.
Yin et al.
E.1 Hyper decoder
Hyper decoder introduced by Yin et al. IK a conditional INR decoder where the spatial coordinates are the inputs to the network and the latent code serves "amplitude modulation" by setting
20

---

Training
M
MM
21 E X
Da
ú(œ1)
Dola, 2) ~
E X
Dg
ú(#N)
Inference
(u, X)
Figure 5: Do reconstructs the solution field one coordinate at a time and stacks the solution at all coordinates for the full field.
W
T
+b
Figure 6: The parameterized dynamics network applies linear transformation the parameter before concatenating it with latent coordinates a
the biases of the network. For a coordinate æ and a latent code a, a hyper decoder with L layers is defined as
20 (xx) x,
2 (xx) = (W() + (i0) (00 @ sin(((())))))) W II,..,, 1]
x4 (x) = ( + b(L) + (a))
where (a) a and 0 (W(I....) ((.....), Wi(....)) (I....1111))) are the parameters of the network. Note that unlike Yin et al. [13] include (1) in the trainable parameters.
E.2 SIREN decoder
is INR MLP with sine activation functions that is shown to be effective in learning continuous functions and solutions PDEs 21 41. this work, similar to Chen al. (2]) we
SIREN 24
also SIREN conditioned on the latent code a, where the latent code is concatenated with the
21

---

input coordinates before being fed into the network (see Fig. 5). For a detailed description of the SIREN architecture, we refer the reader Sitzmann et al. 24
E.3 Dynamics network
We adopt the MLP network with the parameterized Swish activation function used by DINo and extend it to parameterized PDEs by concatenating a linear transformation of the PDE parameters with the input latent coordinates (see Fig. 6). For a latent coordinate a E R and a PDE parameter R", the dynamics network V is defined as:
V/ Φ (a, 3) = MLP(( @ 3,8")
= W bg.
where W3 E R and by E R along with the MLP parameters are the trainable parameters and @ is concatenation. When the PDE is not parameterized, is the only input to the network. The parameterized Swish activation function used in 13 is defined as:
Swish(x)) Sigmoid( Softplus(w)).
where w is a trainable parameter. The output layer of the network is a linear layer of size
F
Experimental settings
F1
Implementation
JAX enables fast and flexible differentiable programming by extending NumPy and SciPy and providing composable function transformations. We implemented 4-ROM in the JAX ecosystem using Equinox for building neural networks (50). Diffrax for numerical integration [51), and Optax
for gradient optimization. As result, any numerical solver implemented in JAX can be easily used for training 4-ROM, and the whole training pipeline, along with the solver step. is compiled and run together. In this work, we experimented with JAX-CFD CFD 277, XLB (28), and Exponax (26) as PDE solvers, while other --a-d solvers can be adopted as well. Furthermore, JAX transformations allow the training pipeline of 4-ROM to support multiple GPUs regardless of the implementation of the PDE solver.
As pointed out in Appendix B.2 used the JAX implementation of SciPy and NumPy linear algebra utilities for differentiable integration with the training pipeline, support of multi-GPU execution, and easy batch operations. Moreover, noted in Appendix B.2 the Jacobian matrix is computed using forward-mode AD, which is accessible through JAX's "jacfwd" function transformation. Finally, the optimizer and ESovv function calls in Algorithms I and are provided by Optax and Diffrax libraries implemented in JAX.
F.2 Compute resources for training
We trained Burgers' and Diffusion models single 48GB A6000 GPU, and KdV, N-S, and LBM models, each on four 40GB A100 GPUs. Training wall-clock times are as follows: 20 minutes for Burgers', 5 hours for Diffusion, 38 hours for N-S, 25 hours for KdV, and 18 hours for LBM.
E33 Hyperparameters
Tables 5 and 6 outline the network and training configurations, respectively. The batch sizes were chosen to be the largest that fit the available memory. In all experiments, used the same decoder configuration for 4-ROM. DINo, PINN-ROM, and CROM, as well as the same dynamics network (or Neural ODE) for R-ROM, DINo, and PINN-ROM.
We observed that a SIREN MLP decoder requires much smaller latent space compared the hyper decoder architecture from DINo. Furthermore, although SIREN is periodic function with respect to the spatial coordinates, similar to the hyper decoder, achieved better reconstruction and forecasting results in the Burgers' and LBM problems, which are not periodic. Furthermore, À = 0.5 worked well in all our experiments except for LBM, where the solver appeared to more sensitive to errors in the reconstructed fields.
22

---

Table 5: Network configuration for each problem. See Section E for detailed description of the architectures.
Latent
D
k
Arch.
Layers
Width
Arch.
Layers
Width
Diffusion
16
Hyper
3
64
MLP
64
Burgers'
SIREN
32
Param. MLP
64
N-S
100
Hyper
3
80
MLP
512
KdV
100
Hyper
80
MLP
256
LBM
4
SIREN
6
64
Param. MLP
64
Table Training hyperparameters for each problem.
Learning rate
Decay rate
Epochs
Warm-up epochs
Y
Diffusion
0.005
0.985
24000
400
0.5
0.1
Burgers'
0.005
0.985
18000
1000
0.5
0.1
N-S
0.005
0.985
24000
100
0.5
0.1
Kdv
0.005
0.985
12000
200
0.5
0.1
LBM
0.001
0.985
20000
2000
0.8
0.1
E4 Error measurement
We reported the root normalized errors in the main text using RNMSE measure similar to the one defined in Eq. 112 For label trajectories u TT:Tas (j 0, ...,, M) in dataset U and the forecasts
u Ti:T's the RNMSE is calculated as:
T2
M
m
RNMSE := Mm - Ti)
(19)
t=
where || |Ix is the Frobenius norm over the spatial axes. As also suggested in [26]. this error measure ensures that the model is evaluated for capturing the change the scale of the solutions over time.
E.5 Training and inference settings
Full grid training, All experiments in Sections 4.2 and 4.3 are conducted by training the models with data on the full solver grids, i.c. X,, Xs. The test grids Xre in Table are all equal to Xs (except for JAD-CROM). Test grids Table are either the full solver grid or randomly sampled subsets of the solver grid with 5% and 2% of its size. In the latter sub-sampled inference case, only the initial conditions f on the sub-sampled grids are employed for the inversion, while the models still reconstruct the predicted solutions X's and are evaluated for recovering the full field. In summary, the evaluation steps for 4-ROM and DINo are as follows.
1. Given initial condition t grid Xte, find 6 by inversion.
2. For the time interval [0, Tiel (containing training and test time steps), integrate V starting with 6 to find all the corresponding "Q...., T
T..
3. Reconstruct the forecast latent coordinates at all time steps on grid Xs find i a
4. Measure the average temporal interpolation ((u ((u....., u Ter
and extrapolation (T aTte)) errors separately (as explained in Appendix E.4b using the labels 1 Ter.
Training on sparse irregular grids. Experiments Section 4.4 are conducted by training models on sub-sampled training grids. For each dataset, fixed training grid X,, is randomly sampled from Xss with 5% and 2% of the original grid size (see Appendix D). This grid is fixed and consistent among all the trajectories and snapshots in dataset. 4-ROM and DINo then trained with the
23

---

data only observed at the sub-sampled training locations. During testing, the models are evaluated for recovering the full grid solution fields similar to the full grid training case.
G On the regularizing effect of X
The hyperparameter À balances the two loss terms in the training objective. Note that while the latent manifold is constructed by Lree. it is regularized by Ldyn in order to enforce the PDE. Consequently, smaller increases the regularization effect of Ldyn on the decoder and the latent manifold. resulting in higher reconstruction loss. At the same time, a very large choice of À (> 0.9) may result in disparity between the actual latent dynamics and what the dynamics network // learns as the latent dynamics. Thus, we suggested the 0.5 to 0.8 range for the choice of A. A larger À may improve the forecast error very early time-steps in exchange for higher error accumulation, while smaller values may deteriorate the reconstructions and result in wrong time derivatives from the solver during the training. In addition, note that as we use normalized errors for both loss functions, the losses would often be in the same order of magnitude when À is close to 0.5 and decrease with similar rates. Below, we further explain our choice of in the case of Burgers' and LBM experiments.
G.1 Burgers' ablation study
To better demonstrate the effect of X, we trained Burgers' with À = 0.1, 0.5, 0.9 and report the results in Table 7 As discussed above, increasing decreases the regularization effect and thus results in smaller reconstruction loss, but fails to minimize the dynamics loss. Conversely, small À increases the reconstruction error in a way that the dynamics learned from the solver may be wrong or noisy.
Table 7: Loss and error values for Burgers' with different choices of X.
Tyrec
Ldyn
Error Ter
Error The
0.9
0.002
0.081
0.137
0.150
0.5
0.019
0.013
0.021
0.028
0.1
0.176
0.011
0,459
0.613
G.2 LBM case
In our experiments, À = 0.5 worked well for all problems, except for LBM, which has a more complex geometry because of the discontinuity the cylinder. Due to the imperfect reconstruction of the cylinder, we observed noise and unexpected artifacts in the time derivatives given by the PDE solver. To alleviate this issue, we increased achieve smaller reconstruction loss, which helped with stable training of the model and improving its generalization.
H Extended results
H.1 Mean squared errors
In the main text, we reported the root normalized errors for our experiments. In cases where the scale of the solutions decays over time (Diffusion, N-S, and KdV), non-normalized error measures such as mean squared error (MSE) may indicate a decreasing error over time, while normalized errors penalize the model for not adapting to the changing scale of the solutions and indicate a growing error. Nevertheless, many other works report results MSE. For completeness, the mean squared emoor for all experiments and both training and test sets are reported here in Tables g and
H.2 Error over time
Figure 7 shows the train and test RNMSEs at every interpolation and extrapolation time step. We note that both 4-ROM and DINo fail in longer forecasting beyond Tis for LBM.
24

---

Table 8: Similar Table experiments.
in the main text but indicating mean squared errors (MSE) for all
Diffusion
Burgers'
Time
[0, Tirl
[Trr Ttel
[0. Tirl
Ter, Tiel
Dataset
cr
U be
U fr
Ute
Utr
Urr
Ute
4-ROM
6.0c-4
8.2e-4
2.9e-5
3.8e-5
2.5c-3
2.6c-3
9.6e-3
1.1e-2
DINo
3. 3.Ge-4
8.8e-4
9.4c-5
1.4e44
2.6c-3
2.7e-3
6.4e-2
1.0c-1
PINN-ROM
5.3c-4
8.2c-4
6.0c-5
8.0e-5
7.3e-2
7.5c-2
1.4++0
1.3e+0
FD-CROM
1.0e-3
1.7e-3
5.1e-3
6.3e-3
9.8e-6
1.2e-5
5.7e-2
8.1e-2
AD-CROM
4.0e-4
9.3e-4
4.5e-4
6.0e-4
1.1e-1
1.le-2
4.2e-1
4.3e-2
JAD-CROM
2.2e-2
3.2e-2
3.2e-2
4.1e-2
8.8e-2
7.2e-2
5.6e-1
5.8e-1
H.3 LBM error over parameters
Figure shows the LBM errors for Reynolds numbers from the training range Brr and extrapolation range Ste- While DINo outperforms 4-ROM the parameter interpolation regime, it fails to extrapolate beyond the training parameters, highlighting the impact of 4-ROM's physics-informed training regularization in generalizing to unseen dynamics. See Fig. 13 for parameter interpolation and extrapolation forecasting examples.
H.4 Seed statistics
We trained 6-ROM and DINo with three different random initializations (same initializations for the two models) for N-S. Average and standard deviation of RNMSEs over the three models are reported in Table [0 Both models achieve a small standard deviation and are stable w.r.t. the initialization.
H.S
Predicting solution dynamics
Figs.
demonstrate predicted fields at multiple times based on unseen initial conditions parameters and provide comparisons with baselines for all PDEs investigated in this study.
3
25

---

Table 9: Similar Table 2 in the main text but indicating mean squared errors (MSE) for all
experiments based on DINo and 6-ROM.
1.4e4
3.4e-3
1.8e-3
3.2c-3
1.8e-3
3.0e-3
Ste
Tae
[T
5.9c-5
2.6e-3
4.8e-4
2.4e-3
4.8e-4
2.1e-3
Ppr
LBM
5.0c-5
5.3e-5
5.0e-5
5.3e-5
5.0e-5
5.3e-5
[0. TorT
8.5c-6
3.6e-7
8.4e-6
3.7e-7
8.5e-6
3.7e-7
Ber
1.7c-2
3.6e-2
1.7e-2
3.Ge-2
5.7e-2
4.2e-2
U1r
[Tr. Tea
= 55%/'s
= -%
Urr
N's = N
1.6c-2
1.lc-2
1.7e-2
1.le-2
5.5e-2
3.1e-2
KdV
Arr = Vs. 1/
Xr = As. 1/4
7.0e-3
2.6e-2
7.0e-3
2.Ge-2
5.5e-2
4.6e-2
Uie
Nr =
[0,
Upr
6.3e-3
1.le-3
6.3e-3
1.2c-3
5.3e-2
3.2e-2
Ute
2.0c-3
3.0e-2
2.0e-3
3.3e-2
2.2e-3
3.2e-2
(r. Tiel
Uir
2.2c-1
5.1c-3
2.3e-1
5.
2.9e-4
5.1e-3
N-S
1.7e-3
1.9c-2
1.7e-3
1.9e-2
1.8e-3
2.0e-2
[0, Tp.
2.4c-4
1.le-4
2.5e-1
1.4e-4
3.5e-4
6.3e-4
Urr
Time
Dataset
4 N-ROM
DINo
4 Q-ROM
DINo
i F-ROMM
DINo
26

---

9-ROM
DINo
PINN-ROSH
FISCROAT
AD-CROM
RNMSD
LAD CDOM
RMMSE
BRROST
Than
DIND
Rol
10-2
PINN-ROM
ED&CROM
AD-CROM
LAD CROM
10
Teet
25
50
100
125
150
175
200
(a) Diffusion
(b) Burgers'
10P
RONSSE
RNAISE
10
O-ROM
O-ROM
DINo
DINo
10
Thain
Test
Tred
30
75
100
125
150
10
(c) N-S
(d) KdV
20:1
8-ROM
DINo
Bc-
50
75
1.00
125
150
175
200
(e) LBM
Figure 7: Error accumulation over time for each dataset. Vertical dashed lines indicate Ter The errors for LBM are plotted beyond the test temporal horizon (The 125) reported in the paper.
0.8
3-ROM
Pry
3t0
0.6
RNMSE
0.4
0.2-
0.0
100
160
200
250
300
Rt
Figure 8: LBM forecasting ((0 - Terl) errors for Reynolds numbers sampled from the training and test intervals.
27

---

Table 10: Seed statistics of R-ROM and DINo trained for N-S with three different random initializa- tions (same initializations for both models).
Time
[0, Terl
TTr, Tie.
IC Set
Uer
U.c
U.r
Ute
Average RNMSE
h-ROM
0.066
0.179
0.134
0.370
DINo
0.033
0.578
0.723
1.554
Standard deviation of RNMSE
4-ROM
5.6e-3
7.3e-3
1.8e-3
6.6e-3
DINo
1.8c-3
7.0e-3
2.5c-2
1.4e-2
28

---

Coadition
Condition
Tip
0.0
0.5
0.0
0.5
1.0
0.0
0.5
0.0
0.2
Ground Trota
Ground Troth
DINo
DINo
AD-CROM
AD-CROM
4-ROM
O-ROX
0.00
0.10
3.02
0.04
0.0220 0500 078
0.01
0.02
Tu 0f DINo
DINo
0.05
0.10
0.05
0.10 6.15
0.05
0.10
0.09
0.04
It of AD-CROM
AD-CROM
0.025 0.050 0.075
0.005 0.010
B.Da
0.305 0.010
1N 0) N-ROM
9-ROOM
Figure 9: Predicted 2D fields for the Diffusion problem based on test initial condition. Comparing results obtained by DINo and f-ROM which are both trained with Nor = 100 and X,, = X'.
29

---

(a) Ttr
(b) Tte
Garand Touls
4-ROM
AD-UROM
LAD-CROM
DIN
(a) Predicted solutions of the Burgers' problem at (a) L = Tor 100 and (b) The 200
u(x, 1)
u(x L)
1.
150
150
100
T+
50-
30
ROM
100
200
7/
Ti
2.5
100
1.5
1.0
0.5
FD-CROM
FD-CROM
Lm
100
8
Ti
50
AD>CROM
AD-CROM
T2
km
150
Th
30
JAD CROM
AD CROM
150
50
PINN-ROM
PINN-ROM
-8
T
50
DINo
DINo
100
200
100
200
(b) Predicted solutions
Figure 10: Predicted dynamics of the Burgers" problem for an unseen test parameter μ 0.0289.
30

---

Condition
The
MAR
Man
U.4
DML
DMN
0.02
Ie- 4-10OM
Iniia Coniiiion
25
The
50
100
The 200
0.0
Goand Trash
DRo
02 D.
01 02 05
611 69
0.1 0.
-)--)--- Difl
0.000 0.050
6-19 0.04
0.01 0.00
2 1000
Figure 11: Predicted 2D fields for the N-S problem based on a test initial condition. Comparing results obtained by DINo and 4-ROM which both trained with N₁, 256 and X,, Xs. Plotted is the velocity magnitude and the corresponding absolute errors.
31

---

Initial Condision
r= Tir 10
Toe 80
Initial Condision
Thr 10
Tie 80
Ground Trath
Grouad TFaat
0.25
0.50
0.25
0.50
0.2
G)
0.50
4)
Figure 12: Predicted 2D fields for the KdV problem based test initial condition. Comparing
results obtained by DINo and 4-ROM which are both trained with N = 512 2 and X,, = Xs-
32

---

Condition
Txe
Condition
Ground Touth
Ground
DINo
DINO
2
9-RO3
3300
0.002
0.08%
0.005
0.050
0.1
0.2
0.1
0.1
le fl DINo
4-- NINO
0,005 0LDLD
0.03
0.00
0.025
0.0L 5.02 0.23
0.02 0.04
6,020 0.000
10 10--8 81 B-ROM
to [o- B-M
(a) 125 C Ber
(b)
Btc
Figure 13: LBM predicted dynamics for two Reynolds numbers Re 125 (left panel) and Re 250 (right panel) respectively from the interpolation and extrapolation parameter intervals. Comparing the magnitude of the velocity and the corresponding absolute emror obtained by DINo and & R-ROM.
33

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We proposed a physics- physics-informed reduced-order model that improves the generalization and temporal stability of the previous data-driven training approaches. In Section 3.1, introduced the dynamics loss, where numerical PDE solver is used for training the latent dynamics of the O-ROM, making it physics-infommed In Section 4 we showed that @-ROM improves the temporal stability for Diffusion, Burgers', Navier-Stokes, and KdV equations (Tables T 2 and improves the generalization to out-of-distribution parameters for Lattice Boltzmann equations (Table 3}. We also claimed that 4-ROM is able to recover the full solution fields and dynamics when trained with sparse observations, which is supported by our experiments in Section 4.4
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize to other settings. It is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We discuss the limitations of the proposed method in Section 5 Throughout the paper, we elaborate on our experimental settings and assumptions (as in Section 2.1 and address potential limitations, such as training costs discussed in Section 3.3 and Appendix
B.2
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-ted system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they seale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
34

---

While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [NA]
Justification: The paper does not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: All hyperparameters, including model and training configurations, are listed in Appendix [F3 Furthermore, the code for generating the datasets and training all models (including baselines) is submitted as supplemental material and will be published.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
35

---

(c) If the contribution is a new model (e.g.. a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The code for training the proposed model and all the baselines is submitted as supplemental material and will be published. For each experiment, bash script included to train and reproduce the experiment. The code for generating all the required datasets is also provided.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.ca
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
liiiiss.cc/pulii//uuu
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We specified all data generation and splitting details Appendix D Dnd Table Furthermore, training details including pseudocodes, loss measures, and optimization details described in Appendix B along with network and training hyperparameters for each dataset in Appendix F3
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
36

---

7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We report seed statistics for the Navier-Stokes equation in Appendix HH.4 including the average and the standard deviation errors over the random initialization of the proposed method and the baseline.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula,
call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report I-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
show in tables
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer:[Yes
Justification: have indicated the compute resources used for this research in Section Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: The research conducted in the paper conforms with NeurIPS code of conduct in every aspect.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
37

---

If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This paper focuses on improving reduced-order physical simulations for real-time applications in science and engineering. While, like any computational tool, the proposed method could misused in unintended contexts, does not directly enable or pose specific risks of negative societal impact. Similarly, the potential positive societal impact lies in accelerating scientific discovery and engineering design by making fast simulations more reliable, accessible and efficient.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative socictal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The data used this work involves physics simulations, as engineering problems, and is generated by open-source computational tools. As such, this works does not involve data or models that pose a high risk of misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
38

---

We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We use open-source libraries for generating all the data and implementing the models, and we cite all the used libraries and their accompanying papers.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, papprssii thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors are encouraged to reach out the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We have documented model implementation and data generation through the pseudo-codes in Appendices and and the data descriptions Appendix D
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This research did not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
39

---

Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This research did not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing. editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: This research does not involve LLMs.
Guidelines:
The answer NA means that the method development in this research does involve L.MSs as any important, original, or non-standard components.
Please refer to our LLM policy (https: /ttp:://////pp...eee cttpeiissscco///oooo /ot/oeeennnee//22///// for what should or should not be described.
40


Reviews:
[{"review_id": "IT5cYYK437", "reviewer": "Reviewer_axmu", "summary": "This work proposes a physics-informed reduced order modeling to simulate fundamental PDEs, including Diffusion, Burgers, Navier-Stokes, KdV, and LBM equations. The approach integrates differentiable solvers into the training process to constrain the latent dynamics and ensure consistency with governing physics. The paper validates the efficacy of the proposed framework for interpolation within the training time horizon and extrapolation beyond training time on unseen initial conditions and parameters. The method is also evaluated under sparse observation settings to demonstrate its robustness.", "questions": "1. This work is dependent on how latent dynamics evolve over time. How could this framework be adapted or modified for the PDEs not reliant on time? Such as elliptic PDEs\n\n2. Operator learning methods have recently shown success for extrapolation tasks. There are no comparisons with recent operator learning methods, such as convolutional neural operators. Also, a comparison with the method PDERefiner would enable to assess the potential of the proposed method in an extrapolation scenario. Could the authors comment on this and explain how Φ-ROM compares to these methods?\n\n3. The discussion mainly focuses on accuracy for interpolation and extrapolation tasks. Since this framework evaluates trajectories over time, how does training time scale with the increasing complexity of PDEs, higher spatial dimensions, or higher-order systems?\n\n4. The training time horizon $T_tr$ and testing time horizon $T_{ex}$ differ for each test case. What are the criteria for selecting these time intervals, and how sensitive are the results to this choice?", "limitations": "Yes, discussed", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\n1. The method presentation is clear, and the training procedure and architecture description are well structured.\n\n2. The authors provided the code, which supports reproducibility and further experimentation.\n\n3. Several PDE test cases with varying levels of complexity are considered, and several comparisons with existing approaches are provided.\n\nWeaknesses:\n\n1. Comparisons with methods focused on extrapolation tasks are missing. Several neural PDE solvers already address the challenge of extrapolation, such as PDE Refiner. A comparison or discussion with respect to those methods is missing in the paper.\n\n2. Training time and computational cost of the proposed method are not discussed or benchmarked, which is essential for assessing practical feasibility.", "was_revised": true, "final_justification": "During the rebuttal phase the authors provided additional experiments comparing the proposed method with FNO variants, and the clarification on scaling, memory, and training time, which strengthen the paper. Hence, I increased the rating to 4.", "rebuttal": "We would like to thank the reviewer for the constructive feedback. Below, we address the questions (Q) and the weaknesses (W) raised in the review.\n## Q1:\nWhile we do believe that our fundamental idea of incorporating differentiable PDE solvers in a neural network for constructing reduced order models can be extended to PDEs with no time-dependence, our current paper is mainly focused on time-*dependent* problems. $\\Phi$-ROM connects the time variance of the solution in the physical space ($d\\mathbf u/dt$) with its counterpart in the latent space ($d\\alpha/dt$) and aims to learn the dynamics of the latent coordinates in time in the reduced latent space. This physics-informed training strategy is hence enabled by a secondary dynamics network, where the correspondence between the solver and the latent manifold is created. However, in a neural ROM for time-*indepdent* PDEs, we imagine that the dynamics network would be replaced by a mapping between two latent spaces, e.g. one latent space corresponding to boundary conditions and another corresponding to the solution space (see CORAL [1] for a closely related work). As the notion of \"temporal dynamics\" is not relevant in time-independent problems, the correspondence between the solver (i.e. physics) and the latent mapping function must be created in another way. As a result, directly extending $\\Phi$-ROM (in its current form) to those scenarios is not trivial to us at this point. Investigating the best strategy to make a correspondence between the differentiable solver and the latent mapping would be a very interesting direction for future research.\n\nIn light of the valid point raised by the reviewer, we will address this in the future work section of our paper.\n\n[1] Serrano, Louis, et al. \"Operator learning with neural fields: Tackling pdes on general geometries.\" Advances in Neural Information Processing Systems 36 (2023): 70581-70611.\n\n## Q2 and W1:\nWe greatly appreciate the reviewer’s feedback. We were hesitant to provide direct comparisons with non-ROM neural surrogates during the preparation of our submission, as we believed such a comparison would not be an apples-to-apples comparison. While it is true that both neural ROMs and other neural surrogates (e.g., neural operators) aim to accelerate the solution of PDEs, they do so in fundamentally different ways. The former evolves the dynamics in a reduced subspace and coordinate system, whereas the latter operates in the full space. This distinction becomes especially significant at scale, particularly for high-resolution and 3D problems. \n\nThat being said, to address the comments raised by the reviewer, we have conducted new experiments to compare $\\Phi$-ROM with two non-ROM neural surrogates by just focusing on the accuracy of the unrolling predictions.\nWe trained autoregressive FNOs with one-step training [1] and pushforward training (FNO-PF) [2] for the Burger's and LBM problems reported in the paper. We used the Apebench library [3] for training the FNO networks. As shown below for both problems, $\\Phi$-ROM outperforms FNO in temporal interpolation and extrapolation tasks. \n\n|| $T_{tr}$ |$T_{te}$|\n|-:|:-:|:-:|\n| $\\Phi$-ROM |0.021|0.028|\n|FNO-1 |0.037  |0.131|\n|FNO-PF |0.049|0.089|\n\n**Table 1.** Burger's unrolling errors for test samples. $\\Phi$-ROM errors are the same as those reported in Table 1 of the paper. We used a 3-layer FNO with 32 hidden channels, 12 Fourier modes, and GELU activation. \n\n|| $T_{tr}$ | $T_{te}$ |\n|-:|:-:|:-:|\n| $\\Phi$-ROM |0.065|0.150|\n|FNO-1 |0.145|4.867|\n|FNO-PF |0.234|3.276|\n\n**Table 2.** LBM unrolling errors for test samples (parameter interpolation). $\\Phi$-ROM errors are the same as those reported in Table 3 of the paper. We used a 4-layer FNO with 64 hidden channels, 16 Fourier modes, and GELU activation.\n\n[1] Li, Zongyi, et al. \"Fourier neural operator for parametric partial differential equations.\" arXiv preprint arXiv:2010.08895 (2020).\n\n[2] Brandstetter, Johannes, Daniel Worrall, and Max Welling. \"Message passing neural PDE solvers.\" arXiv preprint arXiv:2202.03376 (2022).\n\n[3] Koehler, Felix, Simon Niedermayr, and Nils Thuerey. \"APEBench: A benchmark for autoregressive neural emulators of PDEs.\" Advances in Neural Information Processing Systems 37 (2024): 120252-120310.\n\n## Q3 and W2:\nWe understand the scaling concerns regarding the 3D and high-resolution problems. To provide a fuller picture of computational implications of such problems, we profiled $\\Phi$-ROM's training costs for a high-resolution $512\\times512$ 2D KdV problem ($|\\mathcal X_\\mathcal S| = 2^{18}$). \n\n| PDE; Res.; $\\gamma$ | Total $\\nabla_\\theta L_{\\Phi-ROM}$ | Only $\\nabla_\\theta L_{rec}$ | Only $\\nabla_\\theta L_{dyn}$ | Least squares | PDE Solver |\n|--:|:--:|:--:|:--:|:---:|:--:|\n|KdV; $64^2$; 0.1 |   3.3 ms  |      249 $\\mu$s (7.5%)    |      2.8 ms (84%)    |     547 $\\mu$s (16%)   |    98 $\\mu$s (2.9%)  |\n| KdV; $512^2$; 0.1 |  102 ms|       7 ms (6.8%)     |      79 ms (77%)     |     3.7 ms (3.6%)   |   228 $\\mu$s (0.2%)   |\n|           KdV; $512^2$; 0.002 |      16 ms                 |       7 ms (43%)     |       9 ms (56%)     |     563 $\\mu$s  (3.5%)  |   246 $\\mu$s (1.5%)  |\n\n**Table 3.** Time taken by each loss function and the operations inside the dynamics loss for a single snapshot of 2D KdV measured by JAX profiler, with the portion of total training step reported in parentheses. Note that the operations above may contain additional compute overhead, which is not reported. We used a 24GB Nvidia A5000 GPU for these measurements.\n\n- **Training time.** We note that both reconstruction loss and the dynamics loss take significantly more time with increasing grid size. While we cannot control the decoding and reconstruction time present in both loss functions, we can control the hyper-reduction factor $\\gamma$. Note that by reducing $\\gamma$ for the $512^2$ case such that the reduced grid ($\\mathcal X^\\downarrow_{\\mathcal S}$) size is on the same order as the reduced grid for the $64^2$ case, the compute time of $L_{dyn}$ is reduced by about 90% and solving the least sqaures takes about the same time as the low-resolution case. Furthermore, the time taken by the dynamics loss is reduced from 84% of the total time in $64^2$ case to 56% when $\\gamma=0.002$. \n\n| PDE; Resolution; $\\gamma$ | Total $\\nabla_\\theta L_{\\Phi-ROM}$ | Only $\\nabla_\\theta L_{rec}$ | Only $\\nabla_\\theta L_{dyn}$ |\n|--:|:----:|:-:|:--:|\n|KdV; $64^2$; 0.1 |668 MB|376 MB|666 MB|\n|KdV; $512^2$; 0.1 |13106 MB|13106 MB|     1780 MB     |\n|KdV; $512^2$; 0.002 |2854 MB|1834 MB     |     1780 MB     |\n\n**Table 4.** Memory consumption for each loss function for one training snapshot. Note that the required memory for optimizing the total loss is less than the sum of individual losses, as JAX compiles the operations.\n\n- **Memory consumption.** As reported in Table 4, increasing the grid size increases the memory consumption. However, adjusting $\\gamma$ significantly limits the memory consumption by $L_{dyn}$ as well.\n\n- **Adjusting $\\mathbf{\\gamma}$.** As noted, reducing $\\gamma$ helped significantly with the training time and memory consumption. We add that hyper-reduced grids are sampled independently for each snapshot at every epoch. As such, even for small choices of $\\gamma$, the reduced grids would cover the spatial domain given sufficient training time. A lower bound for selection of $\\gamma$ is considering the latent size $k$, and selecting $\\gamma$ such that the size of the reduced grid is at least equal to $k$. Otherwise, the resulting least squares problem would be underdetermined.\n\n- **Solver and PDE complexity.** We conducted the same experiments with the same network size for the Navier-Stokes equation (2nd order PDE) and its FVM-based solver. While we observed similar time profiles, the FVM solver took 35 $\\mu$s ($64^2$) and 66 $\\mu$s ($512^2$), compared to the 98 $\\mu$s and 228 $\\mu$s of the KdV solver, respectively. These times are purely induced by the internals of each numerical solver, as we have not applied any performance improvements to them. While our training time is affected by the timing of the PDE-solver, relying on a differentiable PDE solver does not introduce a computational bottleneck for training $\\Phi$-ROM as is evident in Table 3.\n\n- **Scaling to 3D.** As training a 3D problem would require further hyperparameter tuning regarding the network and latent size, we are not able to train $\\Phi$-ROM for such a problem in a timely manner. Still, we believe the training algorithm should work well with 3D problems with adjustment of $\\gamma$, although it will require more computational resources and time compared to 2D. In addition to adjusting $\\gamma$, careful optimization of the code would also be essential for efficient scaling to 3D. For instance, our two loss functions are computed independently from each other, which means the full field is reconstructed twice in two separate forward passes (thus also requiring two backward passes) of the decoder, while they can be merged into one forward pass.\n\n## Q4: \nWhile we did not follow a specific set of rules for choosing $T_{tr}$ and $T_{te}$, we made some considerations for some datasets. In the case of Navier-Stokes, for instance, the solution is decaying in the long term. As such, we chose $T_{tr}$ such that the training would cover the development of the flow without fully entering the non-turbulent decayed regime (see Figure 11 in the Appendix). In the case of LBM, $T_{tr}$ is selected such that the training data covers the development of vortex shedding in the wake, which is followed by a periodic temporal pattern. Sensitivity to the choice of $T_{tr}$ can be seen through inclusion of such transitions in the flow regime in the training data. We also noted that both $\\Phi$-ROM and DINo failed in long-term temporal extrapolation for LBM, and thus $T_{te}$ covers a short range after $T_{tr}$ (see Figure 7.e, showing LBM error up to $T=200$)."}, {"review_id": "BFQpJEidEp", "reviewer": "Reviewer_MXgS", "summary": "This paper proposes \\Phi-ROM, a physics-informed reduced-order modeling (ROM) framework that leverages differentiable numerical solvers to guide the learning of latent space dynamics for time-dependent and parameterized PDEs. Unlike prior ROM approaches that treat numerical solvers as mere data generators, \\Phi-ROM integrates them into the training loop, ensuring that the learned latent dynamics remain consistent with the discretized governing physics. \n\nThe model architecture builds on an existing framework, using an implicit neural representation (INR) decoder and a dynamics network to evolve latent representations over time. By using automatic differentiation (AD) and hyper-reduction, \\Phi-ROM achieves computational efficiency while being mesh-agnostic, allowing training and inference on different spatial discretizations. The authors demonstrate \\Phi-ROM's performance on a variety of 1D and 2D PDEs using diverse solvers (finite difference, finite volume, and LBM), and show improved generalization to unseen parameters and sparse observations compared to baseline ROM and PINN-style models.", "questions": "- All test cases involve constant-coefficient PDEs. Can the authors demonstrate performance on PDEs with spatially varying diffusivity, anisotropic tensors, or stiff source terms?\n\n- What is the bottleneck in extending this approach to 3D or higher-resolution simulations? Can the authors comment (or been better, demonstrate) an example in 3D (see next point)\n\n- Given the reliance on Jacobians and inversion (even with hyper-reduction), does \\Phi-ROM scale to high-dimensional (e.g., N ~ 10⁶) settings (of mesh size)? What are the memory/runtime implications?\n\n- Can the method support adaptive spatiotemporal sampling during inference? Sampling more near high gradients, for instance?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "*Strength*\n\n- The paper contributes a new category of ROMs that use differentiable solvers, moving beyond traditional projection-based or PINN-augmented methods. The use of differentiable solvers as part of the training process addresses a key limitation of purely data-driven ROM (their divergence from physical constraints during long-time evolution).\n\n- The coupling of hyper-reduction with latent dynamics training is non-trivial and is a potentially novel contribution in the ROM literature.\n\n- The empirical comparisons are extensive and systematic. The authors test their model across various PDE types, solvers, and settings (dense vs. sparse, interpolation vs. extrapolation).\n\n- The mesh-free (or rather the mesh-agnostic) property of the method is particularly appealing for scientific computing pipelines where spatial discretizations evolve over time (e.g., AMR, or moving meshes).\n\n*Weakness*\n- The set of benchmark problems is limited to relatively simple PDEs with constant coefficients and low spatial dimensionality (1D and 2D). These are useful for validating the method but fall short of demonstrating its utility in complex, real-world scenarios, which is arguably the key motivation for such physics-informed ROMs.\n\n- While the components (INRs, differentiable solvers, latent dynamics) are not new individually, their composition in this context is original. Still, more emphasis could be placed on how this method fundamentally differs from DINo or PINN-ROM beyond training strategy.", "was_revised": true, "final_justification": "This is a good paper (as evidenced by the generally positive reviews by all reviewers). The authors have also done a good job of adding additional results and clarifying concerns. \nThis is a good paper, and I keep my positive score.", "rebuttal": "We thank the reviewer for the detailed review and insightful comments. Below, we address the questions (Q) and the weaknesses (W) raised in the review.\n\n## Q1 and W1:\nSupporting parameters that are functions of space is an interesting direction for improvement. We must acknowledge that, in the current form of our $\\Psi$ network, $\\beta$ can be either a scalar or a fixed-size vector, and not a function of space. To extend the concept to spatially varying parameters, one may replace $W_\\beta$ in the dynamics network (see Figure 6 in Appendix E.1) with an encoder, and the encoding of the parameter would then be concatenated with $\\alpha$. The new parameter encoder can simply be a CNN if the parameter is always on a fixed grid, or a mesh-free network, such as transformer-based encoders or PointNets. We believe that this addition to the model should work well with the current training method; however, due to time and resource constraints, we are not able to validate these ideas for this rebuttal.\n\n## Q2 and Q3:\nWe understand the scaling concerns regarding the 3D and high-resolution problems. To provide a fuller picture of computational implications of such problems, we profiled $\\Phi$-ROM's training costs for a high-resolution $512\\times512$ 2D KdV problem ($|\\mathcal X_\\mathcal S| = 2^{18}$). \n\n| PDE; Res.; $\\gamma$ | Total $\\nabla_\\theta L_{\\Phi-ROM}$ | Only $\\nabla_\\theta L_{rec}$ | Only $\\nabla_\\theta L_{dyn}$ | Least squares | PDE Solver |\n|--------------------------:|:-----------------------------------:|:---------------:|:---------------:|:-------------:|:----------:|\n|              KdV; $64^2$; 0.1 |                 3.3 ms                 |      249 $\\mu$s (7.5%)    |      2.8 ms (84%)    |     547 $\\mu$s (16%)   |    98 $\\mu$s (2.9%)  |\n|             KdV; $512^2$; 0.1 |                102 ms                |       7 ms (6.8%)     |      79 ms (77%)     |     3.7 ms (3.6%)   |   228 $\\mu$s (0.2%)   |\n|           KdV; $512^2$; 0.002 |                 16 ms                 |       7 ms (43%)     |       9 ms (56%)     |     563 $\\mu$s  (3.5%)  |   246 $\\mu$s (1.5%)  |\n\n**Table 1.** Time taken by each loss function and the operations inside the dynamics loss for a single snapshot of 2D KdV measured by JAX profiler, with the portion of total training step reported in parentheses. Note that the operations above may contain additional compute overhead, which is not reported. We used a 24GB Nvidia A5000 GPU for these measurements.\n\n- **Training time.** We note that both reconstruction loss and the dynamics loss take significantly more time with increasing grid size. While we cannot control the decoding and reconstruction time present in both loss functions, we can control the hyper-reduction factor $\\gamma$. Note that by reducing $\\gamma$ such that the reduced grid ($\\mathcal X^\\downarrow_{\\mathcal S}$) size is on the same order as the reduced grid for the $64^2$ case, the compute time of $L_{dyn}$ is reduced by about 90% and solving the least sqaures takes about the same time as the low-resolution case. \n\n\n| PDE; Resolution; $\\gamma$ | Total $\\nabla_\\theta L_{\\Phi-ROM}$ | Only $\\nabla_\\theta L_{rec}$ | Only $\\nabla_\\theta L_{dyn}$ |\n|--------------------------:|:-----------------------------------:|:---------------:|:---------------:|\n|              KdV; $64^2$; 0.1 |                668 MB               |      376 MB     |      666 MB     |\n|             KdV; $512^2$; 0.1 |               13106 MB              |     13106 MB    |     1780 MB     |\n|           KdV; $512^2$; 0.002 |               2854 MB               |     1834 MB     |     1780 MB     |\n\n**Table 2.** Memory consumption for each loss function for one training snapshot. Note that the required memory for optimizing the total loss is less than the sum of individual losses, as JAX compiles and optimizes the operations.\n\n- **Memory consumption.** As reported in Table 2, increasing the grid size increases the memory consumption. However, adjusting $\\gamma$ significantly limits the memory consumption by $L_{dyn}$ as well.\n\n- **Adjusting $\\gamma$.** As noted, reducing $\\gamma$ helped significantly with the training time and memory consumption. We note that hyper-reduced grids are sampled randomly and independently for each snapshot at every epoch. As such, even for small choices of $\\gamma$, the reduced grids would cover the spatial domain given sufficient training time. A lower bound for selection of $\\gamma$ is considering the latent size $k$, and selecting $\\gamma$ such that the size of the reduced grid is at least equal to $k$. Otherwise, the resulting least squares problem would be underdetermined.\n\n- **Scaling to 3D.** As training a 3D problem would require further hyperparameter tuning regarding the network and latent size, we are not able to train $\\Phi$-ROM for such a problem in a timely manner. Still, we believe the training algorithm should work well with 3D problems with adjustment of $\\gamma$, although it will require more computational resources and time compared to 2D. In addition to adjusting $\\gamma$, careful profiling and optimization of the code would also be essential for efficient scaling to 3D. For instance, our two loss functions are currently computed independently from each other, which means the full field is reconstructed twice in two separate forward passes of the decoder (thus requiring two backward passes). Instead, both losses can be computed in a single forward pass of the decoder. \n\n## Q4:\nAs the reviewer has nicely highlighted, \"The mesh-free (or rather the mesh-agnostic) property of the method is particularly appealing for scientific computing pipelines where spatial discretizations evolve over time (e.g., AMR, or moving meshes).\" This is absolutely correct because any given training snapshot may live on an arbitrary grid. This capability remains effective at inference time as well because $\\Phi$-ROM evolves the dynamics in the reduced latent space with its own discovered coordinates $\\alpha$. \n\nNote that spatial sampling during inference only happens at the inversion step for computing $\\alpha^0$, and since the decoder is continuous and mesh-free, it can perform both the inversion and reconstruction with any spatial discretization. Furthermore, the dynamics network (and thus $\\Phi$-ROM) is also continuous in time, and in fact, we used Bogacki--Shampine's 3/2 ODE solver with adaptive step size for integration.\n\n## W2:\nWe appreciate the reviewer's recognition of the novelty of our work. To reflect on the reviewer's point about highlighting the novelties beyond the training, we would like to point out that the addition of the PDE parameters to the existing ROM framework is new in this work and was not proposed in those cited works."}, {"review_id": "z9sbibaLOX", "reviewer": "Reviewer_Sn8d", "summary": "The paper proposes a Physics-informed Reduced-Order Model that builds on the DIno framework and directly trains on latent dynamics using the labels from the differentiable PDE solver. Numerical results show the advantages of this method compared to DIno and other physics-informed ROMs.", "questions": "Q1) How well does the method perform compared to non-ROM methods? Including some non-ROM baselines (like Physics-informed neural operator [1] and Physics-informed deep operator network [2]) may help readers get a more comprehensive sense of the position of the proposed method in the broader field of physics-informed machine learning. \n\n[1] Li, Zongyi, et al. \"Physics-Informed Neural Operator for Learning Partial Differential Equations\"\n\n[2] Wang, Sifan, et al. \"Learning the solution operator of parametric partial differential equations with physics-informed DeepONets.\"\n\nQ2) Why we observe smaller forecasting errors in $[T_{tr}, T_{te}]$ compared to in $[0, T_{tr}]$ in Diffusion case of Table 1? Typically, model predictions on the training dataset should have a smaller prediction error compared to the test dataset. \n\nQ3) Could you further explain the reasoning for the choice of $\\lambda$, which depends on the sensitivity of the differentiable solver to the noise and errors in the reconstructed fields?", "limitations": "Yes, adequately described in the paper.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "[Strengths]\n\nS1) Numerical experiments are extensive, including five representative problems and five baselines. \n\nS2) Comprehensive appendix with computation details, dataset description, hyperparameter setups, and visualizations.\n\nS3) The proposed method generalizes much better compared to DINo.\n\n[Weaknesses] \n\nW1) The computation of $\\dot{\\alpha}^{\\star}$ brings restrictions to the methods (requiring differential PDE solvers) and additional computation cost. The tradeoff between accuracy and computational cost for different baselines considered in this paper has not been fully discussed and compared.\n\nW2) There is no guidance on how to choose the reduced latent dimension $k$. The error brought by hyper-reduction for solving $\\dot{\\alpha}^{\\star}$ could be hard to quantify, resulting in learning inaccurate latent dynamics.", "was_revised": true, "final_justification": "I'd like to maintain my score after reviewing the authors' rebuttal. Regarding Q1, although the results of the FNO trained with a data-driven loss are interesting, I think for a fairer comparison with Phi-ROM, FNO should be trained with both a data-driven loss and a PDE-residual loss (i.e., PINO), since Phi-ROM also has access to the solver. Additionally, without seeing quantitative results addressing the identified weaknesses (such as the exact wall-clock time for computing $\\dot{\\alpha}^*$, training times for different methods, and results on the cost–accuracy trade off under different latent dimensions $k$), I would like to keep my rating.", "rebuttal": "We would like to thank the reviewer for the detailed review and the insightful comments. Below, we address the questions (Q) and the weaknesses (W) raised in the review.\n## Q1:\nWe greatly appreciate the reviewer’s feedback. We were hesitant to provide direct comparisons with non-ROM neural surrogates during the preparation of our submission, as we believed such a comparison would not be an apples-to-apples comparison. While it is true that both reduced-order models based on neural networks and other neural surrogates (e.g., neural operators) aim to accelerate the solution of PDEs, they do so in fundamentally different ways. The former evolves the dynamics in a reduced subspace and coordinate system, whereas the latter operates in the full space. This distinction becomes especially significant at scale, particularly for high-resolution and 3D problems. \n\nThat being said, to address the comments raised by the reviewer, we have conducted new experiments to compare $\\Phi$-ROM with two non-ROM neural surrogates by just focusing on the accuracy of the unrolling predictions.\n\nWe trained autoregressive FNOs with one-step training [1] and pushforward training (FNO-PF) [2] for the Burger's and LBM problems reported in the paper. We used the Apebench library [3] for training the FNO networks. As shown below for both problems, $\\Phi$-ROM outperforms FNO in temporal interpolation and extrapolation tasks. \n\n|            | $T_{tr}$ | $T_{te}$ |\n|-----------:|:--------:|:--------:|\n| $\\Phi$-ROM |   0.021  |   0.028  |\n|        FNO-1 |   0.037  |   0.131  |\n|     FNO-PF |   0.049  |   0.089  |\n\n**Table 1.** Burger's unrolling errors for test samples. $\\Phi$-ROM errors are the same as those reported in Table 1 of the paper. We used a 3-layer FNO with 32 hidden channels, 12 Fourier modes, and GELU activation. \n\n|            | $T_{tr}$ | $T_{te}$ |\n|-----------:|:--------:|:--------:|\n| $\\Phi$-ROM |   0.065  |   0.150  |\n|        FNO-1 |   0.145  |   4.867  |\n|     FNO-PF |   0.234  |   3.276  |\n\n**Table 2.** LBM unrolling errors for test samples (parameter interpolation). $\\Phi$-ROM errors are the same as those reported in Table 3 of the paper. We used a 4-layer FNO with 64 hidden channels, 16 Fourier modes, and GELU activation.\n\n[1] Li, Zongyi, et al. \"Fourier neural operator for parametric partial differential equations.\" arXiv preprint arXiv:2010.08895 (2020).\n\n[2] Brandstetter, Johannes, Daniel Worrall, and Max Welling. \"Message passing neural PDE solvers.\" arXiv preprint arXiv:2202.03376 (2022).\n\n[3] Koehler, Felix, Simon Niedermayr, and Nils Thuerey. \"APEBench: A benchmark for autoregressive neural emulators of PDEs.\" Advances in Neural Information Processing Systems 37 (2024): 120252-120310.\n\n## Q2:\nThanks for pointing this out! We believe that the higher error for the training temporal range in the diffusion case is due to a higher reconstruction error at earlier time steps. Note that the initial conditions for the Diffusion equation are generated from 2D Gaussian distributions with one peak of varying location, standard deviation, and magnitude. As the PDE evolves, the blob is *diffused* and flattened over the whole domain. This diffusion process means that the decoder would have an easier job in reconstructing the solutions at later times. Similar phenomena are observed in PINNs, where the network tends to learn the low-frequency parts of the solution more easily. See spectral bias in neural networks [1] and PINNs [2]. \n\nTo further validate this, we perform inversion on two snapshots of the same trajectory at $T=0$ and $T=100$ to examine the reconstruction error independently from the temporal error accumulation. The reconstruction error (RNMSE) at $T=0$ is $0.025$ while it goes down to $0.004$ at $T=100$. \n\n[1] Cao, Yuan, et al. \"Towards understanding the spectral bias of deep learning.\" arXiv preprint arXiv:1912.01198 (2019).\n\n[2] Wang, Sifan, Hanwen Wang, and Paris Perdikaris. \"On the eigenvector bias of Fourier feature networks: From regression to solving multi-scale PDEs with physics-informed neural networks.\" Computer Methods in Applied Mechanics and Engineering 384 (2021): 113938.\n\n## Q3:\nThe hyperparameter $\\lambda$ balances the two loss terms in the training objective. Note that while the latent manifold is constructed by $L_{rec}$, it is regularized by $L_{dyn}$ in order to enforce the PDE. Consequently, smaller $\\lambda$ increases the regularization effect of $L_{dyn}$ on the decoder and the latent manifold, resulting in a higher reconstruction loss. At the same time, a very large choice of lambda (>0.9) may result in a disparity between the actual latent dynamics and what the dynamics network $\\Psi$ learns as the latent dynamics. Thus, we suggested the 0.5 to 0.8 range for the choice of $\\lambda$. A larger $\\lambda$ may improve the forecast error at very early time-steps in exchange for high error accumulation, while smaller values may deteriorate the reconstructions and result in wrong time derivatives from the solver. In addition, we would like to note that as we use normalized errors for both loss functions, the losses would be in the same order of magnitude when $\\lambda$ is close to 0.5 and decrease with similar rates. \n\n### LBM Case\nIn our experiments, $\\lambda=0.5$ worked well for all problems, except for LBM, which has a more complex geometry because of the discontinuity at the cylinder. Due to the imperfect reconstruction of the cylinder, we observed noise and unexpected artifacts in the time derivatives. To alleviate this issue, we increased $\\lambda$ for a smaller reconstruction loss, which helped with training the model and improving its generalization.\n\n### Burger's ablation study\nTo better demonstrate the effect of $\\lambda$, we trained Burger's with $\\lambda=0.1, 0.5, 0.9$ and report the results below. \n\n| $\\lambda$ | $L_{rec}$ | $L_{dyn}$ | Error $[0-T_{tr}]$ | Error $[0-T_{te}]$ |\n|----------:|:---------:|:---------:|:--------------:|:--------------:|\n|       0.9 |   0.002   |   0.081   | 0.137          | 0.150          |\n|       0.5 |   0.019   |   0.013   | 0.021          | 0.028          |\n|       0.1 |   0.176   |   0.011   | 0.459          | 0.613          |\n\n**Table 1.** Loss and error values for Burger's with different choices of $\\lambda$.\n\nAs explained earlier, increasing $\\lambda$ decreases the regularization effect and thus results in smaller reconstruction loss, but fails to minimize the dynamics loss. Conversely, a small $\\lambda$ increases the reconstruction error in a way that the dynamics learned from the solver may be wrong or noisy.\n\n## W1:\nWe would like to highlight that the costs associated with computing $\\dot \\alpha^*$ and the use of solvers are one-time costs that affect only the training time. Furthermore, while $\\Phi$-ROM has a higher training cost per sample compared to DINo and CROM, it greatly improves data efficiency (see Figure 4) and generalization to out-of-distribution parameters (see Table 3). We believe these properties, as well as the improved temporal extrapolation, justify the increased training cost. \n\nWe acknowledge that training $\\Phi$-ROM takes longer than the DINo and CROM baselines. CROM training involves only minimizing an auto-encoder reconstruction loss, which is significantly cheaper than both DINo and $\\Phi$-ROM training. However, CROM moves this training cost to inference by solving the PDE in the physical space during inference. DINo, on the other hand, unrolls the latent dynamics by integration during the training, which $\\Phi$-ROM avoids. Still, DINo's training is computationally cheaper than $\\Phi$-ROM.\n\nFor the KdV problem on a single 24GB Nvidia A5000 GPU, a single training epoch of DINo takes about 11 seconds, while $\\Phi$-ROM takes about 34 seconds. We note that we used different hardware for the experiments in the paper (as noted in Appendix F.2, along with the training times of $\\Phi$-ROM models). \n\n## W2.1: Choice of latent dim $k$\nWe note that apart from the complexity of the PDE, the choice of $k$ highly depends on the decoder architecture and its hyperparameters (e.g. activation function), as the latent manifold is the input to this network. As detailed in Appendix E, we used two different decoders adopted from CROM and DINo, and noticed that the Hyper decoder (from DINo) requires noticeably larger $k$. While we chose the latent dimensions close to those used in the CROM and DINo papers, we relied on a training with only the reconstruction loss (and without minimizing the dynamics loss) to determine the latent dimensions that achieve small reconstruction errors. \n\n## W2.2: Hyperreduction error\nWe understand the reviewer's concern regarding hyper-reduction. We note that hyper-reduced grids are sampled randomly and independently for each snapshot at every epoch. As such, the reduced grids would cover the dynamics of the whole spatial domain given sufficient training time."}, {"review_id": "gTnRhpvShF", "reviewer": "Reviewer_1DUk", "summary": "The paper presents a new network-based method for reconstructing and predicting dynamics governed by spatial-temporal PDEs. It combines two network models: 1) a network that captures the solution manifold with reduced spatial dimensions; 2) a network that learns the temporal evolution of the reduced coordinates on the solution manifold. The training loss of the second network requires a differentiable PDE solver that provides temporal derivatives as supervision. The paper evaluates the proposed method on five 2D PDE problems and compares it with several standard baselines.", "questions": "Regarding the technical method:\n1. The method’s dependence on a differentiable solver seems to arise from the dynamic loss. As Remark 3.2 suggests, one could replace this loss design with the true gradients computed from the underlying PDEs, like in PINNs. To me, this seems to be a better alternative as it avoids the dependency on a differentiable PDE solver and exploits the true time derivatives from the underlying dynamics. So, have you considered this option as suggested in Remark 3.2? Would this alternative achieve comparable or better performance than your proposed method?\n\n2. The hyperparameter lambda is a bit unnatural and mysterious to me. Could you provide an ablation study on the effects of different lambda values and some insights for choosing them in a new problem?\n\n3. The reconstruction and dynamic losses characterize the spatial and temporal errors, respectively. I wonder whether they could be viewed through the lens of some variational or weak forms of the underlying PDEs, which may shed more light on the physical meaning of lambda.\n\nRegarding the experiments:\n\n4. Are there any technical reasons that may prevent the proposed method from being applied to 3D, larger-scale PDE problems?", "limitations": "I appreciate that the limitation section is up front about its dependency on a differentiable solver.\n\nI would suggest the limitation section stress that the current experiments are limited to 2D problems with small grid resolutions.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. This paper’s writing is pretty good. In particular, I like the background and method sections, which achieve a good balance between telling the high-level story and explaining the low-level technical details.\n2. The design of the network method is nice and clean. There are a number of ways of incorporating neural networks into reduced-order modeling, and I feel the overall design of the two network models provides an intuitive way of capturing the reduced-order dynamics. It is also a fairly general framework that can be applied to various PDEs and discretization schemes.\n\nWeaknesses:\n1. The experiments only consider two-dimensional PDE examples with low-resolution discretization schemes, i.e., roughly a hundred by a hundred cells according to the appendix.\n2. The experiments seem to lack comparisons with numerical PDE solvers. The motivation and benefits of applying neural physics over using existing numerical solvers are unclear to me.", "was_revised": true, "final_justification": "I'd like to maintain my score after reading the rebuttal and the other reviews. The rebuttal successfully addressed Questions 1-3 in my original review. The technical content in their answers to Weakness 1/Question 4 and Weakness 2 looks mostly OK, but I hope the rebuttal could fine-tune its wording to avoid overstating the paper's strengths and potential in 3D cases. Some careful text editing would be sufficient.", "rebuttal": "We would like to thank the reviewer for the feedback and the insightful suggestions. Below, we address the questions (Q) and weaknesses (W) raised in the review. \n\n## Q1:\nWe did compare $\\Phi$-ROM with two other physics-informed ROMs, specifically PINN-ROM (as defined in Remark 3.2) and CROM, in Section 4.2, and the results are reported in Table 1 of the paper. As correctly pointed out, using a PINN-inspired loss would avoid the requirement for a differentiable solver. However, as also thoroughly studied in the PINN literature, PINN loss suffers from various failure scenarios when solving PDEs with high-order terms, non-linearity, or stiff source terms [1]. As shown in Table 1, $\\Phi$-ROM and PINN-ROM achieve similar performance in the linear Diffusion problem. However, as shown in Figure 3, PINN-ROM falls short in the non-linear Burger's equation and fails to model the shock. This behaviour is in line with typical failure modes of PINNs reported in the literature. \n\n[1] Krishnapriyan, Aditi, et al. \"Characterizing possible failure modes in physics-informed neural networks.\" Advances in neural information processing systems 34 (2021): 26548-26560.\n\n## Q2:\nThe hyperparameter $\\lambda$ balances the two loss terms in the training objective. Note that while the latent manifold is constructed by $L_{rec}$, it is regularized by $L_{dyn}$ in order to enforce the PDE. Consequently, smaller $\\lambda$ increases the regularization effect of $L_{dyn}$ on the decoder and the latent manifold, resulting in a higher reconstruction loss. At the same time, a very large choice of lambda (>0.9) may result in a disparity between the actual latent dynamics and what the dynamics network $\\Psi$ learns as the latent dynamics. Thus, we suggested the 0.5 to 0.8 range for the choice of $\\lambda$. A larger $\\lambda$ may improve the forecast error at very early time-steps in exchange for high error accumulation, while smaller values may deteriorate the reconstructions and result in wrong time derivatives from the solver. We also like to note that as we use normalized errors for both loss functions, the losses would be in the same order of magnitude when $\\lambda$ is close to 0.5 and decrease at a similar rate. \n\n### LBM Case\nIn our experiments, $\\lambda=0.5$ worked well for all problems, except for LBM, which has a more complex geometry because of the discontinuity at the cylinder. Due to the imperfect reconstruction of the cylinder, we observed noise and unexpected artifacts in the time derivatives. To alleviate this issue, we increased $\\lambda$ for a smaller reconstruction loss, which helped with training the model and improving its generalization.\n\n### Burger's ablation study\nTo better demonstrate the effect of $\\lambda$, we trained Burger's with $\\lambda=0.1, 0.5, 0.9$ and report the results below. \n\n| $\\lambda$ | $L_{rec}$ | $L_{dyn}$ | Error $[0-T_{tr}]$ | Error $[0-T_{te}]$ |\n|-:|:-:|:-:|:-:|:-:|\n|0.9 |0.002|0.081| 0.137| 0.150|\n|0.5 |0.019|0.013| 0.021| 0.028|\n|0.1 |0.176|0.011| 0.459| 0.613|\n\n**Table 1.** Loss and error values for Burger's with different choices of $\\lambda$.\n\nAs explained earlier, increasing $\\lambda$ decreases the regularization effect and thus results in smaller reconstruction loss, but fails to minimize the dynamics loss. Conversely, a small $\\lambda$ increases the reconstruction error in a way that the dynamics learned from the solver may be wrong or noisy.\n\n## Q3:\nWe are inspired by the reviewer's thinking that points to the existing similarities between the variational form of a PDE and the training procedure of $\\Phi$-ROM that incorporates a PDE residual in its dynamics loss. However, it is not immediately trivial to us how a quantity similar to $\\lambda$ would appear in that variational formulation that would help with further interpretation of $\\lambda$ beyond what we have provided above. We believe our coupled and joint loss identifies the PDE in a new coordinate system, one which could potentially be written by composition of the latent space and physical space in an Euler-Lagrange framework. This would be an interesting theoretical investigation for future work.\n    \n## Q4 and W1:\nWe understand the scaling concerns regarding the 3D and high-resolution problems. We believe there are no technical reasons to prevent $\\Phi$-ROM from being applied to 3D and larger-scale PDE problems beyond a natural need for larger compute resources. We will elaborate on this in what follows.\n\nTo provide a fuller picture of computational implications of such problems, we profiled $\\Phi$-ROM's training costs for a high-resolution $512\\times512$ 2D KdV problem ($|\\mathcal X_\\mathcal S| = 2^{18}$). \n\n| PDE; Res.; $\\gamma$ | Total $\\nabla_\\theta L_{\\Phi-ROM}$ | Only $\\nabla_\\theta L_{rec}$ | Only $\\nabla_\\theta L_{dyn}$ | Least squares | PDE Solver |\n|-:|:-:|:-:|:-:|:-:|:-:|\n|KdV; $64^2$; 0.1 |3.3 ms|249 $\\mu$s (7.5%)|2.8 ms (84%)|547 $\\mu$s (16%)|98 $\\mu$s (2.9%)|\n|KdV; $512^2$; 0.1 |102 ms|7 ms (6.8%)|79 ms (77%)|3.7 ms (3.6%)|228 $\\mu$s (0.2%)|\n|KdV; $512^2$; 0.002 |16 ms|7 ms (43%)|9 ms (56%)|563 $\\mu$s  (3.5%)|246 $\\mu$s (1.5%)|\n\n**Table 2.** Time taken by each loss function and the operations inside the dynamics loss for a single snapshot of 2D KdV measured by JAX profiler, with the portion of total training step reported in parentheses. Note that the operations above may contain additional compute overhead, which is not reported. We used a 24GB Nvidia A5000 GPU for these measurements.\n\n- **Training time.** We note that both reconstruction loss and the dynamics loss take significantly more time with increasing grid size. While we cannot control the decoding and reconstruction time present in both loss functions, we can control the hyper-reduction factor $\\gamma$. Note that by reducing $\\gamma$ such that the reduced grid ($\\mathcal X^\\downarrow_{\\mathcal S}$) size is on the same order as the reduced grid for the $64^2$ case, the compute time of $L_{dyn}$ is reduced by about 90% and solving the least sqaures takes about the same time as the low-resolution case. \n\n| PDE; Resolution; $\\gamma$ | Total $\\nabla_\\theta L_{\\Phi-ROM}$ | Only $\\nabla_\\theta L_{rec}$ | Only $\\nabla_\\theta L_{dyn}$ |\n|-:|:-:|:-:|:-:|\n|KdV; $64^2$; 0.1 |668 MB|376 MB|666 MB|\n|KdV; $512^2$; 0.1 |13106 MB|13106 MB|1780 MB|\n|KdV; $512^2$; 0.002 |2854 MB|1834 MB|1780 MB|\n\n**Table 3.** Memory consumption for each loss function for one training snapshot. Note that the required memory for optimizing the total loss is less than the sum of individual losses, as JAX compiles and optimizes the operations.\n\n- **Memory consumption.** As reported in Table 3, increasing the grid size increases the memory consumption. However, adjusting $\\gamma$ significantly limits the memory consumption by $L_{dyn}$ as well.\n\n- **Adjusting $\\gamma$.** As noted, reducing $\\gamma$ helped significantly with the training time and memory consumption. We note that hyper-reduced grids are sampled randomly and independently for each snapshot at every epoch. As such, even for small choices of $\\gamma$, the reduced grids would cover the spatial domain given sufficient training time. A lower bound for selection of $\\gamma$ is considering the latent size $k$, and selecting $\\gamma$ such that the size of the reduced grid is at least equal to $k$. Otherwise, the resulting least squares problem would be underdetermined.\n\n- **Scaling to 3D.** As training a 3D problem would require further hyperparameter tuning regarding the network and latent size, we are not able to train $\\Phi$-ROM for such a problem in a timely manner. Still, we believe the training algorithm should work well with 3D problems with adjustment of $\\gamma$, although it will require more computational resources and time compared to 2D. In addition to adjusting $\\gamma$, careful profiling and optimization of the code would also be essential for efficient scaling to 3D. For instance, our two loss functions are currently computed independently from each other, which means the full field is reconstructed twice in two separate forward passes of the decoder (thus requiring two backward passes). Instead, both losses can be computed in a single forward pass of the decoder. \n\nAs suggested by the reviewer, we will comment on the added computational costs associated with the 3D problems in our limitations section.\n\n## W2:\nWe would like to highlight that all our error metrics are computed by comparing the results of a reduced-order model and its full-order model (i.e. the numerical PDE solver).\n\nWhile numerical solvers offer accurate solutions with convergence guarantees, their application in scenarios that require real-time computation of the solutions is not practical, especially when dealing with high-resolution and 3D problems. Reduced-order models try to solve the computational burden of full-order solvers by operating in a low-dimensional solution manifold. Finally, neural ROMs construct non-linear manifolds by using neural networks, which is essential for modelling of nonlinear PDEs.\n\nCompared to a numerical solver, $\\Phi$-ROM operates in a fixed-size latent space. While the computational cost of a full-order numerical solver would rapidly increase with a high-resolution mesh or 3D problems, the inference speed of $\\Phi$-ROM depends on the significantly smaller latent size. \n\nTo provide a better picture of the difference between time-stepping in the physical space and latent space, we measured the integration time for 100 time-steps of the LBM problem using the numerical solver and $\\Phi$-ROM, both running on a GPU. $\\Phi$-ROM took **1.48 ms** to forecast the latent solution for 100 steps, while solving in the physical space took **19.9 ms**. Note that as the dimension of the problem grows, the gap between the reduced and full order models would grow as well."}, {"review_id": "u2gfviAVGx", "reviewer": "Reviewer_rGnK", "summary": "This paper presents a physics-informed framework for reduced-order modeling (ROM) of time-dependent partial differential equations (PDEs), with a focus on improving generalization and physical consistency. The authors propose Φ-ROM, a novel method that integrates differentiable PDE solvers into the learning pipeline, allowing the latent space dynamics to be directly shaped by the governing physics - this is unlike traditional data-driven ROMs, which often rely solely on high-fidelity simulations for supervision and may suffer from latent space drift. This alignment between latent and physical dynamics enables better forecasting, particularly under extrapolation to unseen initial conditions or parameter regimes. The authors demonstrate the effectiveness of their method across different PDE systems, showing improved long-term prediction accuracy, lower data requirements, and better physical interpretability.", "questions": "Two questions:\n- The authors mentioned \"When trained, the latent coordinates corresponding to a new state are computed by inversion\", \"The\nminimization problem in Eq. (4) is solved iteratively using a gradient-based optimizer.\" -- why didn't the author build an encoder to construct the \\alpha from u, which is probably more common in auto-encoder based ROM literatures?\n- For the description of Fig 1: \"The initial condition and the forecast solution at a target time T∗, can be observed and reconstructed on arbitrary (and possibly different) grids X and X′\" -- why the two grids could ever be different?", "limitations": "Yes. The author explicitly mentioned \"we acknowledge that requiring access to differentiable PDE solvers limits Φ-ROM’s immediate\napplicability\" which is absolutely true. And I don't think there's any major limitations for the method beyond that.", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths: The paper is technically sound, with a clear motivation and rigorous formulation of the proposed Φ-ROM framework. By integrating differentiable solvers into the latent dynamics, the authors bridge a critical gap between physics-informed modeling and data-driven ROMs. The training and evaluation procedures are well-defined, and experimental results appear thorough and convincing. The paper is generally well-written and logically organized. The abstract and introduction clearly articulate the motivation and novelty. Figures and results are presented in a digestible manner.\n\nWeakness: I have two major concerns. 1. While Φ-ROM introduces a physics-informed mechanism by incorporating differentiable solvers into the latent dynamics, it’s unclear whether this approach offers a fundamentally better alternative to PINN-ROM. PINN-ROM directly enforces the true PDEs during training—even if it performs worse empirically—its objective remains aligned with the exact continuous PDE. In contrast, Φ-ROM relies on a discretized solver whose numerical errors (local and global truncation errors, as understood from Taylor expansion analysis) mean the model is optimizing toward a discretization-dependent surrogate. This raises the concern that Φ-ROM may inherit the solver’s biases rather than learn the true physics. The paper does not sufficiently justify why this tradeoff (solver-consistency over PDE-consistency) should be preferred in principle. 2. The Φ-ROM architecture seems to fall into a regime that the original CROM paper explicitly discouraged—namely, building a low-dimensional representation (\\alpha) from already discretized vector fields. CROM advocated modeling the continuous dynamics of the field itself rather than that of a discretized approximation. Yet, this paper claims Φ-ROM outperforms CROM broadly, which introduces a tension: either Φ-ROM violates the key assumptions behind CROM, or CROM’s earlier recommendations do not hold universally. The authors do not discuss this contradiction, nor clarify what might enable Φ-ROM to work better in such a regime. A clearer reconciliation between these two frameworks is warranted.", "was_revised": true, "final_justification": "I have revisited the CROM paper and my confusions are resolved. The authors also successfully addressed my other concerns, I'll raised my rating to reflect the changes.", "rebuttal": "We would like to thank the reviewer for the detailed review and insightful comments. Below, we address the questions (Q) and weaknesses (W) raised in the review.\n## W1:\nWe appreciate the reviewer’s thoughtful feedback regarding the trade-off between solver-consistency and PDE-consistency. As the reviewer is rightly aware, standard numerical solvers are designed to be consistent with the continuous PDE; they converge to the true solution as the mesh is refined and the discretization error vanishes.\n\nIn our view, the fundamental issue with PINN-ROM is that it delegates the task of recovering the continuous PDE solution entirely to the neural network via automatic differentiation. In practice, this makes training significantly harder, because the network must approximate both the solution and the discretization implicitly—often without the benefit of decades of insights embedded in well-tested numerical schemes. This is one reason why PINN-ROMs (and PINNs in general) often struggle in practice, as the reviewer also observes and has pointed out in their comment.\n\nBy integrating a differentiable numerical solver directly into $\\Phi$-ROM, we shift this burden away from the network. The solver provides a reliable, convergent discretization of the PDE, allowing the neural network to focus its capacity on modeling the reduced latent dynamics instead of re-learning the physics from scratch. We agree that this introduces a dependency on the discretization’s properties, but we argue (as also supported by our results presented in the paper) that leveraging a well-understood, consistent solver is a principled and practical way to combine the strengths of numerical analysis and machine learning.\n\n## W2:\nThank you for pointing out the differences between CROM and $\\Phi$-ROM. We would like to raise some points regarding CROM to help explain the apparent contradictions indicated by the reviewer between CROM and $\\Phi$-ROM. First, CROM models the temporal dynamics of the PDE in the reconstructed physical space instead of the latent space. As such, a reduction in CROM is only achieved when *subsampling* is performed in the reconstructed field. Looking at the last table reported on the final page of the CROM paper, it is clear that spatial sub-sampling is absent in most of their reported experiment, meaning that inference has been on the full space, which is effectively similar to solving the full order PDE (i.e. no reduction). Second, CROM adopts two different approaches for modelling the dynamics in the reconstructed physical space: 1) A continuous model of the PDE similar to PINNs (denoted by \"AD-CROM\" in our paper), or 2) the discretized (e.g. finite difference-based) dynamics (denoted by \"FD-CROM\" in our paper). As acknowledged in the CROM paper, their computation of the network gradients using automatic differentiation for continuous physics modelling can yield erroneous gradients (as has also been studied in the PINN literature), and instead, a discretized model is used in such scenarios. It is therefore very important to note that subsampling is **not** possible (or would yield a high truncation error) when the discretized model (\"FD-CROM\") is used for time-stepping. As such, when CROM relies on a discretized model (\"FD-CROM\") due to the failure of their continuous model (\"AD-CROM\"), there is no subsampling reported, which implies no reduction has been achieved, and the PDE is essentially solved in a full-order manner using the discretized physics.\n\nIn our experiments, we considered both versions of CROM, with and without subsampling for AD-CROM. Note that in the Burger's experiment, FD-CROM outperforms all other models for the training temporal horizon by an order of magnitude (see Table 1 and Figure 8.b). This is because FD-CROM is effectively equivalent to solving the full order PDE. However, as continuous modelling yields high gradient error near the shock, AD-CROM and its subsampled version fail to model the PDE (see Figure 3.a). Furthermore, note that CROM is trained by minimizing only the reconstruction loss of an auto-encoder. As there are no other regularizing terms, the decoder fails to generalize to unseen patterns beyond the training time horizon, which explains the relatively poor temporal extrapolation performance of FD-CROM.\n\n## Q1:\nAs the reviewer has suggested, many works use auto-*encoders*, often with fully-connected or CNN encoders. Note that in such designs, even if the decoder network is an INR and mesh-free (as in CROM), the whole encoding and decoding pipeline is **not** mesh-free because of the encoder. As such, during training, all samples must be on the same, prescribed grid $\\mathcal X$ in order to be passed to the encoder, and during inference, the initial condition must again be on the same training grid $\\mathcal X$ for the encoder to work and compute $\\alpha^{t_0}$. Instead, by using inversion (or an auto-*decoder* similar to DINo), we have a completely mesh-free setup that does not impose any assumptions about the data grid and its consistency during training and inference. We address why this is important in our response to the next question (Q2). Nevertheless, we acknowledge that other encoder architectures, such as those based on PointNets and Transformers, may be adopted for a mesh-free setup, and we believe our physics-informed training strategy would work well with those designs too.\n\n## Q2:\nIn our experiments in Sections 4.3 and 4.4 of the paper, we demonstrated how the models can be trained and evaluated on sparse, irregular, and inconsistent observations of the input fields (e.g. the initial conditions) and still recover the full solution fields on the complete grid. In all these experiments, the input grid $\\mathcal X$ and the output grid $\\mathcal X'$ are different. This was possible firstly because of the mesh-free design of the network (as described in our response above to Q1), and secondly because of $\\Phi$-ROM's physics-informed strategy, which enables such training with a very sparse and minimal amount of data (as compared to DINo in Tables 2 and 3 of the paper). In practice, this ability is beneficial for data assimilation tasks, especially when only sensor input observations are available, but the full solution field is required."}]

Author rebuttal:

Meta-review: The paper introduces Φ-ROM, a physics-informed reduced-order modeling framework that integrates differentiable PDE solvers into the training loop, ensuring latent dynamics remain consistent with discretized physics. Reviewers praised the clear methodology, mesh-agnostic design, and extensive experiments across several PDEs, noting strong generalization to unseen parameters. The main limitations are the focus on relatively simple 1D/2D problems and dependence on differentiable solvers, which may hinder scalability to larger or 3D systems. Nonetheless, the consensus is that this represents a solid and original contribution to physics-informed ROMs.
Decision: Accept (poster)