Spectral-Refiner: Accurate Fine-Tuning of Spatiotemporal Fourier Neural Operator for Turbulent Flows
Authors: Shuhao Cao
Francesco Brarda
Rui Peng Li
Yuanzhe Xi
Venue: ICLR 2025 Poster (2025)
Keywords: operator learning
neural operators
Navier-Stokes
PDE
partial differential equations
computational fluid dynamics
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Combining the strengths of Fourier Neural Operators and the wisdom from traditional numerical methods to fine-tune the Spatiotemporal modification of FNO to get 10000 times better accuracy.

Abstract: Recent advancements in operator-type neural networks have shown promising results in approximating the solutions of spatiotemporal Partial Differential Equations (PDEs). However, these neural networks often entail considerable training expenses, and may not always achieve the desired accuracy required in many scientific and engineering disciplines. In this paper, we propose a new learning framework to address these issues. A new spatiotemporal adaptation is proposed to generalize any Fourier Neural Operator (FNO) variant to learn maps between Bochner spaces, which can perform an arbitrary-length temporal super-resolution for the first time. To better exploit this capacity, a new paradigm is proposed to refine the commonly adopted end-to-end neural operator training and evaluations with the help from the wisdom from traditional numerical PDE theory and techniques. Specifically, in the learning problems for the turbulent flow modeled by the Navier-Stokes Equations (NSE), the proposed paradigm trains an FNO only for a few epochs. Then, only the newly proposed spatiotemporal spectral convolution layer is fine-tuned without the frequency truncation. The spectral fine-tuning loss function uses a negative Sobolev norm for the first time in operator learning, defined through a reliable functional-type a posteriori error estimator whose evaluation is exact thanks to the Parseval identity. Moreover, unlike the difficult nonconvex optimization problems in the end-to-end training, this fine-tuning loss is convex. Numerical experiments on commonly used NSE benchmarks demonstrate significant improvements in both computational efficiency and accuracy, compared to end-to-end evaluation and traditional numerical PDE solvers under certain conditions. The source code is publicly available at https://github.com/scaomath/torch-cfd.

Published as a conference paper at ICLR 2025
SPECTRAL-REFINER ACCURATE -INIIG OF SPATIOTEMPORAL FOURIER NEURAL OPERATOR FOR TURBULENT FLOWS
Shuhao Cao"
Francesco Brarda
School of Science and Engineering University of Missouri-Kanse City
Department of Mathematics
Emory University
Rui Peng Li
Yuanzhe Xi
Center for Applied Scientific Computing Lawrence Livermore National Laboratory
Department of Mathematics Emory University
ABSTRACT
Recent udvancements in operator-type neural networks have shown promising re- sults in approximating the solutions of spatiotemporal Partial Differential Equations (PDEs). However, these neural networks often entail considerable training expenses, and may not always achieve the desired accuracy required many scientific and engineering disciplines. In this paper, we propose new learning framework to address these issues. A new spatiotemporal adaptation is proposed to generalize any Fourier Neural Operator (FNO) variant to learn maps between Bochner spaces, which can perform an arbitrary- length temporal super-resolution for the first time. To better exploit this capacity, new paradigm is proposed refine the commonly adopted end-to-end neural operator training and evaluations with the help from the wisdom from traditional numerical PDE theory and techniques. Specifically, in the learning problems for the turbulent flow modeled by the Navier-Stokes Equations (NSE), the proposed paradigm trains FNO only for few epochs. Then, only the newly proposed spatiotemporal spectral convolution layer is fine-tuned without the frequency truncation. The spectral fine-tuning loss function uses a negative Sobolev norm for the first time in operator learning, defined through reliable funetional-type a posteriori error estimator whose evaluation is exact thanks to the Parseval identity. Moreover, unlike the difficult nonconvex optimization problems in the end-to-end training, this fine-tun tenettuning loss is convex. Numerical experiments commonly used NSE benchmarks demonstrate significant improvements in both computational efficiency and accuracy, compared to end-to-end evaluation and traditional numerical PDE solvers under certain conditions. The source code is wwwettttttottttteeeeeeeeeeeeeeee
publicly available at https: / thub.
INTRODUCTION
Recently, Deep Learning (DL) pipelines have proven particularly effective in addressing problems formulated by Partial Differential Equations (PDEs). In this paper, we study the problem of learning Neural Operators (NOs) between infiniteddimensional function spaces (Boullé Townsend, 2023; Kovachki et al., 2023; Azizzadenesheli al., 2024; de Hoop et al., 2022: 2023). Specifically, the problem in consideration is for the Pavier-Stoker Equations (NSE) in the turbulent regime (Re O(103) O(10*)). In computation, the difficulties of solving NSE in this regime roots from its "stiffness" attributed to the nonlinear convection with nearly singular viscous diffusion, well the spatiotemporal nature of being highly transient. For this problem, propose to synergize operator learning with classical numerical PDE methods, complementing one's drawbacks and limitations with the other's strengths.
Corresponning scao@unkc. yuanzhe... x xi@emory.ede edu

---

Published as a conference paper at ICLR 2025
Compromises and drawbacks of traditional numerical methods. To solve NSE, traditional time stepping schemes include Adams- Bashforth/-Mouln or Runge-Kutta (RK) families, e.g., see Canuto et al. (2007, Appendix D), Karniadakis et al. (1991). If one opts to use explicit scheme, or there exists a certain portion of the foreing terms (e.g.. the convection term in NSE) computed via explicit schemes, then small time step sizes must be used: At ~ O (( ((r)))), 2 (Rannacher, 2000, Chapter 4). This necessity is often referred to as the "stiffness" of the PDE, and the threshold or constraint of the time step called the Courant- Friaatr-rFeerrcsseeeww (CFL) condition, e.g., see Johnston & Liu (2004); Wang (2012); He & Sun (2007). CFL poses a sufficient condition on the step size for "stability"", and this requires the step size usually much smaller than the requirement for "consistency". Here, this temporal consistency usually refers to first-order optimal local truncation error, e.g., how the original Butcher tableau is derived for RK (Butcher, 1965). The CFL. puts a threshold for all explicit mime-marching schemes on how fast the local errors in different frequencies can propagate, and thus prevent error accumulation "pollutes" the approximation globally to ensure stability. This means that, for any time-marchin schemes, finer mesh (better spatial consistency) requires the time steps to be much smaller to prevent errors from traveling to neighboring nodes and elements, which greatly increases the computational cost.
Among the end-to-end operator learning for NSE, the common approach is the so-called autoregressive "roll-out". During rolling-out, several snapshots of solutions are concatenated as channels in the input tensor to the neural operator, which outputs an approximated solution at the subsequent time step. This procedure can be repeated recurrently until reaching the model's stability capacity. In contrast to the traditional time marching solvers whose step sizes are restricted by the CFL condition, roll-out can withstand much bigger time step size, NOs used roll- out approach include the original FNO2d (Li et al., 2021), and those in Li et al. (2022); Brandstetter et al. (2023); Gupta & Brandstetter (2023); Fonseca et al. (2023); Zheng et al. (2024). However, for an end-to-end operator learner, the roll-out approach faces the same dilemma of super-resolution capacity in the temporal dimension as a traditional solver: finer time steps cost prohibitively more, while larger time steps does not guarantee stability. To solve this problem. we turn to the prevalent functional analytic framework for studying the convergence and stability of solution trajectories for NSE: the theory of Bochner spaces, e-g.. Aubin-Lions lemma (Lions & Magenes, 2012, Chapter 3), (Temam, 1995, Chapter 2), (Evans, 2022, Chapter 7). Inspired by this theoretical insight, we propose a SeapioTemporal adaption for all Fourier Neural Operator (ST-FNO) variants. ST-FNO now can perform arbitrary zero-shot super-resolution not only spatial dimensions, but also allows the temporal dimension to vary for the first time. For a prototype spatiotemporal operator learner FNO3d in Li et al. (2021), it learns map between fixed number of spatial "snapshots" (product spaces). The newly proposed ST-FNO directly learns a map between Bochner spaces LT V) to L' V) on non-overlapping time intervals To and T. This makes the model become trraecctor-tt--rrjeectto where the initial trajectory is the input of an NO to obtain the output evaluation as an approximation to the subsequent trajectory of the solution. Here V denotes the spatial Hilbert space in which snapshots of the solution at specific time reside, and for detailed notation list, we refer the readers to Appendix A.
Spatiotemporal operator learning.
Limits and the lack of accuracy for neural operators. The NO approach has the potential to bypass various difficulties and compromises of traditional schemes of numerical PDEs. However, in all end-to-end operator learning benchmarks of NSE, even the oatee-of-the-al NOs still fall short in achieving high-accuracy solutions. For example, an end-to-end roll-out approach suffers from error accumulation experimentally, e.g., see Li et al. (2022, Figure 9). To our best knowledge, in terms of the relative difference with the ground truth under the Bochner norm, no NO-only-based operator learning approach can break the barrier of an even 2-digit accuracy. Moreover, to our best knowledge, NO-only approaches have no theoretical stability estimate, such as the error propagation operator contraction. Recently, a remarkable advancement called PDE-Refiner (Lippe al., 2023) learns an extra error correction NO under the Denoising Diffuusion framework (Ho et al., 2020). For single instance, it can get stable long roll-out, and achieve for the first time O(10~8) relative difference with the ground truth after single time marching step, and (11066) in the Bochner norm. Nevertheless, for all models above and their learning frameworks, the optimization is to minimize the difference between the outputs from the NO, namely uv, and the ground truth U.s generated by a traditional numerical PDE solver. The difference with the true solution u under a certain norm is not directly optimized, applications. For difficult PDEs such as the NSE, in a linear time-stepping scheme, the ground truth us may on its own only have 3-digit accuracy in the Bochner norm. Note that this may already
the analytical expression of u is unknown most of the time in real-life
2

---

Published as a conference paper at ICLR 2025
be of optimal order ((tt) by convergence estimates (Heywood & Rannacher, 1986). Therefore, minimizing the difference between u. and üs becomes fruitless if the numerical approximation (ground truth) Ug is not accurate at the first place, as unnecessary computational resources may have been spent to get closer to us.
New hybrid scheme. To address these dilemmas for NSE's approximation, we take an alternative hybrid learning paradigm inspired by traditional numerical methods. Unlike the arduous training of running hundreds of epochs, the newly proposed ST-FNO is trained for only few epochs (e.g.. 10), concluding with the freezing of most model parameters. Then, the newly proposed spatiotemporal spectral convolution, as the last layer, in ST- is fne-tuned with the help of traditional solvers. Note that SF-FNO's ' temporal arbitrary-length inference capacity is attributed to this last layer. During fine-tuming, this layer is relieved from the common frequency truncation FNO has. A traditional solver with a single time step is used to obtain highly accurate approximations for extra field variables to conform with the physics. These extra fields help the evaluation of new loss function in fine-tuning a functional-type a posteriori error estimate measured under a negative Sobolev norm, which is equivalent the variational form-associated norm for the NSE, Note that there shall be no training of extra eror-correcting NOs (Dresdner et al., 2023; Lippe et al., 2023). Moreover, unlike the same nonconvex optimization problem that the extra error- correcting NOs try to tackle as the one in the original end-to-end training. our fine-tuning optimization problem is convex, and allows to achieve high accuracy with a fractional of computational resources when compared with other fine-tuning approaches that optimize the whole model with a physics-informed loss.
Fine-tuning using functional-type a posteriori error estimation. approximated solutions that are closer to the true solution directly, we turn to the a posteriori error estimation techniques. This technique allows computing the error without knowing the true PDE solutions, and has been widely studied for Galerkin-type methods, such as finite element (Ainsworth & Oden, 1993; 1997; Oden al., 1994), as well as for parameterized PDEs (Hesthaven et al., 2016; Patera et al., 2007). Among all types of posteriori error estimation techniques, functional-type posteriori estimator (Repin, 2008) views the error as a functional on the test Sobolev spaces and evaluates accurate representations with the help of an extra dual variable (Ern & Vohralík, 2010). In our hybrid approach, inspired by this, we combine the strengths of NOs and traditional numerical solvers. Using the newly proposed spatiotemporal Discretization-inaare invariant NO, we propose to use negative Sobolev norm in the frequency domain as a functional-type a posteriori error estimation as the fine-tuning loss. Traditionally, the a posteriori error estimation is used to refine local basis, yet this constraint on the purpose of the method leads to inaccurate localized representations for the H --1 error functional. Our new approach is not attached to the local refining requirement (Bonito et al., 2024). As result, the negative Sobolev norm used in the new loss is handily defined through the Gelfand triple in the frequency domain which global. Our method needs no extra dual variable reconstruction problem as in the traditional methods (Em & Vohralík, 2010). Unlike commonly adopted physics-informed operator learning (Li et al., 2024d), NO prediction for other field variables are not needed either. The extra field variables for computing the error are recovered through an -uto-differentiable numerical solver with the NO output of the primal variable (vorticity or velocity). Overall, this practice leads to "refining" the learnable set of the global spectral basis.
To seek highly accurate NO
For a more detailed review of operator learning and further motivation discussion with a much higher degree of mathematical rigor, we refer the readers to Appendix B.
The main contributions of this work are summarized as follows. For the difference between the common roll-out approach versus direct spatiotemporal learning between Bochner spaces, we refer the readers to Figure 1.
Main contributions.
Spatiotemporal Fourier Neural Operator.
We design the first spatiotemporal adaption technique for all FNO variants (ST-FNO) to enable them to learn maps between Bochner spaces.
New hybrid operator learning paradigm. We propose to train and evaluate ST-FNOs using a new strategy, which has sigmificantly improved over the existing methods in speed and accuracy. Only a few epochs of training combined with a spectral-refining fine-tuner yield highly accurate approxi- mations. Extra field variables are obtained by -uto-differnntiable mongle-step (RK4) or multistep (Adam-Bashe Bashforth) solvers, the input of which are neural predictions. The auto-differeniiabiiiyy of the solver makes it possible for the fine-tuning to optimize the parameters in FNO.
Functional a posteriori error estimation. Leveraging the spectral structure of the new Spatio
poral layer in ST-FNO and Parseval identity, the fine-tuning utilizes for the first time the negative
3

---

Published as a conference paper at ICLR 2025
min R ux. Digaluxr),
22 where uxx (10 (0)m))))
STENO
RK1
Rollou
Figure Schematic differences between approaches. 4th-order Runge-Kutta (RK4): small time steps bounded by the CFL condition, de-aliasing filter needed: autoregressive NO rolling-out: using previous evaluation as input repetitively, large time steps. no stability guarantees; Spatiotemporal FNO (ST-FNO) with hybrid fine-tuning: large time steps, yielding arbitrary-length temporal prediction in single evaluation, parallel-in-ting optimization for spectral fine-tuning.
Sobolev norm (functional norm) as loss via posteriori enror estimation. This procedure does not require any ground truth data (c.g.. numerical solution generated by traditional numerical solvers), nor the knowledge of the true solution(s) to the PDE. This new loss proven to be reliable in theory, in the meantime much more efficient in the ablation experiments.
Experimental results. We created native PyTorch port of Google's Computational Fluid Dynamics written in JAX (Kochkov et al., 2021; Dresdner et al., 2023) with enhanced functionality for tensor operations such as the facilitation of fine-tuning for the latent fields, publicly available at https: / gi thub. com/ ccoma wechhthh--eeeeeeeeoo cfd, with scripts to replicate the experiments well as the data generation. The data are available navier-stokes dataset..
2 SPATIOTEMPORAL OPERATOR LEARNING FOR HAVIER-STOKER EQUATIONS
Both drawbacks and advantages traditional numerical methods and NO-based methods in Section 1 play vital roles in shaping our study. We first briefly diseuss the spatiotemporal operator learning problem on Bochner spaces associated with NSE. Then, we detail how to modify generic FNO-based neural architecture obtain an operator learner between Bochner spaces.
2. SPATIOTEMPORAL OPERATOR LEARNING PROBLEM NSE
For 2D NSE, in the the velocity-pressure pressure (V-P) formulation (2.2), the velocity field u(t, [0, x + 2 seen as an element in the Bochner space [[, T). V) where a spatial Sobolev space in which each snapshot at / of the solution uut,-) -) resides. As in Li et al. (2021), we also consider the arricitt-streamffunctio (V-S) formulation (2.1) with periodic boundary condition (PBC). In [0, T Ω2 the scalar- scalar-valued valued vorticity := V and streamfunction & satisfies u rot vb. These two formulations read
VVrrtiity-Streannfuncii
0,w + rot db Vi - vAw = Vxf,
w A = 0.
(2.1)
WeelocityPPressure
0,u (u V)u vAu Vp f.
v 0.
(2.2)
For all analyses in line with the Hilbertian framework, V = H' (T2) for vorticity and H (T2) for velocity, where 2 denotes the unit torus, i.e., Ω (0, 1)2 with component-wise PBC. For fixed forcing function in V", the initial condition is either w00, x) wb(x) or u(0. -) rot (oo((-A)=)) (-))- Here Wp drawn from compactly supported probability measure μ., in which the compactness corre- sponds certain powerlenstrophy spectrum decay law to produce isotropic turbulence (McWilliams,
4

---

Published as a conference paper at ICLR 2025
1984), and we refer the reader to Appendix C for details. Then, we can consider the following map
Tht Xp 11 V/ + 11 V: GH a::= [..t,, ((te, uu:= [w(te+1 :)..................... (tetme) .)]
(2.3)
where fitn, I T. The input and the output the training data for TH are snapshots obtained from a solution trajectory with the same initial condition wo that is drawn from certain compactly supported probability measure μ.. The operator learning task for NSE using a prototype spatiotemporal operator learner, e.g., FNO3d (Li et al., 2021), is to learn this G between fixed number of Cartesian products of spatial Sobolev spaces. Specifically, Tt X2 + )' maps elements in X₁ 11 to elements in V = _1 V. In the task for this prototype "snapshots learner", the number of snapshots E z is fixed. As such, X, and Y' represent spaces of two non-overlapping temporal segments of solution snapshots. Based on these snapshots' discretized approximations from data, these snapshots learners learn an operator Go -V,, where the number of parameters in Go is independent of the spatial discretization size, yet does depend on the number of smapshots f. This dependence in the setting of this task makes not "trajectory-to-triae -otrraecctoor""
In this study, the learning aims to recast (2.3) to a vajectory-to-trajece conforming to the Hilbertian formulation of NSE using Bochner spaces. Using the fact that the weak solutions to (2.1) and (2.2) at given time interval T is in L (T; V), the operator to be learned is
operator learning problem,
G" :L,(tt, V) - L L (fe+1, (e+n,; for V == H (T²) or {u E H (T²) V = 0}. (2.4)
In what follows, we shall present how to modify any common FNO variant to become a Bochner space operator learner that can learn maps from arbitrary-sized discretization in Ixxxxxx to R. xnxnxd
(d = 1 in V-S; d = 2 in V-P). This operator can be trained by standard end-to-end supervised learning pipeline using lower-resolution data, and for inference, the newly proposed spatiotemporal rajectory-to-trrjectoo learner, f, Rt, can all be of variable sizes. In contrast, the snapshot learner FNO3d predicts the subsequent n, snapshots with the first snapshots as input, with both n, and g fixed, and only the spatial resolution n can be much higher than the input-output pairs used in the training. Hereafter we omit the d dimension if no ambiguity arises.
2.2 SPATIOTEMPORAL ADAPTATION OF FOURIER NEURAL OPERATORS
For any FNO variant, such FNO3d (Li al., 2021) or Factorized FNO (Tran al., 2023), with the following meta- meta-architecture: Ge := Q ag KL-1 o ... o d1 Ky P, where P is lifting operator, Q projection operator that does pointwise reduction in the channel dimension, of can be a pointwise universal approximator or simply chosen as nonlinearity. Kj == Ks, is the parametrized spectral convolution that uses spatiotemporal 3D FFT. During training, the operator to be learned is restricted to continuous spatial Sobolev space V in the product spaces is replaced by a finite-dimensional subspace V S ~ x with continuous embeddings {as E [ + X, and {us R g xxxn) - Y The positional encodings p == (t,, 21), E TR3xExxxx represents the underlying spatiotemporal grid, and is concatenated as before feeding to P.
subspaces X' 3 X, - Vn C X. in the sense that the
For spatiotemporal problems, all FNO variants novelly exploit a convenient architectural advantage of operator- valued NNs: the input temporal dimension É is treated as the channel dimension of an image, thus enabling channel mixing as a learnable temporal extrapolation. However, this neat trick coincidentally makes the lifting operator the biggest hurdle for Go become trajectory-to-traacttan operator learner between Bochner spaces, since the channel dimension of P must be hard-coded' and thus depends on the input pair's time steps.
In what follows, we use FNO3d as an example to present three new modifications to FNO3d to become trajectory-to-trajjee trajecccy learner ST-FNO3d, that is, to act as an operator that maps an arbitrary- input to an urbitrary-time-ttg output. We also note that this modification applies to any FNO variant that predicts spatial and temporal dimensions at the same time. These changes are so simple that the -oocooryyt--traectoory modification can serve as drop-in replacements for their snapshot learner counterparts, when the temporal input dimensions are fixed in a task. For the schematic difference, please refer to Figure 2, Figure 3, and Appendix D.2 for more comments.
Line 97 in the original FNO3d source code fourier_3d seeerrddppppp py (link provided for unaltered fork of the master hranch commit de514 in the original FNO GitHub repository).
5

---

Published as a conference paper at ICLR 2025
{ug 1
XXIXXXIN
Figure 2: The FNO3d in Li et al. (2021) is a snapshot learner. spectral convolution layer: ; pointwise nn. Conv3d that works as channel expansion/reduction: pointwise nonlinearity. The TikZ source code to produce this figure is modified from the examples in Iqbal (2018).
duxxxix
Figure 3: The catpitttmmoooll-apppttttttt FNO3d (ST-FNO3d) is now enjjectory-to-traectte learner. D: layer normalization to replace a hard-coded global normalization. Combined with channel mixing, the first spectral convolution layer serves as time-depth-wise separable (global) convolution, after which the time dimension is shrank to a fixed "latent" time dimension through iFFT'ss resampling. I the spatiotemporal spectral convolution layer as the final layer is fine-tuned after the training phase.
Modifications 3D FFT. The first change is in the spatiotemporal FFT from an analytical point of view. Using the spectral convolution K--) in FNO3d as an example, with a slight abuse of notation, denote the latent dimension (width) by d,, let the valued latent representation with de channels be continuously embedded into IT V by vs + v in each channel, and denote : [letI....; letn.), then K(-) "semi-discrete" a sense as follows
DEA / n(((t,)) x)) (s, v)) v(s, y)mm)))
(Kv) = (tv) +E.
(2.5)
where / takes values only in discrete set A, and æ E Ω. Here W E (dd+ is pointwise- xd.)
applied affine linear operator, x E C((A x Q2) x (A x 2); Rd.xd.) and m denotes an approximation to the Lebesgue measure on Ω. Evolving K into spatiotemporal spectral convolution acting Bochner space is straightforward. We change (2.5) slightly by adopting an atom-like measure d(-) in the temporal dimension, which then generalizes to the -rriable-length temporal discretization for any (t, x) E (a, b) x Q2 i.e., one can obtain arbitrarily many snapshots on the interval of interest
(K(v))(1,=) = (W (L, /1." [.
[(((((()0)) ((((1,)), (s. v), y) dm(g))d(ss.
(2.6)
During training using the discretized data, each hidden layer is a map from Rdoxdexxxn gdd.ddxxm xnxn, where d; is a "latent" dimension of time steps that is chosen S ny, The layerwise propa- gation mechanics remains the same: KΦ == F (Rg. (Fu)), where F and FFl denote the spatiotemporal Fourier transform and its inverse, respectively. The global spatiotemporal interaction characterized by the kernel is truncated in terms of modes in the frequency domain at (Timax: kimax). ) such that {(T, k) E 3 |r| < Tmax: [kyl < kmax: 1, are kept. R& is parametrized Timx Xemas Xkmax xd, -tensor. Here (T, k) denotes the coordinate in the spatiotemporal Fourier domain. Note that the integral represented by matrix product with Fx(-, (s, y)) can then be viewed residing in the continuous space the Fourier basis (3.1) can be evaluated at arbitrary point.
Modifications to the lifting operator P. The other major hurdle for FNO3d to become trajectory-to- trajectory is that a global normalization is applied with hard-coded temporal dimension using training

---

Published as a conference paper at ICLR 2025
data. In ST-FNO3d's modification, periodic padding along the temporal dimension is applied to the imput, and then a time-depth separable convolution layer I with variable time steps is used to map an arbitrary number of snapshots to a fixed number of channels d, with fixed latent time steps de- The spatiotemporal positional encodings ps- == ppP. with W, a learnable random projection with periodically padded p as input. This makes Ps match the latent fields' dimensions such that they can be concatenated. The latent fields, concatenated with ps, are then normalized by learnable layer normalization layer Ln(-), instead of global fixed pointwise normalizer for the raw data in FNO3d. Finally, in ST-FNO3d, the lifting operator becomes P = Ln(ps + I(-)) IR fxxxn + Rd. xd, x.xx,
Modifications to the projection operator Q. In FNO3d, the original out projection operator Q maps a tensor in R Xd, Xwxn to another in R. x=x", In the spatiotemporal-ddaplee FNO Q Rd Xd, xdxxnn Xxxn - R==xxxx which maps the last latent representation to match the dimension of the output uss In Q. the key for an arbitrary-length inference is that we compose another spectral convolution Kg after channel reduction. Kg can be thought of as a post-processing layer (also the only layer to be fine-tuned). It takes advantage of the FFT/FFFT's natural super-resolution capacity, especially in the temporal dimension, by gero-padding the latent temporal step dimension (d,) to given arbitrary output time steps. For the necessity of this padding, we refer the reader to Figure 7 for an illustrative example. For the V-P formulation, to impose the divergence- free condition for, S is implemented as an optional non-lcarnable Helmholtz decomposition layer in Q
NEW HYBRID PARADIGM FOR SPATIOTEMPORAL OPERATOR LEARNING
Built upon reasonably accurate approximation, the fine-tuning of ST-FNO is proposed. It is able to yield accuracy on par with traditional numerical methods on the same time horizon, and computational resources used are comparable to the evaluation of NOs. Taking advantage of the efficient ST-FNO zero-shot arbitrary-Iength temporal inference, this new approach does not need thousands of marching steps like traditional methods. Meanwhile, borrows the wisdom from traditional Galerkin methods to improve the accuracy (consistency) of the NO approach, liberating the scheme from trade-offs that the traditional methods must make to ensure stability and convergence.
3.1 A POSTERIORI ERROR ESTIMATION USING A FUNCTIONAL NORM
We shall present the proposed fine-tuning using the V-P formulation in subsequent subsections. Denote an ST-FNO evaluation at a specific im E Tr, = (tz-1.... tttn, in the output time interval u free subspace of S (t= x XS -| CV ;- H (T²), where
(m²) The construction of Q in ST-FNO renders ux
(m) W. where W is the divergence-
S = span {Re ( x - 1/2 S n/2 - 1, n E m S nt/2 1} /R,
(3.1)
where le = 2T (kz)k== and Tm == 2Tm/ 42). Then, temporal continuous approximation
uN := ux((4,) -) can be naturally defined by allowing I vary continuously on T == [4--1 Lp. T + 4p]
thanks to the spectral basis of S, where lp is the temporal periodic padding in Section 2.2. Define
residual functional R(uN)) E (T; V*): & T and for E V
R(uv)= f duuN - (uN V)UN + vruN and R(uy)(v): = (R(ux), v).
(3.2)
R(ux) measures how PDE (2.2) is violated by the finiteddimensional approximation uv, not in a localized/pointwis. fashion, but rather in global way by representing the error based on its inner product against arbitrary v. At specific time f, the functional norm of R (t, E V defined as follows is then excellent measure of the error to be candidate for a loss function in view of Theorem 3.1:
|| [(, llw suP //llll (R(un), v)--
(3.3)
Theorem 3.1 (A posteriori error bound for any fine-tuned approximations, informal version). Let the weak solution to (2.2) be u E (T; V), and deu E LT ). For u that is sufficiently regular, the dual norm of the residual is efficient to estimate the error for any uN
[ $ 1|u- uNlliA(TI)) + 10:(u uN)ll
(3.4)
P Line 205 and 209 in the original FNO3d source code fourier_sd py (link provided for an unaltered fork of the master branch commit d in the original FNO GitHub repository).
7

---

Published as a conference paper at ICLR 2025
Moreover. if and uN are "sufficiently close", then it is reliable to serve enror measure:
||ue u. (x + llu uNlliA(TTaN) & I1(u (tm) )1 + C {f
[[R(ux)(t, :) B dt.
Tm
(3.5)
where Tm == (tms (m+1l and the constants depend on the regularity of the true solution u.
3.2 FINE-TUNING USING NEGATIVE SOBOLEV NORM AS LOSS
The functional norm (3.3) is "global" because it does not have a natural G--like summation form where each summand can be evaluated in a localized neighborhood of grid points. Nevertheless, thanks to the Gelfand triple, and viewing the Fourier transform as an automorphism in the tempered distribution space (e.g., Peetre (1975) and Gel' fand & Shilov (2016, Chapter 3)), the pairing between and V as follows without getting into the intricate natures of the tempered distribution:
define
U.givv' "=" (1 [1++
(k( where é Fe for E V.
(3.6)
The spatial Sobolev space H (T(2) can be alternatively identified using norm || ||s and seminorm Is (e.g., see Ruzhansky& & Turunen (2009, Def. 3.2.2)) as follows for any E
IIf1 Σ (1+ EL+ (19+))))((())) and ffl =
Σ kE (
|x| f(
s $ O.
(3.7)
kEZ
Moreover, we have the subsequent lemma in our specific case.
Theorem 3.2 (Functional norm ">_" negative norm). E
|1/llr (
E :
-111-.
Spatially, we realize a regularized negative Sobolev norm by the Fast Fourier Transform (FFT):
-/ Σ ExEZA (a ( where 2 = (2/n2)2, and a 2 0. (3.8)
With this at hand, (3.3) and the Bochner norm in (3.4) of the residual realized to serve as the loss function in the fine-tuning to "refine" the approximation in the spectral domain
To (uy, (uNu = [[((ux)) (4m, (1mm)))l-m.a... 1)
(3.9)
Parseval identity. In the context of using an optimization algorithm to train an NN as a function approximator, it is known (e.g., Siegel et al. (2023)) that whether the integral in the loss function is accurately computed affects whether the NN can achieve the scientific computing level of ac- curacy. For example, on uniform grid, the accuracy of the weighte spatial f2-norm as the numerical quadrature is highly affected by the local smoothness of the function in considera- tion. Nevertheless, computing the integral the frequency domain yields exponentially convergent approximations (Trefethen & Weideman, 2014, Theorem 3.1) thanks to the Parseval identity.
Why funetional-type norm for the residual evaluation? We note that in "physics-informed" learning of operators, for example, in Li et al. (2024d), the PDE residual evaluated using L norm loss. In the meantime, positive Sobolev norm, which is local in terms of differential operators, is used in Li et al. (2022). To our best knowledge, the - norm has not been applied in either function learning or operator learning problems. In fact, the relation of the Gelfand triple is so simple and elegant, leading to Theorem 3.2, that the functional norm is nothing but an inverse frequency weight the frequency domain that emphasizes the learning of low frequency information. This suits especially well for the learning problem of NSE. Quite contrary to the intuition of FNO variants having the frequency truncation that results error dominating in the tighfrequency part, has been discovered in Lippe et al. (2023) that the error of operator learning for NSE is still dominant the lower end the spectrum. This has been corroborated in our study as well, see Figure 8 and Figure 6. For more detailed and mathematically enriched discussion why functional norm is not widely popular in traditional numerical methods, we refer the reader to Appendix B.
3.3 NEW TRAINING AND SPECTRAL HINE-TUNING PARADIGM
With the trajectory-tor too torracrcty learner and the eтror estimators, we propose a new training- fine- tuning paradigm. Another important motivation of this new paradigm is that, experimentally, all
8

---

Published as a conference paper at ICLR 2025
FNO variants capture the statistical properties of the 2D turbulence (Benzi et al., 1987; Boffetta & Ecke, 2012) quickly in training. Even after a single epoch, the evaluations of ST-FNO converge to a reasonably small neighborhood of the ground truth in terms of the frequency signature of the energy cascade of the flow (see e.g., Figure 8, and Appendix D.3 for more details). Then, 95% of the FLOPPs spent in training only contribute to marginal improvements. These observations motivate a more efficient paradigm than end-to-end, in the meantime not needing to initiate the expensive training of another nonlinear corrector PDE-Refiner (Lippe et al., 2023) in the postprocessing phase.
rethink
The computation of extra fields. In evaluating the loss accurately and relying on this evaluation to apply the gradient method, another key barrier is that the extra field variables in evaluating the residual for the velocity (3.2), one needs to compute 0, Instead of a common approach of using NO to rep- resent the extra field variables (Wen et al., 2022; Brandstetter et al., 2023), we opt apply traditional implicit-explicii (IMEX) numerical solver (.) with a fine time step (Wang, 2012), e.g., ((At))) for 2 2, compute these extra field variables (Line 8 of Algorithm 1), while preserving the compu- tational graph for auto-differentiaaion. We note that, in NSE simulations using traditional numerical solvers, for efficiency, the magnitude of this fine time step never realistic or attainable due time marching, Given the training data with trajectories at {fI..... tt} aiming to predict the trajectories (1++1,...,44+ the new paradigm to train and fine-tune ST-FNO is outlined in Algorithm 1.
Algorithm The new parallel-in-time spectral refining fine-tune strategy small data regime.
Input: ST-FNO Gø.e == To o Ge: time stepping scheme B.y(:: optimizer D(0, Vo())): training dataset: solution trajectories at [f,..... [ as input and at +1 4n; as output.
1: Train the ST-FNO model until the energy signature matches the energy cascade.
2: Freeze + in Ge of ST-FNO Ga.e, keep Qo trainable.
3: Cast all Modull involved and tensors torch. floate and torch. compl ex128 hereafter.
Input: Evaluation dataset: solution trajectories at [t,,.... ty] as input, output time step ?,
4: for m l,... ,
T4 1 do
5: Extract the latent fields u. (m++1) output of Ge at fm+1 and hold them fixed.
(m+1) (0) == uN + (m+))) for all m.
6: By construction of ST-FNO: uN
7: for 1, Iter do
8:
March one step with (At)T using By and approximate the duu follows: (m++)) (0)) ux (m++1) (0)) for all m.
D (m=+1))) (0) =
(m+11 (()))()) for all evaluation time steps. Compute 7 = Σm n (ux (())) Dru (m+t))))))))))) (0))
9:
10:
Apply the optimizer to update parameters in Q: 0 + D(0, vo(n²)).
11:
Forward pass only through Q update uxx"
(m+1) u (t) Qo(u(rrt)) for all m.
Output: A sequence of velocity profiles at corresponding time steps {u mmmi41
(m) 164m
Interpretations of the theoretical results. Theorem 3.1 establishes that the functional norm of the residual R(ux), without accessing u. is a good representation of the error uu in Bochner norms. Estimate (3.4) indicates that reducing the a posteriori error estimator is a necessary condition for reducing the true error. While estimate (3.5) is more delicate in that it is only reliable when uA gets "close" to u. Theorem 3.2 lays the foundation to accurately evaluate this functional norm via FFT. Nevertheless, (3.5) serves as guideline to design the "refining" procedure (lines 1 in Algorithm 1), where the error estimator moves to refine the next time step once the error in the previous one becomes less than given threshold.
4 NUMERICAL EXPERIMENTS
4.1 ILLUSTRATIVE EXAMPLE: TAYLOR- GREEN VORTEX
In this illustrative example, we examine the 2D Anylor-Green vortex model (Taylor & Green, 1937). whose amalytical solution is known and frequently employed a benchmark for evaluating traditional numerical schemes for NSE (Gottlieb & Orszag, 1977). We create a toy train dataset with 10 trajectories on 2562 grid with varying numbers of vortices per wavelength, and the test sample has

---

Published as a conference paper at ICLR 2025
an unseen number of vortices yet still can be resolved up to the Nyquist scale. For details please refer to Appendix C.2. The FNO3d used in this example is scaled down to layer.
4.2 2D ISOTROPIC TURBULENCE
This meta-example contains a series of examples of isotropic turbulence (McWilliams, 1984), featured in various work such as Kochkov et al. (2021); Li et al. (2021). The energy and the enstrophy spectra satisfy the energy law of turbulence first proposed by A. N. Kolmogorov (Kolmogorov, 1941). The data are generated using second-order time-stepping scheme that is proven in theory to preserve the inverse cascade of the energy/enstrophy spectra (Wang, 2012; Gottlieh al., 2012). We consider two cases:
the energy density in wavenumber :: = |k] is fk) ~ (A²2 + --) no drag. fixed force with low wavenumber. (II) The famous decaying turbulence discovered in McWilliams (1984), the initial power spectrum (1:) ~ (7² 2- +
(I) NSE benchmark in Li et al. (2021), == 10-3 wo
( is chosen that the decaying is slow and the enstrophy density evolves to similar energy cascade to the Kolmogorov flow used examples in Kochkov et al. (2021); Lippe et al. (2023); Sun et al. (2023).
Tables 2 and 3 report the results for example (I) and example (II), respectively. In Appendix D.2,Table 5 and Table 6 report the architectural changes and the computational efficiency comparison with roll-out and traditional solver.
5 CONCLUSION AND LIMITATIONS
We designed a new pipeline to train and fine-tune new spatiotemporal modification of FNO to get close to the true solution (not the ground truth generated by the numerical solver) of NSE in the turbulent regime. The evaluation errors in benchmark problems are up to 10 times better than the non-fine-tuned FNO variants trained by a simple end-to-end pipeline. Due the exploitation of the intricate connections with traditional spectral methods, e.g., the optimally learned parameters correspond to a Fourier-Galerkie projection (Kovachki et al., 2021) using the Fourier basis, only FNO-based neural operators benefit from this advancement. How to generalize this new pipeline to other types of operator learners in Appendix B will be worthy of exploration.
Table 1: Results for Taylor- Green vortex, the relative Table 2: Ablation study using forced turbulence. origi- errorsa the final time step, E = Wime WN.
nal example from Li et (2021). = we WN
Fvaluaaton training
After tine-tuning
Evaluntoo training
Anter fine-munng
|e|laa
[R]]-ta
Dellxa
[R]]
ST-ENO3d
[ za
[R]-1e
[
STANOM 1.94 10=1
211x100=
1.24 1
11-
10 & FT 10 H
200x1002 1.27x10-2
282x 10
2
PSARK2 (GT)
5.91x10-1
10
N/A
NA
FT
3.16x10-4
439 139x100-
Table 3: Evaluation metrics of the McWilliams isotropic turbulence example. All models are trained using 64 64 mesh and evaluated on 256 x 256 mesh. Fourier approximation space (3.1) on a 256 x 256 fine grid. r == R(ux) or R( y := HTT All error norms evaluated the final time step. H-1 appending model means the training uses the difference in the H-1 -norm as the loss function. Errors are measured for 32 trajectories in the test dataset. Fine-tuning uses 100 iterations of Adam optimizer line 10.
us uN wis wyy S stands for the
Evaluation after training
After fine-tuning
||
||]]_1, 1.m
|
[|]
-1.
FNO3d 10 ep. L train, no FT
4.32 10
6.09 10-2
N/A
N/A
FNO3d 100 ep. _ train, no FT
3.87 10
10-2
N/A
N/A
ST-FNO3d 10 ep. - train & H FT ST-FNO3d 10 ep. 1,2 train FT ST-FNO3d 10 ep, L train & 11 FT
6.54x
60-
5.71 17
2.24x 10
3.69x10-2
2.35x10-2
1.79 10
9.55 10
2.88 10
1002
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
The research of Brarda and Xi is supported by the National Science Foundation award DMS- DMS-2208412. The work of Li was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DEAC55-07NA22733 was supported by the LLNI-LDRD program under Project No. 24ERD033, Cao is supported part by the National Science Foundation award DMS-2309778.
(LLNL-CONF-87232) and
REFERENCES
Mark Ainsworth and J Tinsley Oden. unified approach posteriori error estimation using element residual methods. Numerische Mathematik, 65:23-50, 1993.
Mark Ainsworth and Tinsley Oden. posteriori error estimation finite element analysis. Computer methods in applied mechanics and engineering, (42(1-2):1-88 1997.
Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi Li, Miguel Liu-Schiaffini, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, pp. 1-9, 2024.
Constantine A Balanis. Advanced Engineering Electroomagnetics. John Wiley & Sons, 2012.
Randolph E. Bank and Bruno D Welfert. posteriori error estimates for the Stokes problem. SIAM Journal on Numerical Analysis, 28(3):591-623, 1991.
Francesca Bartolucci, Emmanuel de Bezenac, Bogdan Raonic, Roberto Molinaro, Siddhartha Mishra, and Rima Alaifari. Representation equivalent neural operators: framework for alias-free operator learning. Advances in Neural Information Processing Systems, 36, 2024.
Benzi, Patarnello, and Santangelo. On the statistical properties of two-dimensional decaying turbulence. Europhysics Letters, 317):811, 1987.
Christine Bernardi, Vivette Girault, and Yvon Maday. Mixed spectral element approximation of the Navier-Stokes equations in the stream-function and vorticity formulation. IMA journal of numerical analysis, 12(4):565-608 1992.
Guido Boffetta and Robert E Boke. Two-dimensional turbulence. Annual review of fluid mechanics, 44:427-451, 2012.
Daniele Boffi, Franco Brezzi, and Michel Fortin. Mixed finite element methods and applications, volume 44. Springer, 2013.
Andrea Bonito, Claudio Canuto, Ricardo H. Nochetto, and Andreas Veeser. Adaptive finite element methods. Acta Numerica, 33:163-485, July 2024.
Nicolas Boullé and Alex Townsend. mathematical guide to operator learning. arXiv preprint arXiv:2312.14688 2023.
Johannes Brandstetter, Rianne van den Berg, Max Welling, and Jayesh Gupta. Clifford neural layers for PDE modeling. In The Eleventh International Conference on Learning Representations, 2023.
Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932-3933 2016.
John Butcher. On the attainable order of runge-kutta methods. Mathematics of computation, 19 (91):408-417, 1965.
Edoardo Calvello, Nikola B Kovachki, Matthew Levine, and Andrew M. Stuart. Continuum attention for neural operators. arXiv preprint arXivv2406.064886 2024.

---

Published as a conference paper at ICLR 2025
Claudio Canuto, M Yousuff Hussaini, Alfio Quarteroni, and Thomas A Zang. Spectral methods: evolution to complex geometries and applications to ffuid dynamics. Springer Science & Business Media, 2007.
Lianghao Cao, Thomas O'Leary-Roseberyy Prashant Jha, Tinsley Oden, and Omar Ghattas. Residual-based error correction for neural operator accelerated ciinnte-dimensional bayesian inverse problems. Journal of Computational Physics, 486: 12104, 2023.
Shuhao Cao. Choose a Transformer: Fourier or Galerkin. In 35th Conference on Neural Information Processing Systems, volume 34. pp. 24924-24948 2021.
Carsten Carstensen and Stefan A Funken. Fully reliable localized error control the fem. SIAM Journal on Scientffic Computing, 21 2164):1465-1484, 1999.
Kathleen Champion, Bethany Lusch, Nathan Kutz, and Steven L Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 1 116(45); 2019.
Yanlai Chen and Shawn Koohy. GPT-PINN: Generative pre-trained physics-informed informed neural networks toward non-intrusive meta-learning of parametric pdes. Finite Elements in Analysis and Design, 228: 104047, 2024.
Chun-Wun Cheng, Jiahao Huang, Yi Zhang, Guang Yang, Carola-Bibiane Schönlieb, and Angelica I Aviles-Rivero. Mamba neural operator: Who wins? transformers vs. state-space models for pdes. arXiv preprint arXiv:2410.02113, 2024.
Alexandre Joel Chorin. Numerical solution of the Navier-Stokes equations. Mathematics of computa- fion, 22(104):745-762 1968.
Maarten V. de Hoop, Daniel Zhengyu Huang, Elizabeth Qian null, and Andrew M. Stuart. The cost-accuracy trade-off in operator learning with neural networks. Journal of Machine Learning, 1 (3):299-341, June 2022.
Maarten V de Hoop, Nikola B Kovachki, Nicholas H Nelsen, and Andrew M Stuart. Convergence rates for learning linear operators from noisy data. SIAM/ASA Journal on Uncertainty Quantification, 11(2):480 513,0-1130000 2023.
Gideon Dresdner, Dmitrii Kochkov, Peter Christian Norgaard, Leonardo Zepeda-Nunez, Jamie Smith, Michael Brenner, and Stephan Hoyer. Learning to correct spectral methods for simulating turbulent flows. Transactions on Machine Learning Research, 2023.
Yiheng Du, Nithin Chalapathi, and Aditi S. Krishnapriyan. Neural spectral methods: Self-supervised learning in the spectral domain. In The Twelftth International Conference Learning Representa- tions, 2024.
Alexandre Emn and Martin Vohralík. A posteriori error estimation based on potential and flux reconstruction for the heat equation. SIAM Journal on Numerical Analysis, 48(1): 223,, 2010.
Lawrence C Evans. Partial differentiall equations, volume 19. American Mathematical Society, 2022.
Vladimir Fanaskov, Alexander Rudikov, and Ivan Oscledets. Neural functional posteriori error estimates. arXiv preprint arXiv:2-402.55585, 2024.
Julian Fischer. A posteriori modeling eтror estimates for the assumption of perfect incompressibiliity in the Navier-Stokes equation. SIAM Journal on Numerical Analysis, 53(5):2178-222) 2015.
Antonio Henrique De Oliveira Fonseca, Emanuele Zappala, Josue Ortega Caro, and David Van Dijk. Continuous spatiotemporal Transformer. In International Conference on Machine Learning, pp. 7343-7365. PMLR, 2023.
L.M. Gel' fand and G.E. Shilov. Generalized Functions, Volume 2: Spaces of Fundamental and Generalized Functions. AMS Chelsea Publishing. American Mathematical Society, 2016. ISBN 9781470426590.
12

---

Published as a conference paper at ICLR 2025
Andrew Gillette, Brendan Keith, and Socratis Petrides. Learning robust marking policies for adaptive mesh refinement. SIAM Journal on Scientific Computing, 46(1)::264-A22 A289, 2024.
Vivette Girault and Pierre-Arnaud Raviart. Finite element methods for Navier-Stoke equations: theory and algorithms, volume 5. Springer Science & Business Media, 2012.
Jonathan Goodman, Thomas Hou, and Eitan Tadmor. On the stability of the unsmoothed fourier method for hyperbolic equations. Numerische Mathematik, 67(1):93-129, 1994.
Somdatta Goswami, Aniruddha Bora, Yue Yu, and George Em Kariiaaaii. Physics-informed deep neural operator networks. In Machine Learning in Modeling and Simulation: Methods and Applications, pp- 219-254. Springer. 2023.
David Gottlieb and Steven A Orszag, Numerical analysis of spectral methods: theory and applications. SIAM, 1977.
Sigal Gottlieb, Florentina Tone, Cheng Wang, Xiaoming Wang, and Djoko Wirosoetisno. Long time stability of classical efficient scheme for two-dimensional Navier-Stokes equations. SIAM Journal on Numerical Analysis, 50(1): (0k1):126-150 2012.
Daniel Greenfeld, Meirav Galun, Ronen Basri, Irad Yavneh, and Ron Kimmel. Learning to optimize multigrid pde solvers. In International Conference on Machine Learning, pp. 2415-2423, PMLR, 2019.
Albert Gu, Tri Dao, Stefano Ermon, Atri Rudra, and Christopher Ré. Hippo: Recurrent memory with optimal polynomial projections. Advances in neural information processing systems, 33: 1474-1487, 2020.
Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavvlet-based operator learning for differential equations. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, 2021.
Jayesh K Gupta and Johannes Brandstetter. Towards multi-spatiotemporalsseeo generalized PDE modeling. Transactions on Machine Learning Research, 2023.
Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 1 115(34):8505-8510 2018.
Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying. Yinpeng Dong, Songming Liu, Ze Cheng, Jian Song, and Jun Zhu. Gnot: A general neural operator Transformer for operator learning. In International Conference on Machine Learning, pp. 12556 12569. PMLR, 2023.
Zhongkai Hao, Chang Su, Songming Liu, Julius Berner, Chengyang Ying, Hang Su, Anima Anandku- mar, Jian Song, and Jun Zhu. DPOT: Auto-regressive denoising operator transformer for large-scale PDE pre-training. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 17616-17635. PMLR, 21-27 Jul 2024.
multigrid. In The Twelfth International Conference on Learning Representations, 2024. Juncai He, Xinliang Liu, and Jinchao Xu. MgNO: Efficient parameterization linear operators via
Yinnian He and Weiwei Sun. Stability and convergence of the Crank-Nicolson/ Addams-Bashton scheme for the dependdnnd 45(2):837-86) 2007.
equations. SIAM Journal on Numerical Analysis,
Jan S Hesthaven, Gianluigi Rozza, Benjamin Stamm, et al. Certified reduced basis methods for parametrized partial differential equations, volume 590. Springer, 2016
John G Heywood and Rolf Rannacher. Finite element approximation of the nonstationary Navier- Stokes problem, part ii: stability of solutions and error estimates uniform in time. SIAM journal on numerical analysis, 23(4):750-777, 1986.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin Advances in Neural Information Processing Systems, volume 33, pp.
Curran Associates, Inc., 2020.
13

---

Published as a conference paper at ICLR 2025
Thomas Y Hou. ulow-up blow-up? unified computational and analytic approach to 3d incompressible euler and Kavier-Stoker equations. Acta Numerica, 18:277-346 2009.
Jun Hu and Pengzhan Jin. A hybrid iterative method based on mionet for pdes: Theory and numerical examples. arXiv preprint arXiv:2402.07156. 2024.
Zheyuan Hu, Nazanin Ahmadi Daryakenari, Qianli Shen, Kenji Kawaguchi, and George Em Kar- niadakis. State-space models are accurate and efficient neural operators for dynamical systems. arXiv preprint M.XXv:2200.032211, 2024.
Ru Huang, Ruipeng Li, and Yuanzhe Xi. Learning optimal multigrid smoothers via neural networks. SIAM Journal on Scientific Computing, (183):119--222) 2022.
Ru Huang, Kai Chang, Huan He, Ruipeng Li, and Yuanzhe Xi. Reducing operator complexity in algebraic multigrid with machine learning approaches. arXiv preprint arXivv2307.07695, 2023.
Haris lqbal. Harisipbal&&/llommurraatttt v1.0.0.
2018.
Hans Johnston and Jian-Guo Liu. Accurate, stable and efficient Navier-Stokes solvers based on explicit treatment of the pressure term. Journal of Computational Physics, 199(1):221-255) 2004.
George Em Karniadakis, Moshe Israeli, and Steven Orszag. High-order splitting methods for the incompressible Navier-Stokes equations. Journal of computational physics, 97(2):414-444 1991.
Aly-Khan Kassam and Lloyd N Trefethen. Fourth-order time-stepping for stiff pdes. SIAM Journal on Scientific Computing, 26(4):1214-1233, 2005.
Ryan Keisler. Forecasting global weather with graph neural networks. arXiv preprint arXiv:2202.07575. 2022.
Georgios Kissas, Jacob H Scidman, Leonardo Ferreira Guilhoto, Victor M Preciado, George J Pappas. and Paris Perdikaris. Learning operators with coupled attention. Journal of Machine Learning Research, 23(215):1-63, 2022.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine bearning-acceleratee computational fluid dynamics. Proceedings of the National Academy of Sciences, 118(1))ee10178411 2021.
Andrey Nikolaevich Kolmogorov. The local structure of turbulence incompressible viscous fluid for very large reynolds. Numbers. In Dokl. Akad. Nauk SSSR, 30:301, 1941.
Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for Fourier neural operators. Journal of Muchine Leurning Research, 22(290):1-76, 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1 2023.
Heinz-Otto Kreiss and Joseph Oliger. Stability of the fourier method. SIAM Journal on Numerical Analysis, 16(3):421-45) 1979.
Hwar C Ku, Thomas D Taylor, and Richard S Hirsh. Pseudospectral methods for solution of the incompressible Navier-Stokes equations. Computers & fluids, 15(2):195-214, 1987.
Samuel Lanthaler, Roberto Molinaro, Patrik Hadorn, and Siddhartha Mishra. Nonlinear reconstruction for operator learning of PDEs with discontinuities. In The Eleventh International Conference on Learning Representations, 2023.
Peter D Lax. Functional analysis volume 55, John Wiley & Sons, 2002.
Sanjiva K Lele. Compact finite difference schemes with spectral-like resolution. Journal of computa- tional physics, 103(1):16-42, 1992.
Zhijie Li, Wenhui Peng, Zelong Yuan, and Jianchun Wang. Long- predictions of turbulence by implicit U-net enhanced Fourier neural operator. Physics of Fluids, 35(7), 2023a.
14

---

Published as a conference paper at ICLR 2025
Zijie Li, Kazem Meidani, and Amir Barati Farimani. Transformer for partial differential equations' operator learning. Transactions on Machine Learning Research, 2023b.
Zijie Li, Dule Shu, and Amir Barati Farimani. Scalable Transformer for PDE surrogate modeling. Advances in Neural Information Processing Systems, 36, 2024a.
Zijie Li, Anthony Zhou, Saurabh Patil, and Amir Barati Farimani. Cafa: Global weather forecasting with factorized attention on sphere. arXiv preprint arXiv:2405.07395, 2024b.
Zijie Li, Saurabh Patil, Francis Ogoke, Dule Shu, Wilson Zhen, Michael Schneier, John R. Buchanan, and Amir Barati Farimani. Latent neural pde solver: A reduced-order modeling framework for partial differential equations. Journal of Computational Physics,
2025.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021.
Zongyi Li, Miguel Liu-Schiaffini, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning chaotic dynamics in dissipative systems. In Advances in Neural Information Processing Systems, 2022.
Zongyi Li, Nikola Kovachki, Chris Choy, Boyi Li, Jean Kossaífi, Shourya Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Arizzadenesheli, et al. fecometry-informe neural operator for large-scale 3d pdes. Advances in Neural Information Processing Systems, 36, 2024c.
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Axizzadenesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1-27, 2024d.
Marten Lienen, Jan Hansen-Palmus. David Lüdke, and Stephan Günnemann. Generative diffusion for 3d turbulent flows. arXiv preprint 2rXiv:2306.01776, 2023.
Jacques Louis Lions and Enrico Magenes. Non-homogeneous boundary value problems and applica- fions: Vol. 1, volume 181. Springer Science & Business Media, 2012.
Phillip Lippe, Bastiaan S. Veeling, Paris Perdikaris, Richard E Turner, and Johannes Brandstet- ter. PDE-refiner: Achieving accurate long rollouts with neural PDE solvers. In Thirty-seventh Conference on Neural Information Processing Systems, 2023.
Xinliang Liu, Bo Xu, Shuhao Cao, and Lei Zhang. Mitigating spectral bias for the multiscale operator learning. Journal of Computational Physics, 506: 1 1 12944, 2024.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 333):218-22,, 2021.
Martine Marion and Temam. Nonlinear Galerkin methods: the finite elements case. Numerische Mathematik, 57:205-226, 1990.
Nick McGreivy and Ammar Hakim. Invariant preservation in machine learned pde solvers via error correction. arXiv preprint 131 2023.
James McWilliams. The emergence of isolated coherent vortices in turbulent flow. Journal of Fluid Mechanics, 146:21-43, 1984.
J Tinsley Oden, Weihan Wu, and Mark Ainsworth. An a posteriori error estimate for finite element approximations of the Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering, 1122-22) 1994.
Steven A. Orszag. Numerical simulation of incompressible flows within simple boundaries. i. Galerkin (spectral) representations. Studies in Applied Mathematics, 50(4):293-327, 1971a.
Steven A. Orszag. On the elimination of aliasing in finite-difference schemes by filtering high- wavenumber components. Journal of the Atmospheric Sciences, 1874)10074 1971b.
15

---

Published as a conference paper at ICLR 2025
Steven A Orszag, Comparison of pseudospectral and spectral approximation. Studies in Applied Mathematics, 51(3):253-259) 1972.
Anthony T Patera, Gianluigi Rozza, et al. Reduced basis approximation and posteriori error estimation for parametrized partial differential equations, 2007.
GS Patterson and Steven A Orszag. Spectral calculations of isotropic turbulence: Efficient removal of aliasing interactions. Physics of Fluids, 14(11):2538, 1971.
Jaak Peetre. On spaces of Triebel-Lizorkin type, Arkiv för Matematik, 13(1):123-130, 1975.
Kai Qi and Jian Sun. Gabor-filtered fourier neural operator for solving partial differential equations. Computers & Fluids, 274:106239, 2024.
Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadenesheli. U-NO: U-shaped neural operators. Transactions on Machine Learning Research, 2023.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Thysics-informee 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.
Rolf Rannacher. Finite element methods for the incompressible Navier-Stokes equations. Springer, 2000.
Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci, Rima Alaifari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning pdes. Advances in Neural Information Processing Systems, 36, 2024.
Sergey Repin. A Posteriori Estimates for Partial Differential Equations. Walter de Gruyter, September 2008. ISBN 10191500
Sergey I Repin. A posteriori error estimation for nonlinear variational problems by duality theory. Journal of Mathematical Sciences, 99:927-935, 2000.
Ricardo Buitrago Ruiz, Tanya Marwah, Albert Gu, and Andrej Risteski. On the benefits of memory for modeling dime-dependent pdes. arXiv preprint arXiv:2409.0331,,, 2024.
Michael Ruzhansky and Ville Turunen. Pseudo -efuddddiffeeentall operators and symmetries: background analysis and advanced topics, volume 2. Springer Science & Business Media, 2009.
Jie Shen. Efficient spectral-galerkin method i. direct solvers of second-and fourth-order equations using legendre polynomials. SIAM Journal Scientific Computing. 15(6):149--1))) 1500, 1994.
Benjamin Shih, Ahmad Peyvan, Zhongqiang Zhang, and George Em Karniadakis. Transformers as neural operators for solutions of differential equations with finite regularity. Computer Methods in Applied Mechanics and Engineering, 434: 17560, 2025.
Dule Shu, Zijie Li, and Amir Barati Farimani. A physics-informed diffusion model for high-fidelity flow field reconstruction. Journal of Computational Physics, 478:1 11972, 2023.
Jonathan W. Siegel, Qingguo Hong, Xianlin Jin, Wenrui Hao, and Jinchao Xu. Greedy training algorithms for neural networks and applications pdes. Journal of Computational Physics, 484: 112084, 2023.
Leslie N Smith and Nicholay Topin. Super-connergeenc:: Very fast training of neural networks using large learning rates. In Artificial intelligence and machine learning for multi-domain operations applications, volume 11006, pp. 369-386. SPIE, 2019.
Zhiqing Sun, Yiming Yang, and Shinjae Yoo. A neural PDE solver with temporal stencil modeling. In International Conference on Machine Learning, pp. 33135-33155. PMLR, 2023.
Eitan Tadmor. Stability analysis of finite difference, pseudospectral and mourier-galerkin approxima- tions for time-dependent problems. SIAM review, 29(4)5525-555, 1987.
16

---

Published as a conference paper at ICLR 2025
Ali Taghibakhshi, Scott MacLachlan, Luke Olson, and Matthew West. Optimization- baase algebraic multigrid coarsening using reinforcement learning, Advances neural information processing systems, 12119-12221110 2021.
Ali Taghibakhshi, Nicolas Nytko, Tareq Uz Zaman, Scott MacLachlan, Luke Olson, and Matthew West. Mg-gon: Multigrid graph neural networks for learning multilevel domain decomposition methods. In International Conference on Machine Learning, pp. 33381-33395. PMLR, 2023.
Geoffrey Ingram Taylor and Albert Edward Green. Mechanism of the production of small eddies from large ones. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences, 1588895)4499-21, 1937.
Roger Temam. Navier-Stokes equations and nonlinear functional analysis. SIAM, 1995.
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized fourier neural operators. In The Eleventh International Conference on Learning Representations, 2023.
Lloyd N Trefethen and JAC Weideman. The exponentially convergent trapezoidal rule. SIAM review, 56(3):385 458, 2014.
Kiwon Um, Robert Brand, Yun Raymond Fei, Philipp Holl, and Nils Thuerey. Solver-in-the-loo. Learning from differentiable physics to interact with iterative pde-solvers. Advances in Neural Information Processing Systems, 33:6111-6122 2020.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
Andreas Veeser and Rüdiger Verfürth. Explicit upper bounds for dual norms of residuals. SIAM journal numerical analysis, 47(3):2387-2400 2009.
Andreas Veeser and Rüdiger Verfürth. Poincaré constants for finite element stars. IMA Journal of Numerical Analysis, 32(1):30-44) 47,0004477000 2012.
Rüdiger Verfürth. A posteriori error estimators for the Stokes equations. Numerische Mathematik, 55 (3):009-325, 1989.
Rüdiger Verfürth. posteriori estimates for nonlinear problems. finite element discretizations of elliptic equations. Mathematics of Computation, 62(206):445-4777 1994.
Rüdiger Verfürth. A posteriori error estimates for finite element discretizations of the heat equation. Calcolo, 40(3):195-212, 2003.
Rüdiger Verfürth, A posteriori error estimation techniques for finite element methods. OUP Oxford, 2013.
Cheng Wang and Jian-Guo Liu. Analysis of finite difference schemes for unsteady navier-stokes equations in vorticity formulation. Numerische Mathematik, 91:543-576, 2002.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informed deeponets. Science advances, 7(40))ea 2021.
Xiaoming Wang. An efficient second order in time scheme for approximating long time statistical properties of the two dimensional Navier-Stokes equations. Numerische Mathematik, 121(4): 753-779, 2012.
E Weinan and Jian-Guo Liu. Projection method i: convergence and numerical boundary layers. SIAM journal on numerical analysis, pp. 1017-1057, 1995.
Gege Wen, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson. U- fno enhanced fourier neural operator-based deep-learning model for multiphase flow. Advances in Water Resources, 163: (44180, 2022.
Haixu Wu, Tengge Hu, Huakun Luo, Jianmin Wang, and Mingsheng Long. Solving nigh-dimensional pdes with latent spectral models. In International Conference Machine Learning, 2023.
17

---

Published as a conference paper at ICLR 2025
Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: A fast transformer solver for pdes on general geometries. In International Conference on Machine Learning, 2024a.
Hao Wu, Huiyuan Wang, Kun Wang, Weiyan Wang. Changan Ye, Yangyu Tao, Chong Chen, Xian- Sheng Hua, and Xiao Luo. Prometheus: Out-of-distribution fluid dynamics modeling with disentangled graph ODE. In Forty-first International Conference Machine Learning, 2024b.
Zipeng Xiao, Zhongkai Hao, Bokai Lin, Zhijie Deng, and Hang Su. Improved operator learning by orthogonal attention. In Forty-first International Conference on Machine Learning, 2024.
Liang Yan and Tao Zhou. An adaptive surrogate modeling based deep neural networks for large- scale bayesian inverse problems. Communications in Computational Physics, 28(5):2100-2205, June 2020.
Jiachen Yang, Tarik Dzanic, Brenden Petersen, Jun Kudo, Ketan Mittal, Vladimir Tomov, Jean- Sylvain Camier, Tuo Zhao, Hongyuan Zha, Tzanio Kolev, et al. Reinforcement learning for adaptive mesh refinement. International Conference on Artificial Intelligence and Statistics, pp. 6014. PMLR, 2023.
Hongkai Zheng, Weili Nic, Arash Vahdat, Kamyar Azizzadenesheli, and Anima Anandkumar. Fast sampling of diffusion models via operator learning. In International Conference on Machine Learning pp. 22339--24000 PMLR, 2023.
Jianwei Zheng, LiweiNo, Ni Xu, Junwei Zhu, XiaoxuLin, and Xiaoqin Zhang. Alias-free Mamba neural operator. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
18

---

Published as a conference paper at ICLR 2025
A
NOTATIONS
Table 4: Notations used in an approximate chronological order and their meaning in this work.
Notation
Meaning
Q= (0, 2x)² or (0, 1)2 the modeling domain IR2 Hilbert spaces defined on the domain St above, f E H : - R
V, H
X.Y
product spaces II H.
II
concatenation, for fj E H. IIieaJs E IsE H viscosity, the inverse of the Reynolds number, strength of diffusion
T
T := S* x where 5 homeomorphic to [o, 1) with the point being satisfies the periodical boundary condition)
(TT)
[u (2))
(T9)
H (T9 := H (T3) x H²(T²)
rot
rot प/ := (0,4) -d,4) is the rotated gradient vector
IP(T:V)
the Bochner spaces containing {u T V JT [[u[t, IlB dt < -.....ppp (o)
(T:V)
Bochner space containing (u T H suPteT [lu(t, < 12) the true wcak solution to the NSE (t, :) E V for any t
u
us
the ground truth data generated at & by the numerical solver in S C V
uN
(m)
neural operator evaluations at Im that can be embedded in S C V the inninner product on H. (u, v) := Sg uv dæ the I7*-norm of u, computed by (3.7). ||u] := ||ullo falls back to Lo-norm the //*-seminorm of u, computed by (3.7), is a norm on HTT(T)) a S b means that 3e independent of asymptotics if any such that S cb = b ++ a S and b dual space of V. contains all continuous functionals ff such that [f('] < ||el| the pairing between v E V and f V", can be identified as (f.v) if E H is continuously embedded in such that |||]| S ||u||y
(u, e) or (u,
|u/l.
hul.
V
(f. v)
4 H.
V E H. E V
compact embeddings by Poincaré inequality, H. = (T - (T?)
A B
A B Σ Eisigg2 A. for A, B e R ab T if both are viewed as column vectors
a 8@0 b
B DETAILED LITERATURE REVIEW AND MOTIVATIONS
Interplay of deep learning and PDEs. PDE solvers are function learners to represent PDE solutions using neural networks (Han et al., 2018; Raissi et al., 2019; Chen & Koohy, 2024). PDE discovery encompasses all the techniques dedicated to identifying and optimizing PDE coefficients from data (Brunton et al., 2016; Champion et al., 2019). Recently, hybrid solver approach has been explored in Greenfeld et al. (2019); Kochkov et al. (2021): Huang al. (2023); Taghibakhshi et al. (2023; 2021); Hu & Jin (2024); Huang et al. (2022). Reinforcement learning has been applied in the field of mesh generation to build more efficient solving pipelines (Yang et al., 2023; Gillette et al., 2024). For neural operators, remarkable outcomes are achieved in diverse applications, such as weather forecasting (Keisler, 2022) and turbulent fluids simulations (Shu et al., 2023; Li et al., 2023a; Lienen et al., 2023; Li et al., 2025), the methods based on neural operators have significantly influenced the progress of the interplay between PDE and deep learning. This success was a natural outcome of several improvements brought to the field, for example, fast solution evaluations, feature very appealing in many engineering applications (Zheng et al., 2023); and the ability to provide mesh-free and resolution-indepeet independent solvers in cases of irregular domains (Hao et al., 2023; Li et al., 2024c).
Neural operators. Looking at neural operator architectures, we can identify two different ap- proaches. The first one creates "basis" (or frames) through nonlinear approximators and aggregates them linearly. Developed architectures that belong to this category include Deep Operator Network (DeepONet) (Lu et al., 2021; Goswami et al., 2023; Wang et al., 2021), wherein aggregation occurs via linear combination. Similarly, Fourier Neural Operator (FNO) (Li et al., 2021), and some of its variants (Gupta et al., 2021; Rahman et al., 2023; Tran et al., 2023; Li et al., 2024c; Wen et al., 2022), aggregate different latent representations through convolution with learnable kernel in the
19

---

Published as a conference paper at ICLR 2025
frequency domain. Additional examples with linear aggregation with U-Net meta- aretitcrci can be found in Raonic et al. (2024); He et al. (2024). Bartolucci et al. (2024) further explores the error analyses for the linear aggregations through the lens of the frame(let) theory. The second ap- proach to designing neural operator architectures aims to obtain the "hasis" through linear projections of the latent representations. In this scenario, the aggregation is nonlinear, for example, using signal-dependent kernel integral. In this category, the most notable example is the scaled dot-product attention operator in Transformer (Vaswani et al., 2017), which builds an instance-dependent kernel. In the context of operator learning applications, Transformer-baas based neural operators have been studied in Cao (2021); Kissas et al. (2022); Li et al. (2023b); Hao et al. (2023); Wu et al. (2023): Li al. (2024a); Fonseca et al. (2023); Li et al. (2024b); Calvello et al. (2024); Liu et al. (2024); Xiao et al. (2024); Hao al. (2024); Wu et al. (2024a); Shih et al. (2025).
also shown in Lanthaler et al. (2023) that nonlinear aggregations outperform its linear counterparts in learning solutions with less regularity. More recently, for spatiotemporal problems, the state-space based neural operators explore an expansion to the natural operator exponential representation of the solution map, see c.g., Zheng et al. (2024); Cheng et al. (2024); Hu et al. (2024); Ruiz et al. (2024).
Numerical methods for NSE. The NSE falls into the category of stiff PDE (system) 0,u Lu (u) + f. where f is the external forcing, L and (()) are linear and nonlinear operator, respectively. Trefethen noted back in Kassam & Trefethen (2005) on the difficulty to design a time-stepping scheme as (-) and L have to be treated differently. There are a long history of numerical methods for NSE we draw inspiration from. Petrov-Galerkin methods have been developed for NSE in Boffi et al. (2013); Girault & Raviart (2012). Nonlinear Galerkin method (Marion & Temam, 1990) inspires us to prove Theorem E.9. Chorin (1968) uses a clever trick to impose the divergence free condition without constructing a divergence-conformie
finite element subspace. Shen (1994) designed various Galerkin methods in the space of orthogonal polynomials. Pioneered by Chorin, Shen, and E, projection methods are among the most popular schemes to solve NSE (see e.g., Weinan & Liu (1995) for a summary), also serves as an inspiration to add the Helmholtz layer. Bernardi et al. (1992) proposed mixed discretization for the vorticity stramffttttt formulation.
Consistency-ssabiiityy trade-offs. In view of the Lax equivalence principle ("consistency"" + "sta- bility" = "convergence" (Lax, 2002, Theorem 8)), the improvement of the stability of the method has trade-off with a method's 'ss consistency at the cost of the approximation capacity. Numerical methods for NSE is the epitome for such trade-off. For example, high-order explicit time stepping schemes offer better local truncation error estimates near boundary layers of the flow (Lele, 1992), yet the lack of stability is more severe and needs higher-order temporal smoothing. On the other hand, implicit schemes can be unconditionally stable for stiff or highly transient NSE with relatively large time steps (1). The stability in implicit schemes becomes much less stringent on the time step, as no CFL condition is required. However, the solution at the next time step requires solving a linear or nonlinear system. Thus, computationally implicit schemes are usually an order of magnitude more expensive compared to explicit schemes. The explcit-explicit (IMEX) stable solver chosen in the fine-tuning postprocessing is from Wang & Liu (2002); Wang (2012).
De-aliasing filter sacrifices consistency for stability. One famous example of the consistency- stability trade-off is the 3/2-rule (also known 2/3-dealiasing filter) for the nonlinear convective term (Orszag, 1971b; Patterson & Orszag, 1971; Hou, 2009; Gottlieb & Orszag, 1977) for pseudo- spectral methods (Orszag, 1971a; 1972). The highest 1 1/3 modes, the inclusion of which contributes to better approximation capacity, are filtered out to ensure long-term stability. Compromises such as the CFL condition and the 3/2-rule must be made for traditional numerical schemes be stable, keeping the balance between stability and accuracy. These constraints apply to traditional numerical methods because any solver has to march consecutive multitude of time steps, which makes the error propagation operator's norm a product of many. Numerical results suggest that for pseudo- spectral spatial discretization with no higher-order Fourier smoothing temporally, the dealiasing filter is indispensable (Tadmor, 1987), as the time marching may experience numerical instability (Kreiss& & Oliger, 1979; Goodman et al., 1994) without it. This is due to the nonlinear interaction in the convective term, which is caused by the amplification of nigh-frequency "aliasing" errors when the underlying solution lacks sufficient smoothness. In this study, the hybrid approach we adopted combines the strengths of NOs and traditional numerical solvers. There is time marching
20

---

Published as a conference paper at ICLR 2025
consecutively for multitude of time steps, which renders the method free of the stability constraints such as the CFL condition and the 3/2-rule.
Why and why not functional-type posteriori error estimation? For the sake of presentation, make some handy assumptions: (1) A V - V denotes a linear operator on the Gelfand triple; (2) uss is obtained through Galerkin methods (Evans, 2022), which is projection onto approximation subspace of S C V; and (3) an evaluation of uN can be continuously embedded into S. Then, we have the following enlightening Pythagorean-type identity for the true solution satisfying Au- =
dintee-iimenional
||u - usl2 2 + llus - uvllla = ||u uvlla since (u us) La (us - uN) E V.
(B.1)
v) = 0, for any E S. This bilinear form (::)a induces a (semi)norm) || lla inheriting the topology from V. Given this orthogonality, minimizing the difference between uN and us becomes fruitless if ||u- uslla laa is relatively big in the first place, and unnecessary computational resources may have been spent to get closer us. Rather, (B.1) indicates that it is more efficient if one design method to reduce the error of | - uxilla directly, while circumventing the fact that u is not accessible.
Through the Riesz representation theorem, us is solved through (u- (u- v),, ==
Speaking of the posteriori enror estimation in traditional PDE discretization, part of the goal is to help the adaptive mesh refinement to get a better local basis. In this regard, the global error functional in negative Sobolev spaces must be approximated using localized 1.2 residuals to indicate where the mesh needs to be refined. There are various compromises for this H -to-to- representation to happen that renders the estimate inaccurate, such as discrete Poincaré constant (Veeser & Verfürth, 2012) or inverse inequalities (Carstensen & Funken, 1999; Veeser & Verfürth, 2009), see also Verfürth (2013, $1.6.2). Functional-ty -tppe a posteriori error estimation (Repin, 2008) consider the error as functional, which is equivalent to error to bilinear form-associated norm as follows
R(us) E v", and (R(us), v) (u us v).: where the weak form of the PDE is
(u. v). f(v) Vu V and (us, v)a = f(u) Ve E
The common approach is using the help from extra "flux" "stress" dual variables (Repin, 2008, $6.4), see also Repin (2000) for how to get an accurate representation of the error functional in L.2. For example, for the Stokes flow, which can be viewed that 0, u = 0 and no nonlinear convection in (2.2), Re ~ O1), Repin (2008, $6.2) estimates error of an H (d? mixed discretization as follows
steady limit of viscous fluid, such
vl|V(u us)ll S ||o + qd - vusl CrlV o a+ +fll.
where Cp is the Poincaré constant of the compact embedding and (ar, a reconstruction field pair. However, the drawback of this approach is that an expensive global minimization problem needs to be solved, e.g., for (ar, q) above. Another main reason to introduce extra field variables that, for finite element methods, the basis functions are local and do not have globally continuous derivative, in that Aus yields singular distributions, whose proper norm is -norm and canmot be evaluated by summing up element- wise L no-moss.. In computation, it has to be replace by V where o E H (div), and is constructed to be closer to the true solution's 's gradient Vu than Vus. Meanwhile, for systems like NSE, the error estimation in Galerkin methods for NSE is further complicated by its saddle point nature from the divergence-free constraint, in that consistency has be tweaked ensure the Ladyzhenskaya- Baeaaaaaa-Bauuua--aaeaee 2012, Chapter III $1 Sec 1.2). Recently, Fanaskov et al. (2024) considers NN as function learner to represent solutions of linear convection-diffsiion equations, yet still falls into the traditional error estimation framework in that an extra flux variable is learned by NN and the error if of typical accuracy of NN function learners (110-3). In contrast, in our study, the need of extra "stress" or "flux" variables to build the residual functional is circumvented as well.
stability condition (Girault & Raviart
Error correction for Bayesian inverse problems. The a posteriori error-correction approaches through sampling for the output of surrogate NNs in Bayesian inverse problems (Yan & Zhou, 2020). More recently, in Cao et al. (2023), the error equation in an inverse problem is approximated by solving for a Galerkin projection of a linearized eтror equation, while in this paper, we leverage the Gelfand triple directly to compute a nonlinear Galerkin projection directly by minimizing the spectral norm.
21

---

Published as a conference paper at ICLR 2025
Sobolev norms in operator learning In the context of function learning problems, Du et al. (2024) applied an L -'pspectral loss for Physics-in informed Neural Networks (PINN) (Raissi et al.. 2019)
and found superior results over mesh- weighted spatial losses. Du et al. (2024) realized the PDE residual as a functional on I,2 and exploited the Parseval identity on this space, which computes su R, v) /1/ We note that this is not the natural space to measure the PDE residual according to the Galerkin formulation of NSE. Using the Poisson equation AA = f as a simpler example, if II', then the residual f Au E H H A is a bounded operator from Hf to H- which guarantees the stability estimate ||u][] S cl while - is unbounded from L to L², The dual space of H' (T²) is the natural space to consider the residual functional. The reason is that the solution to NSE belongs to H' H' (T²) at given snapshot. The dual space of H (T²) happens to be H-1 using the Gelfand triple.
this paper, we compute H (R, u)/ /////==| However, this H---norm a non-localizable funetional norm the spatial domain, and traditional wisdom (FEM or FVM) has to circumvent the direct evaluation of it by localizing the residual, which introduced inaccurate compromises in the residual-baas based error estimation.
A posteriori error estimation for flow problems There is vast amount of literature for the posteriori error estimation for the viscous flow (1/Re > 0) problems under the Hilbertian framework for traditional numerical methods. Our methods draw inspiration from these pioneers and try to address the drawbacks. The most popular a posteriori error estimator for stationary Stokes problem from Verfürth (1989), and it is of residual-type by computing mesh-weighted Lo-norm elementwise, and the singular distribution Aug is represented by the magnitude of flux jumps across the facets in this mesh. In Bank & Welfert (1991), more accurate a posteriori estimation technique is invented for Stokes problem in which local problem is solved to represent the residual functional on a collection of neighboring elements. I,2 residual-type error estimation for stationary NSE is considered in Oden et al. (1994). To our best knowledge, no functional-type a posteriori error estimation has been applied to solve the transient NSE, due its in-efficiency for traditional finite element or finite volume methods. As at every time step. a global nonlinear problem has to be solved if one ought to evaluate the functional accurately using finite element local basis functions, whose computational cost is an order of magnitude higher than implicit Euler methods.
Hybrid methods for turbulent flow predictions. There are quite a few approaches that combine NN-based learners with traditional time-marchin solvers. Notable examples include: solution from Direct Numerical Simulation (DNS) corrected using an NN (Um et al., 2020); interpolation to achieve spatial super-resolution using solutions produced by traditional numerical solvers marching on coarser grids (Kochkov et al., 2021). Sun al. (2023) uses orthogonal polynomial features (HiPPO from Gu et al. (2020)) to achieve temporal super-resolution in addition to the spatial one. McGreivy & Hakim (2023) explores how to incorporate energy-preservin schemes such as upwinding into the ML solvers through traditional numerical methods. Dresdner al. (2023) explores the possibility of using neural networks to correct the result of the traditional time- marching solvers. Qi & Sun (2024) identifies important frequency bands in FNO parametrization and improves FNO's rolling-out performance in NSE benchmarks. Wu et al. (2024b) uses graph-based neural ODE to improve the out-of-distribution performance in fluid dynamics tasks.
C DATA GENERATION
C.I SARTICITY-STRAAMUUUN
FORMULATION
Here we derive the the standard NSE in 3D, u's z-component is O, and u = =u(t, x., y) has only planar dependence
formulation for the convenience of the readers. Consider
{
Ou
(u ))u Vp Vv + vsu
-
at
(C.1)
u 0
Taking V x (.) on both sides, and we assume that the solution is sufficiently regular that the spatial and temporal derivative can interchange, as well as curl and the Laplacian can interchange, one gets the following equation by letting = V x a
aw Ht
2
(v + vAw + Vxf. p
V
X
W Vw w Vu -V V
(C.2)
2
22

---

Published as a conference paper at ICLR 2025
300
450
600
06
Velocity Magnitude
45
04
03
Figure 4: Ground truth streamlines for Taylor-Green vortex example.
Next, one can have w Vu O. thus the equation becomes
dw Vw w Vu +AW+
(C.3)
Ōt
We need introduce streamfunction v recover the velocity V (0, v)). The main simplification comes from the assumption of dependence on x- and y-variable only. Thus, in 2D this becomes rot π 2 rotation of the 2D gradient Vч)). Therefore,
u (u,, W2 = 41 = 0, u2 = -2
and the vector vorticity be equivalently represented by scalar vorticity
w V x u (0, 0, curl(u,, u2)) = (0, 0,42 0,u1) =: (0, 0,w)
The equation becomes nonlinear system for and
Onw + rot
Vw - vAw Vxf,
(C.4)
w+ ト = O.
(C.5)
C.2 MAYLOR-GREEN VORTEX
Taylor-Green vortex serves as one of the most well-known benchmarks for NSE numerical methods, its flow type experiences from laminar, transitional, and finally evolving into turbulent regime. We consider the trajectory before the breakdown phase. We opt to use doubly periodic solution on [0, 2nr²² such that the inflow/outlow occurs on the "boundary" (see Figure 4). The exact solution is given by:
(1, x. = 24t
(ii( cos(Ky) (-cos( sin(ky)) ) or 2x - cos(((z)) sin sin(ray) sin((ar) cos(sey)
where v/ 1/Re chosen as 10 For sample trajectory with x 1, please refer to Figure We apply the pseudo-spectral method with RK2 time stepping for the convection term and Crank- Nicolson for the diffusion term. The dataset has trajectories with ranging from 11, among which K 11 is chosen as the test trajectory.
C.3 2D ISOTROPIC TURBULENCE
The example featured the original FNO paper (Li et al., 2021) can be viewed as a special case of the isotropic turbulence with strongly regularized initial condition, no drag, and a force with low
23

---

Published as a conference paper at ICLR 2025
wavenumber. This is example (I). The example (II) has Seakly-regularian initial condition, small drag, and no external force. The initial condition spectra are chosen such that, after the warmup time, the spectra the comparable with the Kolmogorov flow featured in Kochkov et al. (2021). The similarity is in the sense of the energy cascade formed in the Fourier domain (see Figure and 9). The training data used in both (I) and (II) generated using pseudo-spectral method with de-aliasing filter using an RK4 marching for the convection term and implicit Euler for the diffusion term. The Reynolds number used this example 1000. The time steps for both are At 10 grid, then downsampled to 642 grid. The warmup time is 4.5. The öt for the ST-FNO training and evaluations are At The testing data are generated on 512 grid with a 10- time step then downsampled to 256 There are 1024 trajectories in the training data. The input of model has snapshots from snapshots from / = 4.5 + 10 to 1 = 4.5 + 19 ot. There are 32 trajectories in the testing data. Each of the testing data's s trajectory runs from the same time interval, and has 40 snapshots in the input and 40 snapshots in the output. The initial energy distribution for the training and testing data trajectories are draw from the same energy distribution, see (D.1). For additional experiments using out-of-distribution data with different initial energy distribution and evaluation using trajectories with different Reynolds number than the training data, We refer the reader to Section D.4.
256
4.5 to / =
The output ground truth the subsequent
2.50-4
0.01
00.1
-2. 2.50m4
-1e-6
(b)
(c)
Contours plots of pointwise values of residuals for Example (I). (a): the residual of the ground truth; (b): residual of ST-FNO, where the time derivative in the residual is using the ground truth's; s; (c): the residual after fine-tuning for 10 ADAM iterations where the time derivative is computed using extra- extra-fine-stey fine-step numerical solver.
Figure 5:
10
10
20
30
30
20
0
20
40
20
40
60
(a)
(b)
Figure 6: Pointwise values of residual of ST-FNO in the frequency domain in In scale of Example (I). (a): the residual of ST-FNO before correction, time derivative computed using the ground truth: (b): the residual of ST-FNO after fine-tuning for 10 ADAM iterations where the time derivative is computed using an extra-fine-step numerical solver.
D EXPERIMENTS
D. TRAINING AND EVALUATION
The training torch. optim. OneCycleLR (Smith & Topin, 2019) learning rate strategy with 20% warm-up phase. AdamW W the optimizer with no extra regularization. The learning rate starts and ends with 10-3,, Irmax- The I 10-3 Despite using the negative Sobolev norm is quite
24

---

Published as a conference paper at ICLR 2025
Original on Fnv
1.3
Fi Fnn
-93
-10
60
02 0.4 0.6 06 10 12 1.4
Figure 7: A simple illustration of necessity of padding using a ID non-periodic data. FNO3d cannot handle temporal non-periodic data. The free temporal super-resolution by FFT/iFFT has Runge-like phenomena near the end points without padding. Padding in the ST-FNO3d modification resolves this problem. Since data themselves are intrinsically lower-dimensional, representations, but rather apply the padding only in the last spectral convolution layer to get the output.
opt not to pad the latent
efficient fine-tuning tun-tuning, using training not be the most efficient minimizing the L -o-m due to the optimization being non-convex, we observe that all norms are equivalent "bad" due to nonlinearity of NOs. The result demonstrated obtained from fixing the random number generator seed to be 1127825. All models are trained on a single RTX A4500 or RTX A6000, the inference is done with power limiter on the GPU ensure fairness. The codes to replicate the experiments are open-source and publicly available at https : http//ggttttttteee/ee
D.2 MODELS
Helmholtz layer for the velocity-pressun formulation. For the V-P formulation in a simplified connected convex or periodic domain, is known that exact divergence- free subspace W V for velocity means that the pressure field is not needed. The reason is that in the weak formulation, the pressure is Lagrange multiplier to impose the divergence- Chapter III 1 Section 1.1). Inspired by the postprocessing to eliminate Vp (Ku et al., 1987, eq. (15)) together with the neural Clifford layers (Brandstetter et al., 2023), we add Helmholtz layer S after each application of σj (W,+ + Kÿ)- S performs a discrete Helmholtz decomposition (Girault & Raviart, 2012, Chapter 1 $1 Section 3.1) in the frequency domain to project the latent fields to be solcnoidal.
condition (Girault Raviart, 2012,
Difference with Fourier Neural Operator 3d In view of (2.5), the FFTs in FNO3d transforms are continuous integrals in the spatial dimensions yet discrete sum in the temporal dimension due the channel treatment. In the original FNO3d architecture, despite that explicit restriction imposed on the output time steps, the dam cannot be trivially changed as the data are prepared by applying a pointwise Gaussian normalizer that depends on dout- Denote the lifting operator (channel expansion) in FNO3d by P {a,} + {pr} + IR Xdu xnxn Before the application of / in FNO3d, as is artificially repeated dout times, then P is in a space of linear operators ~ +3))xd This observation makes it independent of the spatial discretization size n n by dependent on d In FNO3d, the channel reduction is a pointwise nonlinear universal approximator, yet in ST-FNO3d, this is linear, which eventually facilitate the convexity of the fine-tuning optimization problem.
D.3 FINE-TUNING
In implementations, choose Itermax 100 and the ADAM optimizer with learning rate 10-! in Algorithm I, i.e., an ADAM is simply run for 100 iterations to update parameters that count only a fraction of spectral layer since the number of parametens accounts for only single channel.
25

---

Published as a conference paper at ICLR 2025
Table 5: The detailed comparison of FNO3d and its spatiotemporal modification. Layer: of spectral convolution layers; ST-FNO has an extra layer of spectral convolution in single channel with skip connection. Modes: (Tmax: mmax). Pre-norm: whether a pointwise Gaussian normalizer is applied for the input data. Eval FLOPs: Giga FLOPs for one evaluation instance. FT: finetuning GFLOPs profiles per ADAM iterations for 256 grid for single instance. A torch. cola type parameter entry counts as two parameters.
Architectures
GFLOPs
layers
channcl/width
modes
activation
pre-norm
Eval
FT
params
FNO3d Example (1)
(5,8)
GELU
13.3
N/A
9.03m
ST-FNO3d Example (1)
20
(5,8)
GELU
27.0
3.28
9.02m
FNO3d Example (II)
10
(5,32)
GELU
Y
28.7
N/A
16.38m
ST-FNO3d Example (II) 4 10
10
(5.32)
GELU
42.5
3.4
16.42m
Table 6: The FLOPs/runtime comparison. FNO3d with a roll-out prediction (10 steps to predict
next 10), ST-FNO3d (spatiotemporal prediction), and a traditional numerical solver's time marching (IMEX 4-th order Runge-Kurt Kutta). Spatial grid size 256 x 256, prediction time interval from = 4.5 to t == 5.5 (stable energy cascade regime), numerical solver At 10 roll-out St 104t, ST-FNO3d öt = 25A Note that ST-FNO3d gets all 40 steps in a single evaluation, the numerical solver has to march 1000 steps, and the roll-out prediction inference is supposed march 10 times yet becomes de-correlated around I z 4.8 after the third step. Runtimes are reported for single instance averaging 20 runs with a 100w GPU power limiter on single RTX 4500.
Architectures
GFLOPs
Runtime (x (x10-2s s)
layers
chunnel
modes
Eval/Step
FT
All
Single eval
All
# parums
FNO3d roll-out Exumple 2
4
10
(5, (2,32)
28.7
N/A
2.9M
4.111.
45.6 4: 6.7
16.38m
ST-FNO3d Example 2 IMEX Runge-Kutta
4
10
(5,32)
1.06
2.51
337.6
22.7 = 4.0
162 -: 4.1
16.42m
TO N/A
N/A
N/A
3.98
N/A
4M
0.449 0.3
450 =
10
To train to stop? Investigating the frequency domain signatures In Stuart's paper the convergence of linear operator learning (de Hoop et al., 2023), the evaluation error scales with number of samples, in the few-data regime, the overfitting comes fast. Part of the reason is that the operator learning problems are in the "small" data regime.
V x u. and Fw where F is applied only in space, then we can define enstrophy spectra as follows:
Denote w(t, x)
ktt =
E
k0tt,))13 and Σ E(t,k):= = | k
B Sk+5k
In Example (II), the enstrophy spectrum decays as O(k-3), which slightly faster than kinetic energy. This is usually referred to as the direct cascade, as opposed to the inverse cascade (kk-5/3 discovered by Kolmogorov. For ST-FNO, it learns the frequency signature of the data after single epoch. See Figure and Figure 9.
D.4 ADDITIONAL EXPERIMENTS
In this subsection, we perform the following additional experiments regarding the possibility of fine-tuning ST-FNO3d for -iitr---distrbbution"" data. All ST-FNO3d models trained for 15 epochs (5 more than the model in Table 3) with the data generated using the McWilliam''s isotropic turbulence data from Example (II) with 10-3, The initial energy distribution satisfies:
[x:)) ~ - (7 (k/ko)))-
(D.1)
The parameter ho in Example (II) roughly means that the initial energy is concentrated at [k| Ps kg- In Table 3 Ko is chosen 4 according to McWilliams (1984). For all examples, the same "bur-iin"
26

---

Published as a conference paper at ICLR 2025
10a
10
100
10a
Enstrophy density
Enstrophy density
10-
10 1
10-2
10-2
10-3
10-3
10-4
10-4
(kFNo(k)
(kF
10-5
10-5
(KT(k)
ECT(K)
10
10 6
O(k-3) 3)
O(k-3
100
10
10"
101
k = ||k]
=
|
(a)
(b)
Figure 8: Enstrophy spectrum density comparison for Example (II) to illustrate the "convergence" of the ST-FNO training. The ST-FNO evaluation is on 256 x 256 grid for a fixed randomly chosen sample, and the training is a 64 64 grid starting from 10 different seeds at different time steps. The error bars are plotted with +/- 10 times the standard deviation from the mean to boost the visibility of the convergence. (a): the comparison after epoch, the average relative L difference with the ground truth 2.150 x 10-1 (b): the comparison after 10 epochs, the average relative I difference with the ground truth is 8.051 10-
time is used 4.5), the training output time steps 10 and evaluation output time steps 40. Please refer to Figure 10 for sample trajectories and Figure 11 for the energy cascade after the "burn-in" time period.
(Ex1) The initial energy is concentrated at ko 2. A sample trajectory be found in Figure 10(a).
(Ex2) Same with above, kn = A sample trajectory can be found in Figure 10(c).
(Ex3) ko = 4 v == 1/Re- Re = 1/5000. In order to resolve the small viscosity, the mesh size becomes
512 512.
(Ex4) ko 8, v L/Re- 1/10000, mesh size is 1024 1024 for the same reason as above.
The evaluation and fine-tuning result be found in Table 7. We found that the paradigm works reasonably well for first three example, yet fails to converge for the last example. For an example plot of the prediction comparison please refer Figure 12. The correlation graph between the ST-FNO3d's predictions and the ground truth can be found in Figure 13, in which the FNO(2+IX Iyd adapted from the original FNO2d code. This roll-out-based model predicts the next snapshot based on 10 given snapshots as input, and the predicted snapshot is concatenated to the current input as the latest snapshot form the input for the next prediction.
Table 7:
Evaluation metrics the McWilliams isotropic turbulence example for "out-of -iiiii-o-ddisiiiuttionnn data. For details please refer Section D.4 ws wy. R(wy)) ) H- - All error norms are evaluated at the final time step. The training is under the 1.--ormm and the fine-tuning is under the -l-ooomm. Symbol s_" means the errors the same as the entry on the row above.
In-distribution evaluation
Out-of-distribttion evaluation
After fine-tuning
Jellza
||])
||E||a
|| -Lm
|| 12
[[]]]_.
(Exl) kg =
5.39 10 1.66
8.52x 10
4.25 2
2.73x 10
5.05 3
(Ex2)
1.46x10
2.31x10
2.06x10-2
5.85x
(Ex3) Re 5-10"
2.94x10-1
2.73x11-2
2.50x10-2
2.96x10-4
(Ex4) Re 10a
4.92x
1.76x10
2.355x10
x
27

---

Published as a conference paper at ICLR 2025
10a
109
109
10
10-2
10-4
10-1
ESPNO
ESENO
ESENO
EGT
EGT
EGT
10
O(K
10-4
O(k 2)
10-6
O(k 3)
100
10'
10P
10
10P
10
(a)
(b)
(c)
109
10
10
10 2
10-4
10
ESPNO
ESPNO
ESPNO
ECT
ECT
ECT
10-6
10-6
10-6
O(k-3)
3(k
((k-3)
10a
10
10P
10
10a 101
(d)
(e)
(f)
Figure 9: Enstrophy spectrum density comparison for Example (II) for two randomly selected
trajectories for given trained ST-FNO. Training and evaluation setups are identical those in Figure 8. (a), (b), (c): the comparison of the spectra of the ST-FNO evaluation and the ground truth after 10 epochs of training at le+1. le+7, and l0+10 of the first trajectory; (d), (e), (f): the comparison for the second trajectory at the same time steps.
E ASSUMPTIONS AND PROOFS
Assumption E.1 (Assumptions for Theorems 3.2, E.9). The The following notions and assumptions are adopted throughout the proof of the three theorems involved, all of which are common in the literature for NSE. While some of them can be proved, we opt for list them here to facilitate the presentation.
(Et) The compact embeddings for Gelfand triple V E H @ V hold, where V H (T²),, and (T²).
2F2) The time interval in consideration is fixed in that we consider the "refining" problem for fixed time interval but with variable time steps.
(Ea) The initial condition uo E V, and L3(T.V")
Lemma E.2 (Skew-symmetry of the trilinear term). For 2, u. U E H²(T²), where T denotes
Ω = (0, 1)2 equipped with component-wise PBC. Denote
u2 = ((x V)u,
((x V)u) dæ. 1
V = 0, then
c(2, u, v) c(z. u.
(E.1)
and specifically v) v)
Proof. We note that this result is common in NSE literature, see e.g., Temam (1995, Part I Sec 2.3), for homogeneous boundary whole space. In this case, it suffices verify that the boundary term
28

---

Published as a conference paper at ICLR 2025
(
(b)
Figure 10: Examples trajectories of the evaluation datasets for additional experiments for Example
(II). (a): ku 2; (b): ko 4. (c): ko
8.
10"
10P
10F
Enenggy spectra density
Emergg spectra density
10
Eivergy spectrra density
10
10
10-2
10:2
10
10
10
10:4
10:
10-5
E(K)
10-4
E(K)
10-3
E(K)
5/3)
Okk-51 5/3)
Ork-514 5/3)
8
10-7
10-6
11
301
10
301
10
k [|
k
||k]
k
[|k]
Figure 11: Energy spectrum densities datasets for additional experiments for Example (II). (a):
ku 2: (b): Ko = 4. (c): ko 8.
vanishes. Here we still provide short argument for the convenience of readers. First by the product rule. and vector calculus identity of V (v u), e.g., Balanis (2012, Appendix II.3.2).
((v x) ((= V)u) v + ((z V)v) u + vV- z.
(E.2)
By Gauss divergence theorem, then obtain:
[ V. (()))))) dæ
(v u)= ds
Σ
/ (v-u)z n, ds,
00
e,con,isisa
with n outer normal with respect to a2, and n, denotes that with respect to edge e,. Now on é1 : f- 0} {y (0, 1)). n, () T w/e, w/e, PBC where eq == for 1) x { (0, and w (z, u, v). The integrals on e and é2 cancel with one another since n2 - Furthermore, if V 0, integrating (F.2) on Ω yields the desired result.
n
29

---

Published as a conference paper at ICLR 2025
(c)
Figure 12: The ground truth versus the ST-FNO3d evaluation/fine-tuiing results the Re 5000 additional experiment for Example (II). (a): the ground truth at 5 (prediction step 20) generated by IMEX RK4; (b): ST-FNO3d prediction. (c): ST-FNO3d fine-tuned prediction.
Correlation with high-res GT
0.98
0.92
FNO (2+1)-D rollout
ST-FNO3D
0.90
ST-FNO3D (finetuned)
ST-FNO8D (super-res)
0.88
ST-FNO8D (super-res, finetuned)
0
10
20
Prediction steps
Figure 13:
Comparison of the correlation with the ground truth (on 1024 1024 grid) generated by IMEX RK4 solver.
Lemma E.3 (Poincaré inequality). V/RR E H//R is compact and the the following Poincaré inequality holds with constant for any V/R
|
(E.3)
Proof. A common textbook proof exploits the equivalence of sequential compactness and compact- ness, once the norm topology is introduced, we prove this in very intuitive way under our special setting: by definition (3.7), if E V/B then 6(0,0) 0) =
Σ
10(k)) <
Σ kE2
|k))))))) = lu/
((2)
E.1 PROOF THEOREM 3.1
To prove Theorem 3.1, we need Lemmas E.4, E.5, and E.6, all of which can be found in classical textbook on analysis and approximations NSE, such as Temam (Temam, 1995), or Girault & Raviart (Girault & Raviart, 2012) for homogeneous Dirichlet boundary condition for the velocity field or slip boundary condition for the (pseudo-)stee Jstress tensor. We shall complement some proofs with simple adaptations to PBC whenever necessary for the convenience of readers.
30

---

Published as a conference paper at ICLR 2025
Before presenting any lemmas, need the following definitions, weighted H' -norm on V H (T²) for a > 0 on fl is defined as
|,,|a,11 == (1w=112> + |o and 1 == ||a
For f E V, define a weighted dual norm as follows:
sup (f. n)
If
(E.4)
vEV ||v|la,1
Lemma E.4 (Contuinity and embedding results for the convection term). For trilinear convection term
(uu V)D, w) - |||||>> |0/1|///r.. and ((u V)v, w) Σ |4|| |0| |l|
(E.5)
Proof. The first result obtained by the Hölder inequality, and the second holds with suite Sobolev embedding results (Temam, 1995, Part I, $2.3, Lemma 2.1).
L
Lemma E.5 (Energy stability of NSE in a bounded domain).
u i (T,)) I (uoll 2 + |If ((A((T))))))
(E.6)
Proof. (E.17) is classical result, see e.g., Temam (1995, Section 3.1).
Lemma E.6 (Fréchet derivative the convection term). Given v, u E V, the linearization of the difference of the convection term reads
v. Vo = u Vu + D(u)(v u) + where D(u) = (5 V)u (u V)
(E.7)
is the Fréchet derivative and ||||| $ w ||u u/11.
Proof. This result is normally used without proof in linearizing NSE, see e.g., Girault & Raviart (2012, Chapter 4 eq. (6.5.2)). Define F(u, G) = T G, then expanding F at v and Vu yields the desired result.
Theorem 2.1 Part (Functional-type a posterior error estimate is efficient (rigorous version)). Consider the Gelfand triple V E H E V', and the weak solution be sufficiently regular, then the dual norm of the residual is efficient to estimate the error:
(T; (T;H) (2.2)
[R(UN)) S [10,(u ux)ll(A(TT + llu - uAIliACTT
(E.8)
The constant depend on the regularity of the true solution u.
Proof. The proof of the lower bound (efficiency) (E.8) follows the skeleton laid out Verfürth (2003, Lemma 4.1) for diffusion and temporal derivative terms, while not needing extend approximation using an affine linear extension operator as Verfürth (2003) does. The treatment of the nonlinear convection term follows Fischer (2015, Lemma for the steady-state NSE (no temporal derivatives), barring the technicality of the divergence-fres condition. Consider at L E T, for any test function H (T²) and V 0, integrating against R(ux) in (3.2) yields
(R(ux), v) =(f d.ur (ux . V)uN + vAuv, v) =(0, (u uv)) (ND(u -u)) + ((u V)u - (ux V)uN,
Integrating by parts for the diffusion term, using the same argument for PBC as in Lemma E.2, and inserting (u V VuuN yield
(R(ww),)) (a,(u u) (vV(u-ux), uN), Vv)
(1)
(I)
(E.9)
+ ((u v)(u - u)). + (((u (((u-u)) V) (),,,)))
(II)
(IV)
31

---

Published as a conference paper at ICLR 2025
Applying the definition of weighted dual norm (E.4) (I) we have
(I) x 1 - (t,) 1/2,-1 1
For (II) we simply have
III <
For (III,, applying Lemma E.4 and the stability estimate in Lemma E.5 for u.
(II) & C/ 1,1 || uwllua/,, || S Ci mA) | un wwwlaaa,,lolll |
For (IV), applying Lemma E.4 and the stability estimate in Lemma E.5 for uwy
(IV) S C|lu uw /ww a uv 1 / C2( uo)l(u uolllaaa,llllllaa
Applying the Cauchy-Schward inequality, has
(R(ux), w)12 = C (10. = (10:(u-uvlel) -ll2 12,-1 + llullluuee uvlliu221) |11./20/2,1
Finally, equipping V with the weighted norm (E.4), and by the definition of L L3(T:V'), we have
|R((u))((((TIV) = 1 Tee |
(R(ux), dt
S C L (110((u (w)((u-ww)(t))ll (1, l/eu uxlBra,)) dt,
where C' C(ug, 14, T).
Theorem 2.1 Part II (Functional-type posterior error estimate is reliable (rigorous version)).
Consider the Gelfand triple V @ H E V", and the weak solution is u E n to (2.2). We assume that (2.1.1) ux and u are sufficiently close in the sense that: there exists E [0,, C) a fixed C E ER+, for J defined in Lemma E.6
[J ux) S >llV(u ux)ll
(E.10)
(2.1.2) u satisfies the Gårding inequality (a weaker coercivity): for any v(t, E V, define
B(v, w; u) == (0.v, w) vVVv, Vw) + ((u w),
and for v in neighborhood such that Assumption (2.1.)) holds, there exists a > 0 such that v,
d u 0| z |10 a||V
(E.11)
dt
The dual norm of the residual then reliable to serve as an error measure the following sense: denote Tm == (fm (m+1)
||ue uN 11 (Tm:H)) + +uu au--wwllerre uNlli I [1(u uN)((mm))) +C cff (R(ux)(tt-)Il) dt.
(E.12)
The constant C = C(v).
Proof. The proof of part II of Theorem 3.1 combines the meta-framework of Verfürth (2003) for time-dependent problems and Verfürth (1994, Section 8) for the stationary NSE. Note that thanks for the construction of divergence-fiee uvy the technicality of applying the Ladyzhenskaya Babušška- Budyzhenskaya-BBauukka Brezzi inf-sup condition avoided. In the meantime, no interpolations are needed to convert the functional norm for the negative Sobolev spaces an L representation.
To prove the upper bound (reliability) (3.5), we simply choose a time step t tm. and let (u (t,),,,) on (4m,;/m+1) (m+1). Note that the discretized approximation' 's evaluation which may
32

---

Published as a conference paper at ICLR 2025
not be the temporal grids im fm+1 naturally defined using the basis in (3.1). Using the error representation in (E.9), we have the (I) and (II) terms (R(uM (R(2)); ; u uNx are
(1) + (II) (Or(u ux): u ux) + (w V -uN), v(u - uNr)
d
(E.13)
Iu u- + v ux)
2 dt
For the (III) and (IV) terms we have
by Lomma
=((+)
(III + (IV) = (uu Vu u uN), u uN) + (((u - uN) V uN a uN).
(E.14)
Note by Lemma E.2, (+) is also:
(+) (1(u -uy) V) u,u u uN):
(E.
since 222,2, O for = u - uN. Consequently, by vector calculus identity, we have
J (= V)u zdx / Vu: (2)0z)
n
Thus, we reach the following error equation:
d d
(R(ux), a uN) Iu uN Vw Vu ((u ux) uN)) dæ.
(E.16)
Now, by Assumption (2.1.2). Poincaré inequality from Lemma E.3,
1 d
d
uu uv/l + vlu ((muvl))11
|u uNl1 B)] (-luuuvv) ux)
2 dl
dl
||u uvl1? + +all - ( -Blla url
dt
B(u
uw-u-uxie
By (E.13), (E.14), and (E.15),
B(u - uN, a uN;u) = (I) + (I) + (II) + (IV) = (R(ux), u-1) unn
< [R(uv))l (i(uw)ll--//a--11 uN/w1/2
s ž R(uN) 11 j/ue uN
As result, we have Integrating from fm to any t E (tms tm+1] yields
llu uvl1 + - uv)i [ di
(E.17)
[" ([R(uN)(4)))3311)) di.
(l, ((lluuuu((mm + (2) (t.. S
Since é (tm: 1m+1] 1] is arbitrary, we have proved the reliability (E.12).
Remark E.7 (Necessity of the "sufficient close" and the regularity conditions). First, the main hurdle prove the reliability that the convection term (*) is not positive definite, one would encounter some difficult in deriving the upper bound as yet moving (*) ) to the right-hand side in (E.16) does not yield a meaningful estimate,
of the error (1) (11) (R(ux), uxy) (+) & (R(ux), uxxx)
We note that Assumption (2.1.2) is equivalent putting threshold ||u(t, 2,,,, which is
convection-doiintted commonly assumed the analysis of numerical methods for convection dominated problems. stronger alternative assumption than Assumption (2.1.2) would be imposing a stronger constraint
33

---

Published as a conference paper at ICLR 2025
in Assumption (2.1.)) Using the "sufficient close" assumption (E.10) the nonlinear term, we have
[()1 S 11 uxx) V) uyullv
- - S (aluuuuxxx
Now if := - O, one would reach a similar estimate as (E.17) by replacing v with v. If no extra assumption y is imposed.
Another way is to consider a Sobolev embedding directly after applying the Hölder inequality, for (E.15) we have
((+) S Vuluu- uvl 2 Vull u uN S PllVullllV(e uux)=));
where p is the constant in the following Sobolev embedding
Iu uNvll S (u (xx)llll))l
Now if = Vul 2 O, one would reach similar estimate as (E.17) by replacing v with i.
E.2 PROOF THEOREM 3.2
Theorem 2.2 (Equivalence of functional norm and negative Sobolev norm), Consider the Gelfand triple V E H @ V", then
Ilfll lfl-1 for H/RR.
(E.18)
Proof. This proof leverages the spectral basis for H without the extra technicality of Schwartz space involved one ought to assume that f E V'. Recall (3.1), and define an infinite version
Sx span { = k 2x
(E.19)
In what follows, shall omit the extra technicality that one has to work on partial sums first, and then considers the convergence in the corresponding norms. We simply take for granted that the differentiation and integration/sum can be interchanged, and directly identify f H/R with its Fourier series in S
f(x) = Σ
kE
Note since f lies in the quotient space, f(0, 0) = 0. Then, the duality pairing can be identified as an L2- product, well the test function 0 H.
(f.u)viv (f.u)
Σ KE2RZA
fk
5(m)e mim
mE2n22
Σ
((K))(()
ME2*73 (0)
1/2 (1=2((()))))))
1/2
Σ
k
(k)))
k
KE2m ((0)
< 1/1-1/hh-
As result, by the definition in (3.7) and the estimate above
(f.v)
If
sup DEV/R.lelV=1
Kf.
sup VEV/R
< 1f
|p|1
34

---

Published as a conference paper at ICLR 2025
To show the other direction, let v (-A) f H' (T²) by the well well-posedness of -A f for E (T²). We have
(-(-)) (-Av, (fi(-A)-1))
Σ 22
f(b)e
E
mE2m2²
Σ e the
Σ m
emm
ke2n22
mE2n22
k ( 12 111
kE2m22 \(0)
Note that by the construction of v above, 10/1 = 11 thus,
(f
(fi(-A)-1)) 1ff
sup .EV |0/1
|0/1
which proves the theorem.
E.3 CONVERGENCE RESULT OF FINE-TUNING
In the following theorem, Theorem E.9, we show that sufficient condition for the optimizer to converge to "reach" a neighborhood of the true solution, thus corroborating the necessity of training. Note in both Theorems 3.1 and E.9, the eπrоr term of uN" to acknowledge that the "initial value" for the predicted trajectory, which is the last snapshot in the
(m) u(tm: -) is present. We have
input trajectory, may have errors. This is the major motivation that we opt to use u which is the input trajectory's 's last snapshot in the skip-connection in Q
(m)+) uN
=u²² + @o(oS(tt)) for = l,... 1
(E.20) In some sense, ST-FNO learns the derivative 0,u's arbitrary-Iength integral. one wants to modify Algorithm 1 in view of Theorem E.9 such that the error control is guaranteed, the following algorithm can be used but loses the "parallel-in-time" nature. We also note that thanks to the spectral convolution in Q being affine linear, showing Theorem E.9 is quite straightforward. as one has to establish the connection between fine-tuning and seeking a nonlinear Galerkin projection in Fourier space (3.1) under the functional norm. Let 0" = arg ming (u)) -ll-1. then Qo- (uN) arg minlres [|uu((+1,)) (u(Im,)) :) + w)|l-
Lemma E.8 (Local strict convexity for the fine-tuning loss). Define ||. ll to be a (dual) graph norm on - (T:: (T())) n L2 (T;; H
( (T²)) where Ta = jt - d, + 8]
||011,,8 = {fr |10(t, (1((,,) + fr 11 Vo) (t. I1 dt
1/2
For u E 1 (Ts; (T2 n (Tx; HT2))) the weak solution to (2.2) on To that is sufficiently smooth, there exists ò, e E Rt, such that on
B(u;e)= = {u L (Ts; L L' (T H (T²)) I|u 01|,,6 S se}.
the functional
J(u(t.))) 1 [R(()) where R(v) = 0,,0 (v-V)u vAv
is strictly convex.
Proof. First, by Theorem 3.2,
1
J(v) = 2
35

---

Published as a conference paper at ICLR 2025
Algorithm 2 An error-guarantee fine-tuning strategy.
Input: ST-FNO Go,e == Qo Ge: time stepping scheme (-); optimizer D(O, Vo(-)); training
dataset: solution trajectories at ,f.....,te" (e) as input and at [I.+1...., :...:...,,,,/ 4n; as output.
1: Train the ST-FNO model until the energy signature matches the inverse cascade.
2: Freeze ④ in Ge of ST-FNO Go e, keep Qe trainable.
3: Cast all nn. Module involved and tensors to torch. float64 and torch. complex122 hereafter.
Input: Evaluation dataset: solution trajectories at [f1,...,44 [g] as input, output time step 74.
4: for m = C...
f4 1 do
5: Extract the latent fields wx (m++) output of Ge at 1m+1 and hold them fixed.
(m+1)) (0) == uN + (0x (m((+))).
6: By construction of ST-FNO: such that uN
7:
March one step with (At)² using By: (0) = (B, - (N)(
8:
<-0 0
9:
while 7/m (O)( D,,x Tol do
10:
Apply the optimizer update parameters in Q: - D(O, Vo(mm)) +1.
Forward pass only through Q to update uN t ux + Qa(uN (m++))
11:
12:
if > Itermax then break
(m++1) uN N
13:
(m) h t+n, Output: A sequence of velocity profiles at corresponding time steps (uN mm=tt1
Then,
J(u + rE) - J(0)
DJ(u;E) == lim T-0
J(U rE)
7
dr
= (DR(v)) (-A)-(R(o))ve))
where
DR(0)E = 0,5 + 1 ( V) (v. V) vÂE the Fréchet derivative DR(v): V V" by Lemma E.6. The Hessian is then
Hess [(;; &, () (DD((v), (- -D(()) (s (- Now, since R(u)= = on Ts. we have Hess (w; E, (DR(U)&, (- DD(()EE)
we assume that DR(u){ e V has its functional norm bounded above and below in B(u u: e). one has for any E
| $ | S Hess J(u; E)
Simply choosing small enough such that v is sufficiently close to in graph norm associated with the PDE to make the coercivity above still true for Hess J(v: E. {) yields the desired local convexity.
Theorem E.9 (Guaranteed convergence of the fine-tuning). In addition to the same assumptions with Theorem 3.1, suppose Lemma E.8 holds for E (0, 1), and given uy be embedded in forad a 0 < é' S c. Denote uN the evaluation in Line 9 of Algorithm 2, and j the iteration Blu; for of optimizer in Line 10, then the fine-tuning using the new loss function (3.9) produces a sequence fuN C B (u;ee) Furthermore, suppose that the optimizer in Line 10 of Algorithm 2 has learning rate converging to O. Then, then the fine-tuning using the new loss function (3.9) produces
(k)
sequence of evaluations converging the best possible approximation (m=+1) S of u(tm-1:-) NN
[m)) starting from u. the sense that for m .,.... E + n4
-2 [1((tm+1) )222 + c2n, (At)7-1 (E.21) N00000 u((mm1 -bl < | u(lm, + cin.n
(m+1)
(m) S l S. First, by Proof. For simplicity, we denote un := u(t),)) .) E V. and uv u x
line 10 fine-tuning algorithm, if one solves the optimization in the functional norm exactly, we have
8" arg min öt (R(Buux(() | where ux(() = uv Qo(ww)
(E.22)
36

---

Published as a conference paper at ICLR 2025
where R(-) is the residual functional computed using the time derivative term from extra-finee fine-step solver's result
R(v) DDv (v
+/u+
Unlike representing the derivative using output from the neural operator, one of the keys of our algorithm is that the error for 8,u D can be explicitly estimated using the framework to develop estimates for truncation emror in traditional time marching schemes for NSE. Due to the choice of the time step, this truncation error will be of higher order. Note for any uv. (o(ww) E S, as a result, solving (E.22) is equivalent to solve the following: denote Ug == (uN gu) for du E S
min (6t) 1/2 (us(us) + (us V)us vAus flls. AuES
Replacing D.(.) by (.(.) we have truncation error term, whose error of order (At)7-1 due to taking the time derivative:
[Dr(us) = (us V)us - vλus - fllss
S [10eus (us V)us vAus fllss + (1D((us) a asullss.
We focus estimating the first term above,
min (8t) 1/2 ||0,us + (us V)us vAus ôuES
flls
(E.23)
By the fact that || Ils inherit the scaling law and the triangle inequality from || ll, it is convex as well in this neighborhood. As result, any gradient- based optimizer with a linearly converging step
(m)) size shall converge to the minimum, achieved at N. with linear convergence rate. Moreover,
(m) (m²)
uu (dus) is the (nonlinear) Galerkin projection of u(lm; E V.
N
37


Reviews:
[{"review_id": "TB8GNqVHx5", "reviewer": "Reviewer_UD1B", "summary": "This paper introduces a new learning framework to improve operator-type neural networks for solving spatiotemporal PDEs. The proposed spatiotemporal adaptation generalizes FNO to enable temporal super-resolution. By refining traditional end-to-end training with insights from numerical PDE theory, the framework trains an FNO briefly and then fine-tunes a new spatiotemporal spectral convolution layer without frequency truncation, using a novel negative Sobolev norm loss function.", "questions": "1. The paper notes a spatial grid size of 256×256 and a prediction time interval from \\( t = 4.5 \\) to \\( t = 5.5 \\), which is quite short for meaningful evaluation. This limited time horizon may not fully showcase the model’s performance, and longer predictions should be considered for comparison. Could it surpass top-performing models such as PDERefiner, TSM, and LI in terms of prediction accuracy?\n\n2. While this method targets turbulent modeling of Navier-Stokes equations, its applicability to other nonlinear or multiphysics PDEs is not discussed. Could Spectral-Refiner be adapted to handle other PDEs, such as heat conduction or elasticity? What modifications or enhancements would be required?\n\n3. The paper uses the H⁻¹ Sobolev norm as a tool for error estimation, yet traditional numerical analysis often questions the suitability of non-local norms. Without access to exact solutions, is this norm sufficiently stable? How well does this error estimation perform with initial conditions that lack global smoothness?\n\n4. This method is suggested to be more efficient and stable than physics-informed neural operators , but the comparative experiments are limited. Could more comprehensive experimental data, including comparisons with popular neural operators like DeepONet, be provided?\n\n5. The paper notes the removal of frequency truncation during fine-tuning, which may enhance the model’s ability to capture low-frequency information. However, does this increase the risk of amplifying high-frequency noise? If so, what measures could control this issue without compromising accuracy?\n\n6. The method involves high-order Fourier transforms and complex error estimations, potentially adding significant computational cost. In practical applications, especially those requiring real-time simulation, can this method meet timing constraints? How does its inference time compare with other methods?\n\n\n7. How different are the training and testing initial conditions?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.The ST-FNO extends FNO capabilities to handle arbitrary temporal and spatial resolutions, improving flexibility and applicability for complex PDEs like NS equations.\n\n2.By designing the ST-FNO as a zero-shot model with arbitrary-length temporal inference, the Spectral-Refiner can adapt flexibly to varying time horizons, making it well-suited for large-scale or long-term simulations without significant retraining.", "weaknesses": "1. Although this paper claims that Spectral-Refiner outperforms traditional numerical methods in accuracy and computational efficiency, it does not specify how the method performs under different resolutions, time steps, and initial conditions.\n\n2. Can the model generalize to varying conditions, such as changes in Reynolds numbers or external forces? If so, please provide correlation curves and energy spectrum curves under these conditions.\n\n3. Model performance should not be evaluated solely based on the final frame; an error propagation curve and error distribution plot are needed to illustrate the model’s overall performance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "cQm5ldUifE", "reviewer": "Reviewer_9ebp", "summary": "This paper introduces the spatiotemporal adaptation technique for all FNO variants (ST-FNO), enabling them to learn mappings between Bochner spaces. The authors propose a novel strategy for training and evaluating ST-FNOs, which surpasses existing methods in both speed and accuracy. Numerical experiments on Navier-Stokes Equations (NSE) benchmarks demonstrate substantial improvements in computational efficiency and accuracy.", "questions": "Can ST-FNO be applied to other equations and still perform effectively?\n\nIn line 65, the phrase “Similar to the traditional time marching solvers” is confusing. Did the authors intend to mean \"Unlike traditional time marching solvers\"?\n\nIn lines 217-218, what is the definition of $q$?\n\nIn lines 1678-1679, there is a '??'. Please correct this.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The concepts of the spatiotemporal adaptation technique and the hybrid operator learning paradigm appear novel.\n- The method is solidly supported by theoretical evidence.\n- Overall, the paper is well-written and clear.", "weaknesses": "- The authors assert that their methods \"yield accuracy comparable to traditional numerical methods, yet with computational resources akin to evaluating NOs.\" This claim would indeed be exciting if validated. However, the paper's experiments appear insufficient in the following aspects:\n\n  - For both FNO and ST-FNO, while they may perform well on the training data distribution, they likely fail when tested with out-of-distribution data. In extreme out-of-distribution scenarios, can ST-FNO still maintain accuracy comparable to traditional numerical methods within a limited fine-tuning timeframe?\n\n  - Table 4 compares the FLOPs of different methods, but the running times for each method are absent. I am concerned about whether the sum of inference and fine-tuning times for ST-FNO is comparable to the inference time for FNO alone.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xCxdYnDuYG", "reviewer": "Reviewer_5tCD", "summary": "This paper presents a novel learning framework that enhances the capabilities of operator-type neural networks, specifically focusing on Fourier Neural Operators (FNOs) for solving spatiotemporal Partial Differential Equations (PDEs). Recognizing the high training costs and variable accuracy of existing models, the authors introduce a spatiotemporal adaptation that allows FNOs to learn mappings between Bochner spaces, enabling arbitrary-length temporal super-resolution. The framework integrates insights from traditional numerical PDE techniques to refine the conventional end-to-end training process. For turbulent flow modeling with the Navier-Stokes Equations (NSE), the approach involves initial training of the FNO for a limited number of epochs, followed by fine-tuning a new spatiotemporal spectral convolution layer without frequency truncation. A unique fine-tuning loss function using a negative Sobolev norm, defined through a functional a posteriori error estimator based on the Parseval identity, is introduced. This loss function simplifies the optimization process, as it is convex, contrasting with the nonconvex challenges typical of end-to-end training. The proposed method significantly improves both computational efficiency and accuracy in numerical experiments on standard NSE benchmarks, outperforming traditional numerical PDE solvers and end-to-end evaluations.", "questions": "Could you elaborate on the motivation for learning maps between Bochner spaces? How does this choice influence the model’s performance and its ability to handle temporal super-resolution?\n\nGiven the model’s complexity, how interpretable are the results?\n\nHow sensitive is the model to hyperparameter choices, especially those related to the spatiotemporal spectral convolution layer? Are there guidelines for selecting optimal hyperparameters for new applications?", "limitations": "", "rating": 5, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper introduces a novel learning framework that effectively combines traditional numerical PDE techniques with operator-type neural networks, enhancing the capabilities of Fourier Neural Operators (FNOs) in solving spatiotemporal PDEs. The ability to generalize FNO variants to learn maps between Bochner spaces and perform arbitrary-length temporal super-resolution is a significant advancement for dynamic environments. By proposing a training strategy that involves limited initial epochs followed by fine-tuning, the authors address the challenges of high training costs while maintaining performance. This approach can make training more accessible and efficient.\nThe introduction of a fine-tuning loss function that is convex simplifies the optimization process, reducing the complications often associated with nonconvex problems in neural network training. The framework demonstrates significant improvements in computational efficiency and accuracy through numerical experiments on commonly used Navier-Stokes Equation benchmarks, providing strong empirical support for the proposed methods.", "weaknesses": "As with many neural network-based approaches, the interpretability of the results may be challenging.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nWGG5h84Xq", "reviewer": "Reviewer_rmBG", "summary": "This paper according to me addresses limitations in current neural operator approaches for modeling turbulent flows governed by the Navier-Stokes equations. The key contributions are:\n1. A new spatiotemporal adaptation (ST-FNO) of Fourier Neural Operators to enable learning maps between Bochner spaces, allowing arbitrary-length temporal predictions.\n2. A novel training-fine-tuning paradigm that combines limited epochs of end-to-end training with targeted fine-tuning of a spectral convolution layer.\n3. A new loss function for fine-tuning based on a functional-type a posteriori error estimator using a negative Sobolev norm, which is reliably evaluated through Parseval's identity.\n4. Empirical demonstration of significant improvements in both computational efficiency and accuracy compared to existing methods.\nOverall I like the paper! It's important for a lot of problems that require FNO to predict multiple temporal steps, especially in high Reynolds number ranges for NS.", "questions": "Just some questions:\n1. Computational efficiency trade-offs: Can you provide a more detailed comparison of the computational costs of your method (including the fine-tuning phase) versus traditional CFD methods and pure neural network approaches? Specifically, how does the accuracy-to-compute-time ratio compare across these methods for different problem scales?\n2. Error accumulation in multi-step predictions: For multi-step predictions, how does the error accumulate over time compared to traditional numerical methods? Is there a point at which the accuracy degrades significantly, and if so, how might this be mitigated?\n3. Spectral layer fine-tuning: You mention fine-tuning only the last spectral convolution layer. Have you experimented with fine-tuning other layers or using a different architecture for this layer? How sensitive is the method to the design of this final layer?\n4. Arbitrary-length temporal predictions: Can you elaborate on how your ST-FNO handles arbitrary-length temporal predictions? Specifically, how does the performance change as the prediction length increases beyond what was used in training? Are there any limitations on how far into the future the model can reliably predict?\n5. Handling of boundary conditions: How does your method handle different types of boundary conditions, especially for the vorticity-streamfunction formulation? Can it adapt to changing boundary conditions without retraining?\n6. Temporal super-resolution: Can your method perform temporal super-resolution, i.e., predict at a higher temporal frequency than the training data? If so, how accurate are these interpolated predictions?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. Theoretical foundation: I like that the paper provides a rigorous mathematical analysis of discretization mismatch errors and proves the reliability of the proposed error estimator.\n2. Novel architecture: The ST-FNO adaptation enables flexible spatiotemporal predictions, addressing a key limitation of existing neural operators. This is extremely useful in various timestepping problems.\n3. Efficient training paradigm: The proposed approach of limited training followed by targeted fine-tuning offers a computationally efficient alternative to standard end-to-end training.\n4. Strong empirical results: The method demonstrates significant improvements in accuracy and efficiency on challenging turbulent flow benchmarks.\n5. Open-source implementation: The authors provide code to reproduce their results, enhancing reproducibility and potential impact.", "weaknesses": "1. Limited scope: The method is primarily designed for rectangular domains and uniform grids, which may limit applicability to more complex geometries.\n2. Dependency on numerical solvers: The fine-tuning process relies on traditional numerical solvers for computing extra field variables, which may introduce computational overhead.\n3. Assumption sensitivity: The theoretical guarantees rely on certain assumptions about the closeness of solutions, which may not always hold in practice, I'm not sure how much it does in practice.\n4. Long-term stability analysis: While short-term performance is demonstrated, a more extensive investigation of long-term stability for time-dependent problems would be valuable. This could include chaos indicators or Lyapunov exponent analysis for longer time horizons, particularly for highly turbulent regimes.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)