Spectral-Refiner: Accurate Fine-Tuning of Spatiotemporal Fourier Neural Operator for Turbulent Flows

## Abstract

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 Spatiotemporal Fourier Neural Operator (SFNO) that learns maps between Bochner spaces, and a new learning framework to address these issues. This new paradigm leverages wisdom from traditional numerical PDE theory and techniques to refine the pipeline of commonly adopted end-to-end neural operator training and evaluations.
Specifically, in the learning problems for the turbulent flow modeling by the Navier-Stokes Equations (NSE), the proposed architecture initiates the training with a few epochs for SFNO, concluding with the freezing of most model parameters.
Then, the last linear spectral convolution layer is fine-tuned without the frequency truncation. The optimization uses a negative Sobolev norm for the first time as the loss in operator learning, defined through a reliable functional-type a posteriori error estimator whose evaluation is almost exact thanks to the Parseval identity.
This design allows the neural operators to effectively tackle low-frequency errors while the relief of the de-aliasing filter addresses high-frequency errors. Numerical experiments on commonly used benchmarks for the 2D NSE demonstrate significant improvements in both computational efficiency and accuracy, compared to end-to-end evaluation and traditional numerical PDE solvers.

## 1 Introduction

In recent years, there has been a rapid increase in the number of remarkable results obtained with Deep Learning (DL).
Notably, DL has proven particularly effective in addressing problems formulated by Partial Differential Equations (PDEs).
In this paper, we study the problem of learning operators between infinite-dimensional function spaces [, , , , ]. We propose to synergize operator learning with classical numerical PDE methods, complementing one’s drawbacks and limitations with the other’s strengths.
The problem in consideration is using Navier-Stokes Equations (NSE) to model turbulent flows. For fluid, a relatively large Reynolds number ($\mathcal{O}(10^{3})$ to $\mathcal{O}(10^{4})$) drives the flow to become turbulent, mainly attributed to the nonlinear convection. In a computational framework, the difficulties in simulating NSE roots from its “stiffness” and the spatiotemporal nature of being transient.
In the operator learning problem associated, the mapping between the spaces containing trajectories of solutions is approximated through operator-valued neural networks (NNs), or neural operators (NOs), 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.

Figure: Figure 1: Schematic difference between traditional time-marching solvers, “rollout” neural operators, and SFNO. Connection widths approximately represent computational cost (not drawn to scale). 4th-order Runge-Kutta (RK4): small time steps bounded by the CFL condition, de-aliasing filter needed; autoregressive NOs: using previous evaluation as input, large time steps, no stability guarantees; Spatiotemporal FNO with hybrid fine-tuning: large time steps, arbitrary-lengthed ground truth as input to give evaluations, fine-tuning stability inheriting from traditional solvers, no de-aliasing filter preserves all higher modes, parallel-in-time optimization. See Appendix A for notations, and Section 1.1/Appendix B for detailed discussions.
Refer to caption: x1.png

The prevalent functional analytic framework for studying the convergence and stability of these solution trajectories is the theory of Bochner spaces (e.g., Aubin-Lions lemma [, Chapter 3], [, Chapter 2], [, Chapter 7]).
This theoretical insight motivates this paper’s new design of neural operators. This new neural architecture that we design in this paper is based on the state-of-the-art FNO3d, introduced initially in []. While FNO3d learns maps between a fixed number of spatial “snapshots” (product spaces), the new Spatiotemporal Fourier Neural Operator (SFNO) proposes to learn a map between Bochner spaces $L^{2}(\mathcal{T}_{0};\mathcal{V})$ to $L^{2}(\mathcal{T};\mathcal{V})$ on non-overlapping time intervals $\mathcal{T}_{0}$ and $\mathcal{T}$ directly.
The learned operator provides increased flexibility for the temporal variable, allowing its dimension to vary, in addition to an arbitrary spatial discretization size. Here $\mathcal{V}$ denotes the spatial Hilbert space in which snapshots of the solution at a specific time reside, and for a detailed notation list, we refer the readers to Appendix A.

#### Main contributions.

The main contributions of this work are summarized as follows.

- •
Spatiotempoeral Neural Operator. We design the first neural operator SFNO that can learn maps between Bochner spaces.
- •
New hybrid paradigm. We propose to train and evaluate SFNO using a new strategy in spatiotemporal predictions, which has significantly improved over the existing methods in speed and accuracy. Unlike the arduous training of running the hundreds of epochs, we train SFNO only a few epochs (e.g., $10$), and then the *last* layer of SFNO is fine-tuned with the help of traditional solvers to obtain highly accurate approximations for extra field variables to conform with the physics.
- •
Functional a posteriori error estimation. Leveraging the spectral neural architectures in SFNO, Parseval identity, and Gelfand’s triple, we train and fine-tune the models using a new loss through the *a posteriori* error estimation in the negative Sobolev norm (functional norm). This procedure does not require any ground truth data (e.g., numerical solution of the underlying PDEs generated by traditional numerical solvers), nor the true solution(s) to the PDE. This new loss is proven to be reliable in theory, in the meantime much more efficient in the ablation experiments.
- •
Experimental results. We created a native PyTorch port of Google’s Computational Fluid Dynamics written in JAX [, ] with enhanced functionality for tensor operations such as the facilitation of fine-tuning for the latent fields, publicly available at [https://github.com/scaomath/torch-cfd](https://github.com/scaomath/torch-cfd), with scripts to replicate the experiments. The data are available at
[https://huggingface.co/datasets/scaomath/navier-stokes-dataset](https://huggingface.co/datasets/scaomath/navier-stokes-dataset)
.

### 1.1 Related Work and Motivations

Here we summarize the studies that motivate us to propose this new approach. For a more detailed review of operator learning and its analysis and applications, and further motivation discussion with a much higher degree of mathematical rigor, we refer the readers to Appendix B.

#### Compromises and drawbacks of traditional numerical methods.

The NO approach has the potential to bypass various difficulties and compromises of traditional schemes of numerical PDEs. In traditional numerical methods that use explicit time stepping, or there exists a certain portion of the forcing terms computed by explicit schemes, such as Adams-Bashforth or Runge-Kutta (RK) families (e.g., see [, Appendix D], []), the “stiffness” refers to the necessity for extremely small time step sizes ($\Delta t\sim\mathcal{O}(h^{\alpha})$, $\alpha\geq 1$, see e.g., [, Chapter 4]). This requirement, usually referred to as the Courant-Friedrichs-Lewy (CFL) condition (e.g., see [, , ]), poses a “stability” constraint on the step size to be much smaller than the requirement for “consistency”. Here, this temporal consistency usually refers to a first-order optimal local truncation error, e.g., how the original Butcher tableau is derived for RK [].
The CFL puts a threshold on how fast the local errors in different frequencies can propagate and prevent the accumulation “pollutes” the approximation globally to ensure stability. This means that, for finer mesh (better spatial consistency), the time steps have to be much smaller to prevent errors from traveling to neighboring nodes and elements.

#### De-alias filter sacrifices accuracy to maintain stability.

One famous example of the consistency-stability trade-off is the $3/2$-rule (also known as $2/3$-dealiasing dealiasing filter) for the nonlinear convective term [, , , ] for pseudo-spectral methods [, ]. The highest $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 to be stable, keeping the balance between stability and accuracy. These constraints apply to traditional numerical methods because any solver has to march a consecutive multitude of time steps, which makes the error propagation operator’s norm a product of many.
In this study, we explore a “hybrid” approach that combines the strengths of NOs and traditional numerical solvers.
In our approach, there is no time marching consecutively for a multitude of time steps, which renders the method free of the stability constraints such as the CFL condition and the $3/2$-rule.

#### Spatiotemporal operator learning.

Among the end-to-end spatiotemporal operator learners, the common approach is an autoregressive “roll-out”, similar to the traditional time marching solvers, yet with a much bigger time step size not restricted by the CFL condition. For example, in the original FNO2d [] and [, , , ], several snapshots of solutions $\{u^{(i)}\}_{i=1}^{m}$ are concatenated as the input to the NO, which outputs $u^{(m+1)}$ at the subsequent time step. This procedure can be repeated recurrently until reaching the final time step $t_{M}$. However, this type of approach suffers from error accumulation experimentally (e.g., see [, Figure 9]).
To our best knowledge, NO-only approaches have no theoretical stability estimate, such as the error propagation operator is a contraction. Recently, a remarkable advancement called PDE-refiner [] learns an error correction NN under the Denoising Diffusion framework [] to get a stable long roll-out and to achieve for the first time $\mathcal{O}(10^{-8})$ relative difference with the ground truth after a single time marching step. For the difference between the “roll-out” approach versus direct spatiotemporal learning between Bochner spaces, see Figure 1.

#### How can NO evaluations get closer to the true solutions of PDEs?

In all the models above and their learning frameworks, the optimization is to minimize the difference between the outputs from the NO, namely $u_{\mathcal{N}}$, and the ground truth $u_{\operatorname{\mathcal{S}}}$, generated by a traditional numerical PDE solver.
The difference with the true solution $u$ under a certain norm is not directly optimized, since for real-life applications the analytical solution $u$ is unknown. For difficult PDEs such as the NSE, in a linear time stepping scheme, the ground truth $u_{\operatorname{\mathcal{S}}}$ may on its own only have 3-digit accuracy in Bochner norm. Note that this may already be of optimal order $\mathcal{O}(\Delta t)$ by convergence estimates []; Moreover, in the end-to-end operator learning benchmarks of NSE, even the state-of-the-art NOs still fall short in achieving high-accuracy solutions. To our best knowledge, no NO-only-based operator learning approach can break the barrier of a 2-digit accuracy in terms of the relative difference with $u_{\operatorname{\mathcal{S}}}$ under the Bochner norm. Nevertheless, we have the following enlightening Pythagorean-type identity: for the true solution $u$ satisfying $A(u)=f$

$$ $\|u-u_{\operatorname{\mathcal{S}}}\|_{a}^{2}+\|u_{\operatorname{\mathcal{S}}}- u_{\mathcal{N}}\|_{a}^{2}=\|u-u_{\mathcal{N}}\|_{a}^{2}\quad\text{ since }(u-u _{\operatorname{\mathcal{S}}})\perp_{a}(u_{\operatorname{\mathcal{S}}}-u_{ \mathcal{N}})\in\mathcal{V}.$ (1.1) $$

For the sake of presentation, we make some handy assumptions to get this identity: (1) $A:\mathcal{V}\to\mathcal{V}^{\prime}$ denotes a linear operator on the Gelfand triple;
(2) $u_{\operatorname{\mathcal{S}}}$ is obtained through Galerkin methods [], which is a projection onto a finite-dimensional approximation subspace of $\operatorname{\mathcal{S}}\subset\mathcal{V}$; and (3) an evaluation of $u_{\mathcal{N}}$ can be continuously embedded into $\operatorname{\mathcal{S}}$. Through the Riesz representation theorem, $u_{\operatorname{\mathcal{S}}}$ is solved through

$$ $(u-u_{\operatorname{\mathcal{S}}},v)_{a}:=\langle A(u)-A(u_{\operatorname{ \mathcal{S}}}),v\rangle=0,\quad\text{ for any }v\in\operatorname{\mathcal{S}}.$ (1.2) $$

This bilinear form $(\cdot,\cdot)_{a}$ induces a (semi)norm $\|\cdot\|_{a}$ inheriting the topology from $\mathcal{V}$.
Given this orthogonality, minimizing the difference between $u_{\mathcal{N}}$ and $u_{\operatorname{\mathcal{S}}}$ becomes fruitless if $\|u-u_{\operatorname{\mathcal{S}}}\|_{a}$ is relatively big in the first place, and unnecessary computational resources may have been spent to get closer to $u_{\operatorname{\mathcal{S}}}$. Rather, (1.1) indicates that it is more efficient if one can design a method to reduce the error of $\|u-u_{\mathcal{N}}\|_{a}$ directly, while circumventing the fact that $u$ is not accessible.

#### Functional-type a posteriori error estimation.

The *a posteriori* error estimation, computing the error without knowing the true PDE solutions, has been widely studied for Galerkin-type methods, such as finite element [, , ], as well as for parameterized PDEs [, ]. Among all types of *a posteriori* error estimation techniques, functional-type *a posteriori* estimator [] views the error as a functional on the test Sobolev spaces and evaluates accurate representations with the help of a dual variable [].
In our study, for a mesh-invariant method such as neural operators, we propose to use a negative Sobolev norm in the frequency domain as a functional-type *a posteriori* error estimation as loss. It leads to “refining” the learnable set of the global spectral basis, thus not be bounded by many constraints needed for the local adaptive mesh refinement [] to inaccurate localized representations for the $H^{-1}$ error functional. The negative Sobolev norm used in our new loss is handily defined through the Gelfand triple in the frequency domain, no extra flux reconstruction is needed either.

## 2 Spatiotemporal Operator Learning for Navier-Stokes Equations

Both drawbacks and advantages of traditional numerical methods and NO-based methods in Section 1.1 play vital roles in shaping our study.
We first briefly discuss the spatiotemporal operator learning problem on Bochner spaces associated with NSE. Then, we detail how we modify the Fourier Neural Operator 3D (FNO3d) architecture to obtain the first operator learner between Bochner spaces.
Next, we shall rethink the efficacy of traditional end-to-end training strategies for NOs. In the last part of this section, a new training-and-fine-tuning paradigm is proposed to yield accuracy on par with traditional numerical methods, yet with computational resources comparable to the evaluation of NOs. This new approach does not need hundreds of thousands of marching steps like traditional numerical methods. We shall see that the new method leverages the efficiency in NO evaluations. At the same time, it also 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 numerical methods must make to ensure stability and convergence.

### 2.1 Operator learning problem for NSE

For 2D NSE, the velocity field $\bm{u}(t,\bm{x}):[0,T]\times\Omega\to\mathbb{R}^{2}$ is seen as an element in the Bochner space $L^{p}([0,T],\mathcal{V})$ where $\mathcal{V}$ is a spatial Sobolev space in which each snapshot at $t$ of the solution $\bm{u}(t,\cdot)$ resides.
The operator learning task for NSE using FNO3d [] is to learn a map $G^{\mu}$ between a fixed number of Cartesian products of spatial Sobolev spaces.
Specifically, $G^{\mu}:\mathcal{X}\rightarrow\mathcal{Y}$ maps elements in $\mathcal{X}=\mathop{\mathchoice{\textstyle\prod}{\prod}{\prod}{\prod}}
\nolimits_{i=1}^{\ell}\mathcal{V}$ to elements
in $\mathcal{Y}=\mathop{\mathchoice{\textstyle\prod}{\prod}{\prod}{\prod}}
\nolimits_{i=1}^{n_{t}}\mathcal{V}$, where $\ell\in\mathbb{Z}^{+}$ is *fixed* in this learning task, while $n_{t}$ can vary. As such, $\mathcal{X}$ and $\mathcal{Y}$ represent spaces of two non-overlapping temporal segments of solution snapshots. As in [], the NSE variant considered here is the vorticity-streamfunction (V-S) formulation (2.1) with periodic boundary condition (PBC). We also consider the velocity-pressure (V-P) formulation (2.2). In $\Omega\times(0,T]$, for vorticity $\omega:=\nabla\!\times\!\bm{u}$, and
streamfunction $\psi$, these two formulations read

$$ $\displaystyle(\text{Vorticity-Streamfunction})$ $\displaystyle\partial_{t}\omega+\operatorname{\mathbf{rot}}\psi\cdot\nabla{ \omega}-\nu\Delta{\omega}=\nabla\!\times\!\bm{f},\quad\omega+\Delta\psi=0.$ (2.1) $\displaystyle(\text{Velocity-Pressure})$ $\displaystyle\partial_{t}\bm{u}+(\bm{u}\cdot\nabla){\bm{u}}-\nu\Delta{\bm{u}}+ \nabla p=\bm{f},\quad\nabla\cdot\bm{u}=0.$ (2.2) $$

For all analyses in line with the Hilbertian framework, $\mathcal{V}=H^{1}(\mathbb{T}^{2})$ for vorticity and $\bm{H}^{1}(\mathbb{T}^{2})$ for velocity, where $\mathbb{T}^{2}$ denotes the unit torus, i.e., $\Omega=(0,1)^{2}$ with a component-wise PBC. For a fixed forcing function in $\mathcal{V}^{\prime}$, the initial condition is either $\omega(0,\bm{x})=\omega_{0}(\bm{x})$ or $\bm{u}(0,\cdot)=\operatorname{\mathbf{rot}}\bigl{(}(-\Delta)^{-1}\omega_{0}
\bigr{)}$. Here $\omega_{0}$ is drawn from a compactly supported probability measure $\mu$, in which the compactness corresponds to certain power/enstrophy spectrum decay law to produce isotropic turbulence [] (refer to Appendix C for details). Then, the map $\mathcal{G}:\mathop{\mathchoice{\textstyle\prod}{\prod}{\prod}{\prod}}
\nolimits_{i=1}^{\ell}\mathcal{V}\to\mathop{\mathchoice{\textstyle\prod}{\prod
}{\prod}{\prod}}\nolimits_{i=1}^{n_{t}}\mathcal{V}$ is as follows:

$$ $\mathcal{G}:\mathbf{a}:=\bigl{[}\omega(t_{1},\cdot),\dots,\omega(t_{{\ell}}, \cdot)\bigr{]}\mapsto\mathbf{u}:=\bigl{[}\omega(t_{\ell+1},\cdot),\dots,\omega (t_{\ell+n_{t}},\cdot)\bigr{]},$ (2.3) $$

where $t_{\ell+n_{t}}\leq T$, and the input and the output are snapshots obtained from a solution trajectory with the same initial condition $\omega_{0}$.
Based on these snapshots’ discretized approximations from data, FNO3d builds a parametrized map $\mathcal{G}_{\theta}:\mathcal{X}\to\mathcal{Y}$, with $\theta\in\mathbb{C}^{p}\times\mathbb{R}^{q}$, here $p$ is independent of the spatial discretization size, but $q$ does depend on the number of input snapshots $\ell$. The operator learning problem is then to seek $\mathcal{G}_{\theta^{\dagger}}$ through a standard end-to-end supervised learning pipeline to represent the distribution of the input-output pairs of $\mathcal{G}$ well enough. During testing (or inference), the FNO3d predicts the subsequent $n_{t}$ snapshots with the first $\ell$ snapshots as input, both of which may feature a much higher spatial resolution than the input-output pairs used in the training phase.

In this study, the learning aims to recast (2.3) to conform to the Hilbertian formulation of NSE:

$$ $\widetilde{\mathcal{G}}:L^{2}(t_{1},t_{\ell};\mathcal{V})\to L^{2}(t_{\ell+1}, t_{\ell+n_{t}};\mathcal{V})\;\text{ for }\;\mathcal{V}:=H^{1}(\mathbb{T}^{2}) \text{ or }\{\bm{v}\in\bm{H}^{1}(\mathbb{T}^{2}):\nabla\cdot\bm{v}=0\},$ (2.4) $$

since weak solutions to (2.1) and (2.2) at a given time interval $\mathcal{T}$ is in $L^{2}(\mathcal{T};\mathcal{V})$.
In what follows, we shall present how to design an operator learner to map arbitrary-sized discretization in $\mathbb{R}^{\ell\times n\times n\times d}$ to $\mathbb{R}^{n_{t}\times n\times n\times d}$ ($d=1$ in V-S; $d=2$ in V-P).
Henceforth, $\ell$, $n_{t}$, $n$ are all of variable sizes, and we omit the $d$ dimension if no ambiguity arises.

### 2.2 Spatiotempoeral Fourier Neural Operator

SFNO $\mathcal{G}_{\theta}$ features the same meta-architecture of FNO3d [] as follows.

$$ $\mathcal{G}_{\theta}:=\widetilde{Q}\circ S\circ\sigma_{L}\circ K_{L-1}\circ \cdots\circ S\circ\sigma_{1}\circ K_{0}\circ\widetilde{P}.$ (2.5) $$

In (2.5), $\widetilde{P}$ is a lifting operator, $\widetilde{Q}$ is a projection operator that does a pointwise reduction in the channel dimension,
$\sigma_{j}$ can be a pointwise universal approximator or simply chosen as a nonlinearity, and $S$ is an optional non-learnable Helmholtz layer (see Appendix D.2). $K_{j}:=K_{\phi_{j}}$ is the parametrized spectral convolution integral operator of Fredholm type.
During training of $\mathcal{G}_{\theta}$, the operator to be learned $\widetilde{\mathcal{G}}$ is restricted to finite-dimensional subspaces $\mathcal{X}\supset\mathcal{X}_{n}\to\mathcal{Y}_{n}\subset\mathcal{X}$, in the sense that the continuous spatial Sobolev space $\mathcal{V}$ in the product spaces is replaced by a finite-dimensional subspace $\mathcal{V}\supset\operatorname{\mathcal{S}}\simeq\mathbb{R}^{n\times n}$ with continuous embeddings $\{\mathbf{a}_{\operatorname{\mathcal{S}}}\in\mathbb{R}^{\ell\times n\times n}
\}\hookrightarrow\mathcal{X}$, and
$\{\mathbf{u}_{\operatorname{\mathcal{S}}}\in\mathbb{R}^{n_{t}\times n\times n}
\}\hookrightarrow\mathcal{Y}$. The positional encodings $\mathbf{p}_{\operatorname{\mathcal{S}}}:=(t_{i},\bm{x}_{j})_{i=1}^{\ell}\in
\mathbb{R}^{3\times\ell\times n\times n}$ represents the spatiotemporal grid, and is concatenated to $\mathbf{a}_{\operatorname{\mathcal{S}}}$ before feeding to $P$.

The spectral convolution integral operator $K(\cdot)$ in FNO3d features a slightly different architecture than most presentations in the literature that are meant for spatial tensor-to-tensor maps in NOs.
Slightly abusing of notation, denote the latent dimension (width) by $d_{v}$, let the vector-valued latent representation with $d_{v}$ channels be continuously embedded into $\mathop{\mathchoice{\textstyle\prod}{\prod}{\prod}{\prod}}\nolimits_{i=1}^{n_{
t}}\mathcal{V}$ by $\mathbf{v}_{\operatorname{\mathcal{S}}}\hookrightarrow\mathbf{v}$ in each channel, and
denote $\Lambda:=[t_{\ell+1},\dots,t_{\ell+n_{t}}]$, then $K(\cdot)$ is “semi-discrete” in a sense as follows

$$ $(K\mathbf{v})(t,\bm{x}):=W\mathbf{v}(s,\bm{y})+\sum\nolimits_{s\in\Lambda}\int _{\Omega}\kappa\bigl{(}(t,\bm{x}),(s,\bm{y})\bigr{)}\mathbf{v}(s,\bm{y}) \mathrm{d}m(\bm{y}),$ (2.6) $$

where $t$ takes values only on $\Lambda,$ and $\bm{x}\in\Omega$. Here $W\in\mathbb{R}^{(d_{v}+1)\times d_{v}}$ is a pointwise-applied affine linear operator, $\kappa\in C\bigl{(}(\Lambda\times\Omega)\times(\Lambda\times\Omega);\mathbb{R}
^{{d_{v}\times d_{v}}}\bigr{)}$, and $m$ denotes an approximation to the Lebesgue measure on $\Omega$.
Evolving $K$ into a spatiotemporal spectral convolution acting on Bochner space is straightforward. We change (2.6) slightly by adopting of an atom-like measure $\delta(\cdot)$ in the temporal dimension, which then generalizes to the variable-length temporal discretization of snapshots for any $(t,\bm{x})\in(a,b)\times\Omega$

$$ $(\widetilde{K}(\mathbf{v}))(t,\bm{x}):=W\mathbf{v}(s,\bm{y})+\int_{a}^{b}\int_ {\Omega}\kappa\bigl{(}(t,\bm{x}),(s,\bm{y})\bigr{)}\mathbf{v}(s,\bm{y})\, \mathrm{d}m(\bm{y})\mathrm{d}\delta(s).$ (2.7) $$

During training using the discretized data, each hidden layer is a map from $\mathbb{R}^{d_{v}\times d_{t}\times n\times n}\to\mathbb{R}^{d_{v}\times d_{t}
\times n\times n}$, where $d_{t}$ is a “latent” dimension of time steps that is chosen $\leq n_{t}$.
SFNO proceeds to parametrizes $\kappa(\cdot,\cdot)$ in spacetime the same way FNO3d does $\widetilde{K}_{\phi}\bm{v}=W\bm{v}+\mathcal{F}^{-1}\left(R_{\phi}\cdot(
\mathcal{F}\bm{v})\right)$,
where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the spacetime 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 $(\tau_{\max},k_{\max})$, such that $\{(\tau,\bm{k})\in\mathbb{Z}^{3}:\left\lvert\tau\right\rvert\leq\tau_{\max},
\left\lvert\bm{k}_{j}\right\rvert\leq k_{\max},j=1,2\}$ are kept. $R_{\phi}$ is then parametrized as a $\mathbb{C}^{\tau_{\max}\times k_{\max}\times k_{\max}\times d_{v}\times d_{v}}$-tensor. Here $(\tau,\bm{k})$ denotes the coordinate in the spatiotemporal Fourier domain.
Note that the integral represented by matrix product with $\mathcal{F}\kappa(\cdot,(s,\bm{y}))$ can then be viewed as residing in the continuous space as the Fourier basis (2.8) can be evaluated at arbitrary point.

#### Modifications to the lifting and the projection operator.

FNO3d novelly exploits a convenient architectural advantage of operator-valued NNs: the input temporal dimension $\ell$ is treated as the channel dimension of an image, thus enabling channel mixing as a learnable temporal extrapolation.
Coincidentally, this neat trick makes the lifting operator the biggest hurdle for $\mathcal{G}_{\theta}$ to become an operator learner between Bochner spaces, since $P$’s dimension must be hard-coded (^1^11Line 99, 214, and 217 in [fourier_3d.py at gh/neuraloperator/neuraloperator](https://github.com/neuraloperator/neuraloperator/blob/master/fourier_3d.py#L99)) and thus depends on the input pair’s time steps.
In what follows, we shall present a simple modification to this operator to enable $\mathcal{G}_{\theta}$ in (2.5) to act as an operator that maps an arbitrary-time-step input to an arbitrary-time-step output. We implement a time-depth separable convolution layer $\mathcal{I}$ with variable time steps and periodic padding along the temporal dimension to map an arbitrary number snapshots to a fixed number of channels $d_{v}$.
This layer is largely inspired by the stability and error analysis [, ] for Chorin-Temam projection methods for NSE [, ].
The latent fields, concatenated with the positional encodings (periodically padded), are then composed with a standard pointwise convolution operator $P\in\mathbb{R}^{(d_{v}+3)\times d_{v}}$.
In the SFNO for V-S formulation, the lifting operator is $\widetilde{P}:=P\circ\mathcal{I}:\mathbb{R}^{\ell\times n\times n}\to\mathbb{R
}^{d_{v}\times d_{t}\times n\times n}$ where $d_{t}$ is the latent time steps. The out projection operator $\widetilde{Q}:\mathbb{R}^{d_{v}\times d_{t}\times n\times n}\to\mathbb{R}^{n_{
t}\times n\times n}$ maps the last latent representation to match the dimension of the output $\mathbf{u}_{\operatorname{\mathcal{S}}}$.
In $\widetilde{Q}$, we compose another spectral convolution $K_{L}$ of Fredholm integral of the second kind after channel reduction. $K_{L}$ also pads the temporal frequencies to match the output time steps, as well as mod the constant field in the frequency domain to yield a solution in $H^{1}(\mathbb{T}^{2})/\mathbb{R}$.
These changes are so simple that the new implementations can serve as drop-in replacements for their counterparts in FNO3d when the temporal input dimensions are fixed. For more details, we refer the readers to Appendix D.2.

### 2.3 A posteriori error estimation under a functional norm

We shall present the proposed fine-tuning using the V-P formulation in subsequent subsections.
Denote an SFNO evaluation at a specific $t_{m}\in\mathcal{T}_{n_{t}}:=\{t_{\ell+1},\cdots,t_{\ell+n_{t}}\}$ in the output time interval as $\bm{u}_{\mathcal{N}}^{(m)}$.
The construction of $\widetilde{Q}$ in SFNO renders $\bm{u}_{\mathcal{N}}^{(m)}\in\mathcal{W}$, where $\mathcal{W}$ is the divergence-free subspace of $\operatorname{\mathcal{S}}|_{t=t_{m}}\times\operatorname{\mathcal{S}}|_{t=t_{m
}}\subset\mathcal{V}:=\bm{H}^{1}(\mathbb{T}^{2})$, where

$$ $\operatorname{\mathcal{S}}:=\operatorname{span}\left\{\mathfrak{Re}\left(e^{i \tau_{m}t}e^{i\bm{k}\cdot\bm{x}}\right):\;-{n}/{2}\leq k_{j}\leq{n}/{2}-1,-{n_ {t}}/{2}\leq m\leq n_{t}/2-1\right\}/\mathbb{R},$ (2.8) $$

for $\bm{k}=2\pi(k_{j})_{j=1,2}$ and $\tau_{m}=2\pi m/(T-t_{l})$. Then, a temporal continuous approximation $\bm{u}_{\mathcal{N}}:=\bm{u}_{\mathcal{N}}(t,\cdot)$ can be naturally defined by allowing $t$ vary continuously on $\mathcal{T}:=[t_{\ell+1}-t_{p},T+t_{p}]$ thanks to the spectral basis of $\operatorname{\mathcal{S}}$, where $t_{p}$ is the temporal periodic padding in Section 2.2. Define residual functional $R(\bm{u}_{\mathcal{N}})\in L^{2}(\mathcal{T};\mathcal{V}^{\prime})$: at $t\in\mathcal{T}$ and for $\bm{v}\in\mathcal{V}$

$$ $R(\bm{u}_{\mathcal{N}}):=\bm{f}-\partial_{t}\bm{u}_{\mathcal{N}}-(\bm{u}_{ \mathcal{N}}\cdot\nabla){\bm{u}_{\mathcal{N}}}+\nu\Delta{\bm{u}_{\mathcal{N}}} ,\;\text{ and }\;R(\bm{u}_{\mathcal{N}})(\bm{v}):=\left\langle R(\bm{u}_{ \mathcal{N}}),\bm{v}\right\rangle.$ (2.9) $$

$R(\bm{u}_{\mathcal{N}})$ measures how PDE (2.2) is violated by the finite-dimensional approximation $\bm{u}_{\mathcal{N}}$, not in a localized/pointwise fashion, but rather in a global way by representing the error based on its inner product against arbitrary $\bm{v}$.
At a specific time $t$, the functional norm of $R(\bm{u}_{\mathcal{N}})(t,\cdot)\in\mathcal{V}^{\prime}$ defined as follows is then an excellent measure of the error to be a candidate for a loss function in view of Theorem 2.1:

$$ $\|R(\bm{u}_{\mathcal{N}})(t,\cdot)\|_{\mathcal{V}^{\prime}}:=\sup\nolimits_{ \bm{v}\in\mathcal{V},\|\bm{v}\|_{\mathcal{V}}=1}\left|\left(R(\bm{u}_{\mathcal {N}}),\bm{v}\right)\right|.$ (2.10) $$

###### Theorem 2.1 ( A posteriori error bound for any fine-tuned approximations, informal version) .

Let the weak solution to (2.2) be $\bm{u}\in L^{2}(\mathcal{T};\mathcal{V})$, and $\partial_{t}\bm{u}\in L^{2}(\mathcal{T};\mathcal{V}^{\prime})$. Assume $\bm{u}$ be sufficiently regular, then the dual norm of the residual is efficient to estimate the error for any $\bm{u}_{\mathcal{N}}$:

$$ $\|R(\bm{u}_{\mathcal{N}})\|_{L^{2}(\mathcal{T};\mathcal{V}^{\prime})}^{2} \lesssim\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}_{L^{2}(\mathcal{T};\mathcal{V})}+ \|\partial_{t}(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}_{L^{2}(\mathcal{T};\mathcal{ V}^{\prime})}.$ (2.11) $$

Moreover, if $\bm{u}$ and $\bm{u}_{\mathcal{N}}$ are “sufficiently close”, then it is reliable to serve as an error measure:

$$ $\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}_{L^{\infty}(\mathcal{T}_{m};\mathcal{H})}+ \|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}_{L^{2}(\mathcal{T}_{m};\mathcal{V})}\leq \bigl{\|}(\bm{u}-\bm{u}_{\mathcal{N}})(t_{m},\cdot)\bigr{\|}^{2}_{\mathcal{V}} +C\int_{\mathcal{T}_{m}}\|R(\bm{u}_{\mathcal{N}})(t,\cdot)\|_{\mathcal{V}^{ \prime}}^{2}\,\mathrm{d}t.$ (2.12) $$

where $\mathcal{T}_{m}:=(t_{m},t_{m+1}]$, and the constants depend on the regularity of the true solution $\bm{u}$.

###### Theorem 2.1 ( A posteriori error bound for any fine-tuned approximations, informal version) .

### 2.4 Fine-tuning using negative Sobolev norm as loss function

The functional norm (2.10) is “global” because it does not have a natural $\ell^{2}$-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., see [] and [, Chapter 3]), we can define the pairing between $\mathcal{V}$ and $\mathcal{V}^{\prime}$ as follows without getting into the intricate natures of the tempered distribution:

$$ $\langle f,g\rangle_{\mathcal{V},\mathcal{V}^{\prime}}\;\text{``}\!=\!\text{''} \int_{\mathbb{Z}^{2}}(1+|\bm{k}|^{2})^{-s}\hat{f}(\bm{k})\overline{\hat{g}(\bm {k})}(1+|\bm{k}|^{2})^{s}\mathrm{d}\bm{k},\text{ where }\hat{v}:=\mathcal{F}v \text{ for }v\in\mathcal{V}^{\prime}.$ (2.13) $$

The spatial Sobolev space $H^{s}(\mathbb{T}^{2})$ can be alternatively identified using norm $\|\cdot\|_{s}$ and seminorm $|\cdot|_{s}$ (e.g., see [, Def. 3.2.2]) as follows for any $s\in\mathbb{R}$

$$ $\|f\|_{s}^{2}:=\sum_{\bm{k}\in\mathbb{Z}^{2}}(1+|\bm{k}|^{2})^{s}|\hat{f}(\bm{ k})|^{2},\;\text{ and }\;|f|_{s}^{2}:=\sum_{\bm{k}\in\mathbb{Z}_{n}^{2} \backslash\{\bm{0}\}}|\bm{k}|^{2s}\bigl{|}\hat{f}(\bm{k})\bigr{|}^{2},\;s\neq 0.$ (2.14) $$

Moreover, we have the subsequent lemma in our specific case.

###### Theorem 2.2 (Functional norm “ $\simeq$ ” negative norm) .

If $f\in\mathcal{H}:=L^{2}(\mathbb{T}^{2})/\mathbb{R}$, $\|f\|_{\mathcal{H}^{\prime}}=|f|_{-1}$.

Spatially, we realize a regularized negative Sobolev norm by the Fast Fourier Transform (FFT):

$$ $\|f\|_{-1,\alpha,n}^{2}:=\sum\nolimits_{\bm{k}\in\mathbb{Z}_{n}^{2}\backslash \{\bm{0}\}}(\alpha+|\bm{k}|^{2})^{-1}\bigl{|}\hat{f}(\bm{k})\bigr{|}^{2},\text { where }\mathbb{Z}_{n}^{2}:=(\mathbb{Z}/n\mathbb{Z})^{2},\text{ and }\alpha \geq 0.$ (2.15) $$

With this norm at hand, (2.10) and the Bochner norm in (2.11) of the residual are realized to serve as the loss function in the fine-tuning

$$ $\eta_{m}(\bm{u}_{\mathcal{N}},\partial_{t}\bm{u}_{\mathcal{N}}):=\|R(\bm{u}_{ \mathcal{N}})(t_{m},\cdot)\|_{-1,\alpha,n}.$ (2.16) $$

###### Theorem 2.2 (Functional norm “ ≃ similar-to-or-equals \simeq ≃ ” negative norm) .

#### Parseval identity.

In the context of using an optimization algorithm to train an NN as a function approximator, it is known (e.g., []) that whether the integral in the loss function is accurately computed affects whether the NN can achieve the scientific computing level of accuracy.
For example, on a uniform grid, the accuracy of the mesh-weighted spatial $\ell^{2}$-norm as the numerical quadrature is highly affected by the *local* smoothness of the function in consideration.
Nevertheless, computing the integral in the frequency domain yields exponentially convergent approximations [, Theorem 3.1] thanks to the Parseval identity.

#### Why functional-type norm for the residual evaluation?

We note that in “physics-informed” learning of operators, for example, [], the PDE residual is evaluated using $L^{2}$-norm as loss. In the meantime, *positive* Sobolev norm, which is local in terms of differential operators, is used in []. To our best knowledge, the functional norm has not been applied in either function learning problems, such as PINN [], or operator learning problems. In fact, the relation of the Gelfand triple is so simple and elegant, leading to Theorem 2.2, that the functional norm is nothing but an inverse frequency weight in 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 high-frequency part, it has been discovered in [] that the error of operator learning for NSE is still dominant in the lower end of the spectrum. This has been corroborated in our study as well, see Figure 5 and Figure 4. For a more detailed and mathematically enriched discussion on why functional norm is not widely popular in traditional numerical methods, we refer the reader to Appendix B.

### 2.5 New paradigm

With the SFNO and the error estimators, we propose a new training-fine-tuning paradigm. Another important motivation of this new paradigm is that, experimentally, either SFNO or FNO3d captures the statistical properties of the 2D turbulence [, ] quickly in training. Even after a single epoch, the evaluations of SFNO converge to a tight neighborhood of the ground truth in terms of the frequency signature of the inverse cascade of Kolmogorov flow (for details, please refer to Appendix D.3). Then, $95\%$ of the FLOPs spent in training has marginal improvements, which motivates us to rethink a more efficient paradigm than end-to-end, in the meantime not needing to initiate the expensive training of another nonlinear corrector as PDE-refiner [] in the postprocessing phase.

#### 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 (2.9),
one needs to compute $\partial_{t}\bm{u}_{\mathcal{N}}$.
Instead of a common approach of using NO to represent the extra field variables [, ], we opt to apply *a traditional numerical solver $\mathcal{G}_{\alpha}(\cdot)$ with an extremely fine time step*, e.g., $\mathcal{O}(\Delta t^{\alpha})$ for $\alpha\geq 2$, to compute these extra field variables (Line 10 of Algorithm 1), while preserving the computational graph for auto-differentiation. We note that, in NSE simulations using traditional numerical solvers, for efficiency, the magnitude of this fine time step is never realistic or attainable due to time marching. Given the training data with trajectories at $\{t_{1},\dots,t_{{\ell}}\}$ aiming to predict the trajectories at $\{t_{\ell+1},\dots,t_{\ell+n_{t}}\}$, the new paradigm for SFNO is outlined in Algorithm 1.

Figure: Algorithm 1 The new parallel-in-time training-fine-tuning strategy in small data regime.

#### Interpretations of the theoretical results.

Theorem 2.1 establishes that the functional norm of the residual $R(\bm{u}_{\mathcal{N}})$, without accessing $\bm{u}$, is a good representation of the error $\bm{u}-\bm{u}_{\mathcal{N}}$ in Bochner norms. Estimate (2.11) indicates that reducing the *a posteriori* error estimator is a *necessary* condition for reducing the true error. While estimate (2.12) is more delicate in that it is only reliable when $\bm{u}_{\mathcal{N}}$ gets “close” to $\bm{u}$. Theorem 2.2 lays the foundation to accurately evaluate this functional norm via FFT. Nevertheless, (2.12) serves as a guideline to design the “refining” procedure (lines 13 in Algorithm 1), where the error estimator moves to refine the next time step once the error in the previous one becomes less than a given threshold.

## 3 Numerical Experiments

### 3.1 Illustrative example: Taylor-Green vortex

In this illustrative example, we examine the 2D Taylor-Green vortex model [], whose analytical solution is *known* and frequently employed as a benchmark for evaluating traditional numerical schemes for NSE []. We create a toy train dataset with 10 trajectories on a $256^{2}$ grid with varying numbers of vortices per wavelength, and the test sample has an unseen number of vortices yet still can be resolved up to the Nyquist scale. For details please refer to Appendix C.2. The sizes FNO3d and SFNO models are also scaled down with $L=1$ in (2.5).

**Table 1: Results for Taylor-Green vortex example $\varepsilon:=\omega_{\text{true}}-\omega_{\mathcal{N}}$, the errors are reported at the final time step.**
|  | Evaluation after training | After fine-tuning |  |  |
| --- | --- | --- | --- | --- |
|  | $\|{\varepsilon}\|_{L^{2}}$ | $\|R\|_{-1,n}$ | $\|{\varepsilon}\|_{L^{2}}$ | $\|R\|_{-1,n}$ |
| FNO3d | $1.84\text{\times}{10}^{-01}$ | $5.40\text{\times}{10}^{-01}$ | N/A | N/A |
| SFNO | $1.94\text{\times}{10}^{-01}$ | $2.18\text{\times}{10}^{-01}$ | $1.24\text{\times}{10}^{-07}$ | $3.21\text{\times}{10}^{-07}$ |
| PS+RK2 (GT) | $5.91\text{\times}{10}^{-6}$ | $1.16\text{\times}{10}^{-05}$ | N/A | N/A |

**Table 2: Results for forced turbulence, small numbers of vortices. $\varepsilon:=\omega_{\operatorname{\mathcal{S}}}-\omega_{\mathcal{N}}$, original example from [].**
|  | Evaluation after training | After fine-tuning |  |  |
| --- | --- | --- | --- | --- |
|  | $\|{\varepsilon}\|_{L^{2}}$ | $\|R\|_{-1,n}$ | $\|{\varepsilon}\|_{L^{2}}$ | $\|R\|_{-1,n}$ |
| FNO3d 10 ep | $1.31\text{\times}{10}^{-2}$ | $1.30\text{\times}{10}^{-2}$ | N/A | N/A |
| SFNO 10 ep + $L^{2}$ FT | $2.08\text{\times}{10}^{-02}$ | $1.27\text{\times}{10}^{-02}$ | $2.82\text{\times}{10}^{-04}$ | $2.78\text{\times}{10}^{-05}$ |
| SFNO 10 ep + $H^{-1}$ FT | – | – | $3.16\text{\times}{10}^{-04}$ | $4.59\text{\times}{10}^{-07}$ |

### 3.2 2D Kolmogorov turbulence

This meta-example contains a series of examples of isotropic turbulence [], featured in various work such as [, ]. The energy and the enstrophy spectra satisfy the energy law of turbulence first proposed by A. N. Kolmogorov []. The data are generated using a second-order time-stepping scheme that is proven in theory to preserve the inverse cascade of the energy/enstrophy spectra [, ]. We consider two cases:

- (I)
FNO3d benchmark, $\nu=10^{-3}$, $\omega_{0}\sim\mathcal{N}(0,(-\Delta+\tau^{2}I)^{-\alpha/2})$, the energy density in wavenumber $k:=|\bm{k}|$ is $f(k)\sim(k^{2}+\tau^{2})^{-\alpha}$, no drag, fixed force with low wavenumber.
- (II)
The famous decaying turbulence discovered by the McWilliams [], the initial power spectrum is chosen that the decaying is slow and the enstrophy density evolves to the inverse cascade of Kolmogorov flows featured in [, , ].

Tables 2 and 3 report the results for example (I) and example (II), respectively.

**Table 3: Evaluation metrics of the McWilliams isotropic turbulence example. All models are trained using on $64\times 64$ mesh and evaluated on $256\times 256$ mesh. $\varepsilon:=\bm{u}_{\operatorname{\mathcal{S}}}-\bm{u}_{\mathcal{N}}$ or $\omega_{\operatorname{\mathcal{S}}}-\omega_{\mathcal{N}}$. $\operatorname{\mathcal{S}}$ stands for the Fourier approximation space (2.8) on a $256\times 256$ fine grid. $r:=R(\bm{u}_{\mathcal{N}})$ or $R(\omega_{\mathcal{N}})$. $\mathcal{Y}:=H^{-1}(\mathbb{T}^{2})$. $\|\cdot\|_{*}$: Bochner norm with convection. All error norms are evaluated at 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 in line 12.**
|  | Evaluation after train | After fine-tuning |  |  |
| --- | --- | --- | --- | --- |
|  | $\|{\varepsilon}\|_{L^{2}}$ | $\|R\|_{-1,n}$ | $\|{\varepsilon}\|_{L^{2}}$ | $\|R\|_{-1,n}$ |
| FNO3d 10 epochs $+L^{2}$ train | $2.44\text{\times}{10}^{-02}$ | $6.09\text{\times}{10}^{-02}$ | N/A | N/A |
| FNO3d 10 epochs $+H^{-1}$ train | $6.52\text{\times}{10}^{-02}$ | $2.78\text{\times}{10}^{-01}$ | N/A | N/A |
| FNO3d 100 ep $+L^{2}$ train | $4.32\text{\times}{10}^{-02}$ | $2.44\text{\times}{10}^{-02}$ | N/A | N/A |
| SFNO 10 epochs $+H^{-1}$ train & $H^{-1}$ FT | $6.54\text{\times}{10}^{-02}$ | $6.19\text{\times}{10}^{-02}$ | $5.71\text{\times}{10}^{-03}$ | $2.24\text{\times}{10}^{-05}$ |
| SFNO 10 epochs $+L^{2}$ train & FT | $3.69\text{\times}{10}^{-02}$ | $2.35\text{\times}{10}^{-02}$ | $1.79\text{\times}{10}^{-03}$ | $9.55\text{\times}{10}^{-05}$ |
| SFNO 10 epochs $+L^{2}$ train & $H^{-1}$ FT | – | – | $2.88\text{\times}{10}^{-04}$ | $4.02\text{\times}{10}^{-06}$ |

## 4 Conclusion and Limitations

We designed a new pipeline to train and fine-tune an FNO variant to get close to the true solution (not ground truth generated by numerical solver) of NSE. The evaluation errors in benchmark problems are up to $10^{5}$ times better than FNO3d trained by a simple end-to-end pipeline. However, we acknowledge that our method is currently limited to non-factorized parametrized FNO, since many intricate connections with traditional numerical methods are exploited, e.g., the optimally learned parameters correspond to a Fourier-Galerkin projection []. How to generalize this new pipeline to other types of operator learners in Appendix B will be worthy of exploration.

## Acknowledgments and Disclosure of Funding

The research of Brarda and Xi is supported by the National Science Foundation award 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 DEAC52-07NA27344 and was supported by the LLNL-LDRD program under Project No. 24ERD033. Cao is supported in part by the National Science Foundation award DMS-2309778.

## Appendix A Notations

**Table 4: Notations used in an approximate chronological order and their meaning in this work.**
| Notation | Meaning |
| --- | --- |
| $\Omega$ | $\Omega=(0,2\pi)^{2}$ or $(0,1)^{2}$ the modeling domain in $\mathbb{R}^{2}$ |
| $\mathcal{V},\mathcal{H}$ | Hilbert spaces defined on the domain $\Omega$ above,<br>$f\in\mathcal{H}:\Omega\to\mathbb{R}$ |
| $\mathcal{X},\mathcal{Y}$ | product spaces $\mathop{\mathchoice{\textstyle\prod}{\prod}{\prod}{\prod}}\nolimits_{j\in<br>\Lambda}\mathcal{H}$ |
| $||$ | concatenation, for $f_{j}\in\mathcal{H}$, $||_{j\in\Lambda}f_{j}\in\mathop{\mathchoice{\textstyle\prod}{\prod}{\prod}{<br>\prod}}\nolimits_{j\in\Lambda}\mathcal{H}$ |
| $\nu$ | viscosity, the inverse of the Reynolds number, strength of diffusion |
| $\mathbb{T}^{2}$ | $\mathbb{T}^{2}:=\mathbb{R}^{2}/\mathbb{Z}^{2}$. |
| $H^{1}(\mathbb{T}^{2})$ | $\{u\in H^{1}(\Omega):\,u\text{ satisfies the periodical boundary condition}\}$ |
| $\bm{H}^{1}(\mathbb{T}^{2})$ | $\bm{H}^{1}(\mathbb{T}^{2}):=H^{1}(\mathbb{T}^{2})\times H^{1}(\mathbb{T}^{2})$ |
| $L^{p}(\mathcal{T};\mathcal{V})$ | the Bochner spaces containing $\{u:\mathcal{T}\!\to\!\mathcal{V}\,\big{|}\int_{\mathcal{T}}\|u(t,\cdot)\|_{<br>\mathcal{V}}^{p}\,\mathrm{d}t<+\infty\}$, $p\in[1,\infty)$ |
| $L^{\infty}(\mathcal{T};\mathcal{V})$ | Bochner space containing<br>$\{u:\mathcal{T}\!\to\!\mathcal{H}\,\big{|}\operatorname{ess\,sup}_{t\in<br>\mathcal{T}}\|u(t,\cdot)\|_{\mathcal{V}}<+\infty\}$ |
| $\bm{u}$ | the true weak solution to the NSE $\bm{u}(t,\cdot)\in\mathcal{V}$ for any $t$ |
| $\bm{u}_{\operatorname{\mathcal{S}}}^{(l)}$ | the ground truth data generated at $t_{l}$ by the numerical solver in $\operatorname{\mathcal{S}}\subset\mathcal{V}$ |
| $\bm{u}_{\mathcal{N}}^{(m)}$ | neural operator evaluations at $t_{m}$ that can be embedded in $\operatorname{\mathcal{S}}\subset\mathcal{V}$ |
| $(u,v)$ or $(\bm{u},\bm{v})$ | the $L^{2}$-inner product on $\mathcal{H}$, $(u,v):=\int_{\Omega}uv\,\mathrm{d}\bm{x}$ |
| $\|u\|_{s}$ | the $H^{s}$-norm of $u$, computed by (2.14), $\|u\|:=\|u\|_{0}$ falls back to $L^{2}$-norm |
| $|u|_{s}$ | the $H^{s}$-seminorm of $u$, computed by (2.14), is a norm on $H^{s}(\mathbb{T}^{2})/\mathbb{R}$ |
| $\lesssim$ | $a\lesssim b$ means that $\exists c$ independent of asymptotics if any such that $a\leq cb$ |
| $\simeq$ | $a\simeq b\Leftrightarrow a\lesssim b$ and $b\lesssim a$ |
| $\mathcal{V}^{\prime}$ | dual space of $\mathcal{V}$, contains all continuous functionals $f$ such that $|f(v)|\lesssim\|v\|_{\mathcal{V}}$ |
| $\langle f,v\rangle$ | the pairing between $v\in\mathcal{V}$ and $f\in\mathcal{V}^{\prime}$, can be identified as $(f,v)$ if $f\in\mathcal{H}$ |
| $\mathcal{V}\hookrightarrow\mathcal{H}$ | $\mathcal{V}$ is continuously embedded in $\mathcal{H}$ such that $\|u\|_{\mathcal{H}}\lesssim\|u\|_{\mathcal{V}}$ |
| $\mathcal{V}\Subset\mathcal{H}\Subset\mathcal{V}^{\prime}$ | compact embeddings by Poincaré inequality, $\mathcal{H}=L^{2}(\mathbb{T}^{2})$, $\mathcal{V}=H^{1}(\mathbb{T}^{2})$ |
| $\bm{A}:\bm{B}$ | $\bm{A}:\bm{B}=\sum_{1\leq i,j\leq 2}A_{ij}B_{ij}$ for $\bm{A},\bm{B}\in\mathbb{R}^{2\times 2}$ |
| $\bm{a}\otimes\bm{b}$ | $\bm{a}\bm{b}^{\top}$ if both are viewed as column vectors |

## Appendix B Detailed Literature Review and Motivations

#### Interplay of deep learning and PDEs.

PDE solvers are function learners to represent PDE solutions using neural networks [, , ].
PDE discovery encompasses all the techniques dedicated to identifying and optimizing PDE coefficients from data [, ]. Recently, hybrid solver approach has been explored in [, , , , , , ]. Reinforcement learning has been applied in the field of mesh generation to build more efficient solving pipelines [, ]. For neural operators, remarkable outcomes are achieved in diverse applications, such as weather forecasting [] and turbulent fluids simulations [, , ], 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, a feature very appealing in many engineering applications []; and the ability to provide mesh-free and resolution-independent solvers in cases of irregular domains [, ].

#### Neural operators.

Looking at neural operator architectures, we can identify two different approaches. 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) [, , ], wherein aggregation occurs via linear combination. Similarly, Fourier Neural Operator (FNO) [], and some of its variants [, , , , ], aggregate different latent representations through convolution with a learnable kernel in the frequency domain.
Additional examples with a linear aggregation with a U-Net meta-architecture can be found in [, ]. [] further explores the error analyses for the linear aggregations through the lens of the frame(let) theory.
The second approach to designing neural operator architectures aims to obtain the “basis” through linear projections of the latent representations.
In this scenario, the aggregation is nonlinear, for example, using a signal-dependent kernel integral. In this category, the most notable example is the scaled dot-product attention operator in Transformer [], which builds an instance-dependent kernel. In the context of operator learning applications, Transformer-based neural operatos have been studied in [, , , , , ]. It is also shown in [] that nonlinear aggregations outperform its linear counterparts in learning solutions with less regularity.

#### Numerical methods for NSE.

The NSE falls into the category of a stiff PDE (system) $\partial_{t}u=Lu+N(u)+f$, where $f$ is the external forcing, $L$ and $N(\cdot)$ are a linear and a nonlinear operator, respectively. Trefethen noted back in [] on the difficulty to design a time-stepping scheme as $N(\cdot)$ 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 [, ]. Nonlinear Galerkin method [] inspires us to prove Theorem E.9. [] uses a clever trick to impose the divergence free condition without constructing a divergence-conforming finite element subspace.
[] 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., [] for a summary), also serves as an inspiration to add the Helmholtz layer. [] proposed a mixed discretization for the vorticity-streamfunction formulation.

#### Consistency-stability trade-offs.

In view of the Lax equivalence principle (“consistency” $+$ “stability” $\implies$ “convergence” [, Theorem 8]), the improvement of the stability of the method has a trade-off with a method’s 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 [], 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 $\mathcal{O}(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.

#### De-aliasing filter sacrifices consistency for stability.

Numerical results suggest that for pseudo-spectral spatial discretization with no higher-order Fourier smoothing temporally, the dealiasing filter is indispensable [], as the time marching may experience numerical instability [, ] without it. This is due to the nonlinear interaction in the convective term, which is caused by the amplification of high-frequency “aliasing” errors when the underlying solution lacks sufficient smoothness.

#### Why and why not functional-type a posteriori error estimation?

Speaking of the *a posteriori* error 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 $L^{2}$ residuals to indicate where the mesh needs to be refined.
There are various compromises for this $H^{-1}$-to-$L^{2}$ representation to happen that renders the estimate inaccurate, such as discrete Poincaré constant [] or inverse inequalities [, ], see also [, \origS1.6.2]. Functional-type *a posteriori* error estimation [] consider the error as a functional, which is equivalent to error to bilinear form-associated norm as follows

$$ $R(u_{\operatorname{\mathcal{S}}})\in\mathcal{V}^{\prime},\text{ and }\langle R (u_{\operatorname{\mathcal{S}}}),v\rangle=(u-u_{\operatorname{\mathcal{S}}},v) _{a},$ $$

where the weak form of the PDE is

$$ $(u,v)_{a}=f(v)\;\forall v\in\mathcal{V}\text{ and }(u_{\operatorname{\mathcal{ S}}},v)_{a}=f(v)\;\forall v\in\operatorname{\mathcal{S}}\subset\mathcal{V}.$ $$

The common approach is using the help from extra “flux” or “stress” dual variables ([, \origS6.4], see also []) for how to get an accurate representation of the error functional in $L^{2}$. For example, for the Stokes flow (steady-state viscous fluid, letting $\partial_{t}\bm{u}=0$ and no nonlinear convection in (2.2), $\nu\sim\mathcal{O}(1)$), [, \origS6.2] estimates error of an $\bm{H}(\mathrm{div})$–$L^{2}$ mixed discretization as follows

$$ $\nu\|\nabla(\bm{u}-\bm{u}_{\operatorname{\mathcal{S}}})\|\leq\|\bm{\sigma}+q \bm{I}-\nu\nabla\bm{u}_{\operatorname{\mathcal{S}}}\|+C_{P}\|\nabla\cdot\bm{ \sigma}+\bm{f}\|,$ $$

where $C_{P}$ is the Poincaré constant of the compact embedding and $(\bm{\sigma},q)$ is a reconstruction field pair. However, the drawback of this approach is that an expensive global minimization problem needs to be solved, e.g., for $(\bm{\sigma},q)$ above. Another main reason to introduce extra field variables is that, for finite element methods, the basis functions are local and do not have a globally continuous derivative, in that $\Delta\bm{u}_{\operatorname{\mathcal{S}}}$ yields singular distributions, whose proper norm is $H^{-1}$-norm and cannot be evaluated by summing up element-wise $L^{2}$-norms. In computation, it has to be replace by $\nabla\cdot\bm{\sigma}$ where $\bm{\sigma}\in\bm{H}(\mathrm{div})$, and is constructed to be closer to the true solution’s gradient $\nabla\bm{u}$ than $\nabla\bm{u}_{\operatorname{\mathcal{S}}}$.
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 to be tweaked to ensure the Ladyzhenskaya-Babuška-Brezzi stability condition [, Chapter III \origS1 Sec 1.2].
Recently, [] considers NN as a function learner to represent solutions of linear convection-diffusion 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 $\mathcal{O}(10^{-3})$. In contrast, in our study, the need of extra “stress” or “flux” variables to build the residual functional is circumvented as well.

#### A posteriori error estimation for flow problems

There is a vast amount of literature for the *a posteriori* error estimation for the viscous flow ($1/\text{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 is from [], and it is of residual-type by computing a mesh-weighted $L^{2}$-norm elementwise, and the singular distribution $\Delta\bm{u}_{\operatorname{\mathcal{S}}}$ is represented by the magnitude of flux jumps across the facets in this mesh. In [], a more accurate *a posteriori* estimation technique is invented for Stokes problem in which a local problem is solved to represent the residual functional on a collection of neighboring elements. $L^{2}$ residual-type error estimation for stationary NSE is considered in []. To our best knowledge, no functional-type *a posteriori* error estimation has been applied to solve the transient NSE, due to 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 marching solvers. Notable examples include: solution from Direct Numerical Simulation (DNS) corrected using an NN [];
interpolation to achieve spatial super-resolution using solutions produced by traditional numerical solvers marching on coarser grids [].
orthogonal polynomial features (HiPPO from []) to achieve temporal super-resolution in addition to the spatial one []. [] explores how to incorporate energy-preserving schemes such as upwinding into the ML solvers through traditional numerical methods.
[] identifies important frequency bands in FNO parametrization and improves FNO’s rolling-out performance in NSE benchmarks.

## Appendix C Data Generation

### C.1 Vorticity–Streamfunction Formulation

Here we derive the vorticity-streamfunction formulation for the convenience of the readers. Consider the standard NSE in 3D, $\bm{u}$’s $z$-component is 0, and $\bm{u}=\bm{u}(t,x,y)$ has only planar dependence

$$ $\left\{\begin{aligned} \dfrac{\partial\bm{u}}{\partial t}+(\bm{u}\cdot\nabla) \bm{u}&=-\dfrac{1}{\rho}\nabla p+\nu\Delta\bm{u}+\bm{f}.\\ \hfill\nabla\cdot\bm{u}&=0\end{aligned}\right.$ (C.1) $$

Taking $\nabla\!\times\!(\cdot)$ 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 $\bm{\omega}=\nabla\times\bm{u}$

$$ $\frac{\partial\bm{\omega}}{\partial t}+\nabla\times\nabla\left(\frac{\bm{u}^{2 }}{2}\right)+\bm{u}\cdot\nabla\bm{\omega}-\bm{\omega}\cdot\nabla\bm{u}=-\nabla \times\left(\frac{1}{\rho}\nabla p\right)+\nu\Delta\bm{\omega}+\nabla\!\times \!\bm{f}.$ (C.2) $$

Next, one can have $\bm{\omega}\cdot\nabla\bm{u}=0$, thus the equation becomes

$$ $\frac{\partial\bm{\omega}}{\partial t}+\bm{u}\cdot\nabla\bm{\omega}=\bm{\omega }\cdot\nabla\bm{u}+\nu\Delta\bm{\omega}+\nabla\!\times\!\bm{f}.$ (C.3) $$

We need to introduce a streamfunction $\psi$ to recover the velocity $\bm{u}=\nabla\!\times\!(0,0,\psi)$. The main simplification comes from the assumption of dependence on $x$- and $y$-variable only. Thus, in 2D this becomes $\bm{u}=\operatorname{\mathbf{rot}}\psi$ (a $\pi/2$ rotation of the 2D gradient $\nabla\psi$). Therefore,

$$ $\bm{u}=(u_{1},u_{2},0)\implies u_{1}=\partial_{y}\psi,u_{2}=-\partial_{x}\psi$ $$

and the vector vorticity can be equivalently represented by a scalar vorticity $\omega$

$$ $\bm{\omega}=\nabla\!\times\!\bm{u}=(0,0,\operatorname{curl}(u_{1},u_{2}))=(0,0 ,\partial_{x}u_{2}-\partial_{y}u_{1})=:(0,0,\omega)$ $$

The equation becomes a nonlinear system for $\omega$ and $\psi$:

$$ $\displaystyle\partial_{t}\omega+\operatorname{\mathbf{rot}}\psi\cdot\nabla{ \omega}-\nu\Delta{\omega}$ $\displaystyle=\nabla\!\times\!\bm{f},$ (C.4) $\displaystyle\omega+\Delta\psi$ $\displaystyle=0.$ (C.5) $$

### C.2 Taylor-Green Vortex

Taylor-Green vortex serves as one of the most well-known benchmarks for NSE numerical methods, as its flow type experiences from laminar, to transitional, and finally evolving into turbulent regime. We consider the trajectory before the breakdown phase. We opt to use a doubly periodic solution on $[0,2\pi)^{2}$ such that the inflow/outflow occurs on the “boundary” (see Figure 2). The exact solution is given by:

$$ $\displaystyle\bm{u}(t,x,y)=e^{-2\kappa^{2}\nu t}\begin{pmatrix}\sin(\kappa x) \cos(\kappa y)\\ -\cos(\kappa x)\sin(\kappa y)\end{pmatrix}\text{ or }e^{-2\kappa^{2}\nu t} \begin{pmatrix}-\cos(\kappa x)\sin(\kappa y)\\ \sin(\kappa x)\cos(\kappa y),\end{pmatrix}$ $$

where $\nu=1/\text{Re}$ is chosen as $10^{-3}$.
For a sample trajectory with $\kappa=1$, please refer to Figure 2.
We apply the pseudo-spectral method with RK2 time stepping for the convection term and Crank-Nicolson for the diffusion term. The dataset has 11 trajectories with $\kappa$ ranging from $1$ to $11$, among which $\kappa=11$ is chosen as the test trajectory

Figure: Figure 2: Ground truth streamlines for Taylor-Green vortex example.
Refer to caption: x2.png

### C.3 Kolmogorov flow

The example featured in the original FNO paper [] can be viewed a strongly-regularized initial condition, no drag, force with low wavenumber. Our example (II) has weakly-regularized initial condition, small drag, force with mid wavenumber. This resembles the example featured in the Jax-CFD paper [].

Figure: (a)
Refer to caption: x3.png

Figure: (a)
Refer to caption: x6.png

## Appendix D Experiments

### D.1 Training and evaluation

The training uses torch.optim.OneCycleLR [] learning rate strategy with a $20\%$ warm-up phase. AdamW is the optimizer with no extra regularization. The learning rate starts and ends with $10^{-3}\cdot lr_{\max}$. The $lr_{\max}=10^{-3}$. The result demonstrated is obtained from fixing the random number generator seed.
All models are trained on a single RTX A4500 or RTX A6000. The codes to replicate the experiments are open-source and publicly available.
(^2^22Please download the code in the anonymous repository and follow the instructions in README.md.)

Despite using the negative Sobolev norm is quite efficient in fine-tuning, using it in training not be the most efficient in minimizing the $L^{2}$-norm due to the optimization being nonconvexity, we observe that all norms are equivalent “bad” due to nonlinearity of NOs.

### D.2 Models

#### Helmholtz layer for the velocity-pressure formulation.

For the V-P formulation in a simplified connected convex or periodic domain, it is known that an exact divergence-free subspace $\mathcal{W}\subset\mathcal{V}$ for velocity means that the pressure field is not needed. The reason is that in the weak formulation, the pressure is a Lagrange multiplier to impose the divergence-free condition [, Chapter III 1 Section 1.1].
Inspired by the postprocessing to eliminate $\nabla p$ [, eq. (15)] together with the neural Clifford layers [], we add a Helmholtz layer $S$ after each application of $\sigma_{j}\circ(W_{j}+K_{j})$. $S$ performs a discrete Helmholtz decomposition [, Chapter 1 \origS1 Section 3.1] in the frequency domain to project the latent fields to be solenoidal. For details on the implementation, please refer to Appendix D.2.

#### Difference with Fourier Neural Operator 3d

In view of (2.6), the FFTs in FNO3d transforms are continuous integrals in the spatial dimensions yet a discrete sum in the temporal dimension. Despite that no explicit restriction imposed on the output time steps, the $d_{\text{out}}$ cannot be trivially changed as the data are prepared by applying a pointwise Gaussian normalizer that depends on $d_{\text{out}}$.
Denote the lifting operator (channel expansion) in FNO3d by $P:\{\mathbf{a}_{h}\}\oplus\{\mathbf{p}_{h}\}\to\mathbb{R}^{d_{v}\times d_{
\text{out}}\times n\times n}$. Before the application of $P$ in FNO3d, $\mathbf{a}_{h}$ is artificially repeated $d_{\text{out}}$ times, then $P$ is in a space of linear operators $\simeq\mathbb{R}^{(d_{\text{in}}+3)\times d_{v}}$. This observation makes it independent of the spatial discretization size $n\times n$ by dependent on $d_{\text{in}}$. In FNO3d, the channel reduction is a pointwise nonlinear universal approximator, yet in SFNO, this is linear.

SFNO temporal dimension is a super-resolution, while FNO produces trajectories in future time steps.

The choice of latent time marching draws inspiration from the inverted bottleneck architecture in Transformers and ConvNeXt [] to select the temporal global basis.

**Table 5: The detailed comparison of FNO3d and SFNO. Layer: # of spectral convolution layers; SFNO has an extra layer of spectral convolution in a single channel with skip connection. Modes: $(\tau_{\max},k_{\max})$. 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 a $256^{2}$ grid for a single instance. A torch.cfloat type parameter entry counts as two parameters. Models are not profiled for Example 3.1 as they are not full-fledged models.**
|  | Architectures | GFLOPs | # params |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | layers | channel/width | modes | activation | pre-norm | Eval | FT |  |
| FNO3d Example 1 | $1$ | 4 | $(2,64)$ | ReLU | N | – | – | 1.05m |
| SFNO Example 1 | $1+\frac{1}{4}$ | 4 | $(2,64)$ | ReLU | N | – | – | 1.11m |
| FNO3d Example 2 (I) | $4$ | 20 | $(5,8)$ | GELU | Y | $13.3$ | N/A | 9.03m |
| SFNO Example 2 (I) | $4+\frac{1}{20}$ | 20 | $(5,8)$ | GELU | N | $27.0$ | $3.28$ | 9.02m |
| FNO3d Example 2 (II) | $4$ | 10 | $(5,32)$ | GELU | Y | 28.7 | N/A | 16.38m |
| SFNO Example 2 (II) | $4+\frac{1}{10}$ | 10 | $(5,32)$ | GELU | N | 42.5 | $3.4$ | 16.42m |

### D.3 Fine-tuning

In implementations, we choose $\texttt{Iter}_{\max}=100$ in Algorithm 1, i.e., an ADAM is simply run for 100 iterations to update parameters that count only a fraction of a spectral layer since there is only 1 channel of weights.

#### To train or to stop? Investigating the frequency domain signatures

In Stuart’s paper on the convergence of linear operator learning [], 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.

Denote $\omega(t,\bm{x}):=\nabla\!\times\!\bm{u}$, and $\hat{\omega}:=\mathcal{F}\omega$ where $\mathcal{F}$ is applied only in space, then we can define enstrophy spectra as follows:

$$ $\mathcal{E}(t,k)=\sum_{k-\delta k\leq|\bm{k}|\leq k+\delta k}|\hat{\omega}(t, \bm{k})|^{2},\;\text{ and }\;\sum_{k}\mathcal{E}(t,k):=\int_{\Omega}|\omega|^{ 2}\mathrm{d}\bm{x}$ $$

The enstrophy spectrum which decays slightly faster than kinetic energy but still behaves like an inverse cascade. For the original FNO3d, it learns the frequency signature of the data after a single epoch. See Figure 5.

Figure: (a)
Refer to caption: x8.png

## Appendix E Assumptions and Proofs

###### Assumption E.1 (Assumptions for Theorems 2.1 , 2.2 , E.9 ) .

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.

- ($E_{1}$)
The compact embeddings for Gelfand triple $\mathcal{V}\Subset\mathcal{H}\Subset\mathcal{V}^{\prime}$ hold, where $\mathcal{V}=\bm{H}^{1}(\mathbb{T}^{2})$, and $\mathcal{H}=\bm{L}^{2}(\mathbb{T}^{2})$.
- ($E_{2}$)
The time interval $\mathcal{T}$ in consideration is fixed in that we consider the “refining” problem for a fixed time interval but with variable time steps.
- ($E_{3}$)
The initial condition $\bm{u}_{0}\in\mathcal{V}$, and $\bm{f}\in L^{2}(\mathcal{T};\mathcal{V}^{\prime})$.

###### Lemma E.2 (Skew-symmetry of the trilinear term) .

For $\bm{z},\bm{u},\bm{v}\in\bm{H}^{1}(\mathbb{T}^{2})$, where $\mathbb{T}^{2}$ denotes $\Omega:=(0,1)^{2}$ equipped with component-wise PBC. Denote

$$ $c(\bm{z},\bm{u},\bm{v}):=\bigl{(}(\bm{z}\cdot\nabla)\bm{u},\bm{v}\bigr{)}=\int _{\Omega}((\bm{z}\cdot\nabla)\bm{u})\cdot\bm{v}\,\mathrm{d}\bm{x}.$ $$

If $\nabla\cdot\bm{z}=0,$ then

$$ $c(\bm{z},\bm{u},\bm{v})=-c(\bm{z},\bm{v},\bm{u}),$ (E.1) $$

and specifically $c(\bm{z},\bm{v},\bm{v})=0$.

###### Proof.

We note that this result is common in NSE literature, see e.g., [, Part I Sec 2.3], for homogeneous boundary or whole space. In this case, it suffices to verify that the boundary term vanishes. Here we still provide a short argument for the convenience of readers.
First by the product rule, and a vector calculus identity of $\nabla\cdot(\bm{v}\cdot\bm{u})$ (see e.g., [, Appendix II.3.2])

$$ $\displaystyle\nabla\cdot\left((\bm{v}\cdot\bm{u})\bm{z}\right)=$ $\displaystyle\left((\bm{z}\cdot\nabla)\bm{u}\right)\cdot\bm{v}+\left((\bm{z} \cdot\nabla)\bm{v}\right)\cdot\bm{u}+(\bm{v}\cdot\bm{u})\nabla\cdot\bm{z}.$ (E.2) $$

By Gauss divergence theorem, we then obtain:

$$ $\int_{\Omega}\nabla\cdot\left((\bm{v}\cdot\bm{u})\bm{z}\right)\,\mathrm{d}\bm{ x}=\int_{\partial\Omega}(\bm{v}\cdot\bm{u})\bm{z}\cdot\bm{n}\,\mathrm{d}s=\sum _{e_{i}\subset\partial\Omega,1\leq i\leq 4}\int_{\partial\Omega}(\bm{v}\cdot \bm{u})\bm{z}\cdot\bm{n}_{i}\,\mathrm{d}s,$ $$

with $\bm{n}$ is outer normal with respect to $\partial\Omega$, and $\bm{n}_{i}$ denotes that with respect to edge $e_{i}$. Now on $e_{1}:=\{x=0\}\times\{y\in(0,1)\}$, $\bm{n}_{1}=(-1,0)^{\top}$, $\bm{w}|_{e_{1}}=\bm{w}|_{e_{2}}$ by PBC where $e_{2}:=\{x=1\}\times\{y\in(0,1)\}$ and $\bm{w}\in\{\bm{z},\bm{u},\bm{v}\}$. The integrals on $e_{1}$ and $e_{2}$ cancel with one another since $\bm{n}_{2}=-\bm{n}_{1}$.
Furthermore, if $\nabla\cdot\bm{z}=0$, integrating (E.2) on $\Omega$ yields the desired result.
∎

###### Lemma E.3 (Poincaré inequality) .

$\mathcal{V}/\mathbb{R}\Subset\mathcal{H}/\mathbb{R}$ is compact and the following Poincaré inequality holds with constant $1$:
for any $v\in\mathcal{V}/\mathbb{R}$

$$ $\|v\|\leq|v|_{1}$ (E.3) $$

###### Proof.

A common textbook proof exploits the equivalence of sequential compactness and compactness, once the norm topology is introduced, we prove this in a very intuitive way under our special setting: by definition (2.14), if $v\in\mathcal{V}/\mathbb{R}$ then $\hat{v}(0,0)=0$

$$ $\|v\|_{0}^{2}=\sum_{\bm{k}\in\mathbb{Z}^{2}\backslash\{\bm{0}\}}|\hat{v}(\bm{k })|^{2}\leq\sum_{\bm{k}\in\mathbb{Z}^{2}\backslash\{\bm{0}\}}|\bm{k}|^{2}|\hat {v}(\bm{k})|^{2}=|v|_{1}^{2}.$ $$

∎

###### Assumption E.1 (Assumptions for Theorems 2.1 , 2.2 , E.9 ) .

###### Lemma E.2 (Skew-symmetry of the trilinear term) .

###### Proof.

###### Lemma E.3 (Poincaré inequality) .

###### Proof.

### E.1 Proof of Theorem 2.1

To prove Theorem 2.1, we need Lemmas E.4, E.5, and E.6, all of which can be found in classical textbook on analysis and approximations of NSE, such as Temam [], or Girault & Raviart [] for homogeneous Dirichlet boundary condition for the velocity field or slip boundary condition for the (pseudo-)stress tensor. We shall complement some proofs with simple adaptations to PBC whenever necessary for the convenience of readers.

Before presenting any lemmas, we need the following definitions, a weighted $\bm{H}^{1}$-norm on $\mathcal{V}=\bm{H}^{1}(\mathbb{T}^{2})$ for $\alpha>0$ on $\Omega$ is defined as

$$ $\|\bm{v}\|_{\alpha,1}:=\left(\|\bm{v}\|_{L^{2}}^{2}+\|\alpha\nabla\bm{v}\|^{2} _{L^{2}}\right)^{1/2},\text{ and }|\bm{v}|_{\alpha,1}:=\|\alpha\nabla\bm{v}\|_ {L^{2}}.$ $$

For $\bm{f}\in\mathcal{V}^{\prime}$, define a weighted dual norm as follows:

$$ $\|\bm{f}\|_{\alpha^{-1},-1}:=\sup_{\bm{v}\in\mathcal{V}}\frac{\langle\bm{f}, \bm{v}\rangle}{\|\bm{v}\|_{\alpha,1}}$ (E.4) $$

###### Lemma E.4 (Contuinity and embedding results for the convection term) .

For trilinear convection term

$$ $\bigl{(}(\bm{u}\cdot\nabla)\bm{v},\bm{w}\bigr{)}\leq\|\bm{u}\|_{L^{4}}\,|\bm{v }|_{1}\|\bm{v}\|_{L^{4}}.\quad\text{ and }\quad\bigl{(}(\bm{u}\cdot\nabla)\bm{ v},\bm{w}\bigr{)}\leq|\bm{u}|_{1}\,|\bm{v}|_{1}|\bm{w}|_{1}$ (E.5) $$

###### Proof.

The first result is obtained by the Hölder inequality, and the second holds with suite Sobolev embedding results [, Part I, \origS2.3, Lemma 2.1].
∎

###### Lemma E.5 (Energy stability of NSE in a bounded domain) .

$$ $\|\bm{u}\|_{L^{\infty}(\mathcal{T};\mathcal{H})}^{2}\leq\|\bm{u}_{0}\|_{ \mathcal{H}}^{2}+\|\bm{f}\|_{L^{2}(\mathcal{T};\mathcal{V}^{\prime})}^{2}.$ (E.6) $$

###### Proof.

(E.17) is a classical result, see e.g., [, Section 3.1].
∎

###### Lemma E.6 (Fréchet derivative of the convection term) .

Given $\bm{v},\bm{u}\in\mathcal{V}$, the linearization of the difference of the convection term reads

$$ $\bm{v}\cdot\nabla\bm{v}=\bm{u}\cdot\nabla\bm{u}+D(\bm{u})(\bm{v}-\bm{u})+\bm{r }\text{ where }D(\bm{u})\bm{\xi}:=\bigl{(}\bm{\xi}\cdot\nabla\bigr{)}\bm{u}+( \bm{u}\cdot\nabla)\bm{\xi},$ (E.7) $$

is the Fréchet derivative
and $\|\bm{r}\|\lesssim\|\bm{u}-\bm{v}\|\,\|\bm{u}-\bm{v}\|_{1}$.

###### Proof.

This result is normally used without proof in linearizing NSE, see e.g., [, Chapter 4 eq. (6.5.2)]. Define $F(\bm{u},\bm{G}):=\bm{u}^{\top}\bm{G}$, then expanding $F$ at $\bm{v}$ and $\bm{G}=\nabla\bm{u}$ yields the desired result.
∎

###### Theorem 2.1 Part I (Functional-type a posterior error estimate is efficient (rigorous version)) .

Consider the Gelfand triple $\mathcal{V}\Subset\mathcal{H}\Subset\mathcal{V}^{\prime}$, and the weak solution $\bm{u}\in L^{2}(\mathcal{T};\mathcal{V})\cap L^{\infty}(\mathcal{T};\mathcal{H})$ to (2.2) be sufficiently regular, then the dual norm of the residual is efficient to estimate the error:

$$ $\|\mathcal{R}(\bm{u}_{\mathcal{N}})\|_{L^{2}(\mathcal{T};\mathcal{V}^{\prime}) }^{2}\lesssim\|\partial_{t}(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}_{L^{2}(\mathcal {T};\mathcal{V}^{\prime})}+\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}_{L^{2}(\mathcal {T};\mathcal{V})}.$ (E.8) $$

The constant depend on the regularity of the true solution $\bm{u}$.

###### Proof.

The proof of the lower bound (efficiency) (E.8) follows the skeleton laid out in [, Lemma 4.1] for diffusion and temporal derivative terms, while not needing to extend approximation using an affine linear extension operator as [] does. The treatment of the nonlinear convection term follows [, Lemma 7] for the steady-state NSE (no temporal derivatives), barring the technicality of the divergence-free condition.
Consider at $t\in\mathcal{T}$, for any test function $\bm{v}\in\bm{H}^{1}(\mathbb{T}^{2})$ and $\nabla\cdot\bm{v}=0$,
integrating $\bm{v}$ against $R(\bm{u}_{\mathcal{N}})$ in (2.9) yields

$$ $\displaystyle\langle R(\bm{u}_{\mathcal{N}}),\bm{v}\rangle=$ $\displaystyle\bigl{(}\bm{f}-\partial_{t}\bm{u}_{\mathcal{N}}-(\bm{u}_{\mathcal {N}}\cdot\nabla){\bm{u}_{\mathcal{N}}}+\nu\Delta{\bm{u}_{\mathcal{N}}},\bm{v} \bigr{)}$ $\displaystyle=$ $\displaystyle\bigl{(}\partial_{t}(\bm{u}-\bm{u}_{\mathcal{N}}),\bm{v}\bigr{)}- \bigl{(}\nu\Delta(\bm{u}-\bm{u}_{\mathcal{N}}),\bm{v}\bigr{)}$ $\displaystyle+\bigl{(}(\bm{u}\cdot\nabla){\bm{u}}-(\bm{u}_{\mathcal{N}}\cdot \nabla){\bm{u}_{\mathcal{N}}},\bm{v}\bigr{)}.$ $$

Integrating by parts for the diffusion term, using the same argument for PBC as in Lemma E.2, and inserting $(\bm{u}\cdot\nabla){\bm{u}_{\mathcal{N}}}$ yield

$$ $\displaystyle\langle R(\bm{u}_{\mathcal{N}}),\bm{v}\rangle=$ $\displaystyle\underbrace{\bigl{(}\partial_{t}(\bm{u}-\bm{u}_{\mathcal{N}}),\bm {v}\bigr{)}}_{(\text{I})}+\underbrace{\bigl{(}\nu\nabla(\bm{u}-\bm{u}_{ \mathcal{N}}),\nabla\bm{v}\bigr{)}}_{(\text{II})}$ (E.9) $\displaystyle+\underbrace{\bigl{(}(\bm{u}\cdot\nabla)(\bm{u}-\bm{u}_{\mathcal{ N}}),\bm{v}\bigr{)}}_{(\text{III})}+\underbrace{\left(\left((\bm{u}-\bm{u}_{ \mathcal{N}})\cdot\nabla\right){\bm{u}_{\mathcal{N}}},\bm{v}\right)}_{(\text{ IV})}.$ $$

Applying the definition of a weighted dual norm (E.4) on $(\text{I})$ we have

$$ $(\text{I})\leq\left\|\partial_{t}(\bm{u}-\bm{u}_{\mathcal{N}})(t,\cdot)\right \|_{\nu^{-1/2},-1}\|\bm{v}\|_{\nu^{1/2},1}.$ $$

For $(\text{II})$ we simply have

$$ $(\text{II})\leq\|\nu^{1/2}\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|_{L^{2}}\|\nu^{ 1/2}\nabla\bm{v}\|_{L^{2}}.$ $$

For $(\text{III})$, applying Lemma E.4 and the stability estimate in Lemma E.5 for $\bm{u}$:

$$ $(\text{III})\leq C\|\bm{u}\|_{\nu^{-1},1}\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{\nu^ {1/2},1}\|\bm{v}\|_{\nu^{1/2},1}\leq\nu^{-1}C_{1}(t,\bm{u}_{0})\|\bm{u}-\bm{u} _{\mathcal{N}}\|_{\nu^{1/2},1}\|\bm{v}\|_{\nu^{1/2},1}.$ $$

For $(\text{IV})$, applying Lemma E.4 and the stability estimate in Lemma E.5 for $\bm{u}_{\mathcal{N}}$:

$$ $(\text{IV})\leq C\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{\nu^{1/2},1}\|\bm{u}_{ \mathcal{N}}\|_{\nu^{-1},1}\|\bm{v}\|_{\nu^{1/2},1}\leq\nu^{-1}C_{2}(t,\bm{u}_ {0})\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{\nu^{1/2},1}\|\bm{v}\|_{\nu^{1/2},1}.$ $$

Applying the Cauchy-Schwarz inequality, one has

$$ $|\langle R(\bm{u}_{\mathcal{N}}),\bm{v}\rangle|^{2}\leq C\left(\left\|\partial _{t}(\bm{u}-\bm{u}_{\mathcal{N}})(t,\cdot)\right\|_{\nu^{-1/2},-1}^{2}+\|\bm{u }-\bm{u}_{\mathcal{N}}\|_{\nu^{1/2},1}^{2}\right)\|\bm{v}\|_{\nu^{1/2},1}^{2}.$ $$

Finally, equipping $\mathcal{V}^{\prime}$ with the weighted norm (E.4), and by the definition of $L^{2}(\mathcal{T};\mathcal{V}^{\prime})$, we have

$$ $\displaystyle\|R(\bm{u}_{\mathcal{N}})\|_{L^{2}(\mathcal{T};\mathcal{V})}^{2}=$ $\displaystyle\;\int_{\mathcal{T}}\sup_{\bm{v}\in\mathcal{V}}\frac{|\langle R( \bm{u}_{\mathcal{N}}),\bm{v}\rangle|^{2}}{\|\bm{v}\|_{\nu^{1/2},1}^{2}}\, \mathrm{d}t$ $\displaystyle\leq$ $\displaystyle\;C\int_{\mathcal{T}}\left(\left\|\partial_{t}(\bm{u}-\bm{u}_{ \mathcal{N}})(t,\cdot)\right\|_{\nu^{-1/2},-1}^{2}+\|\bm{u}-\bm{u}_{\mathcal{N }}\|_{\nu^{1/2},1}^{2}\right)\mathrm{d}t,$ $$

where $C=C(\bm{u}_{0},\nu,\mathcal{T})$.
∎

###### Theorem 2.1 Part II (Functional-type a posterior error estimate is reliable (rigorous version)) .

Consider the Gelfand triple $\mathcal{V}\Subset\mathcal{H}\Subset\mathcal{V}^{\prime}$, and the weak solution is $\bm{u}\in L^{2}(\mathcal{T};\mathcal{V})\cap L^{\infty}(\mathcal{T};\mathcal{H})$ to (2.2).
We assume that

- $($2.1.1$)$
$\bm{u}_{\mathcal{N}}$ and $\bm{u}$ are sufficiently close in the sense that: there exists $\gamma\in[0,C)$ for a fixed $C\in\mathbb{R}^{+}$, for $\bm{J}$ defined in Lemma E.6
$\|\bm{J}(\bm{u},\bm{u}_{\mathcal{N}})\|_{\mathcal{V}^{\prime}}\leq\gamma\|
\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|.$
(E.10)
- $($2.1.2$)$
$\bm{u}$ satisfies the Gårding inequality (a weaker coercivity): for any $\bm{v}(t,\cdot)\in\mathcal{V}$, define
$B(\bm{v},\bm{w};\bm{u}):=(\partial_{t}\bm{v},\bm{w})+\nu(\nabla\bm{v},\nabla
\bm{w})+\bigl{(}(\bm{u}\cdot\nabla)\bm{v},\bm{w}\bigr{)},$
and for $\bm{v}$ in a neighborhood such that Assumption $($2.1.1$)$ holds, there exists $\alpha,\beta>0$ such that $\alpha-\beta\geq\nu$,
$B(\bm{v},\bm{v};\bm{u})+\beta\|\bm{v}\|^{2}\geq\frac{\mathrm{d}}{\mathrm{d}t}
\|\bm{v}\|^{2}+\alpha\|\nabla\bm{v}\|^{2}.$
(E.11)

The dual norm of the residual is then reliable to serve as an error measure in the following sense: denote $\mathcal{T}_{m}:=(t_{m},t_{m+1}]$

$$ $\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}_{L^{\infty}(\mathcal{T}_{m};\mathcal{H})}+ \|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}_{L^{2}(\mathcal{T}_{m};\mathcal{V})}\leq \bigl{\|}(\bm{u}-\bm{u}_{\mathcal{N}})(t_{m},\cdot)\bigr{\|}^{2}_{\mathcal{V}} +C\int_{\mathcal{T}_{m}}\|R(\bm{u}_{\mathcal{N}})(t,\cdot)\|_{\mathcal{V}^{ \prime}}^{2}\,\mathrm{d}t.$ (E.12) $$

The constant $C=C(\nu)$.

###### Proof.

The proof of part II of Theorem 2.1 combines the meta-framework of [] for time-dependent problems and [, Section 8] for the stationary NSE. Note that thanks for the construction of divergence-free $\bm{u}_{\mathcal{N}}$, the technicality of applying the Ladyzhenskaya-Babušška-Brezzi inf-sup condition is avoided. In the meantime, no interpolations are needed to convert the functional norm for the negative Sobolev spaces to an $L^{2}$ representation.

To prove the upper bound (reliability) (2.12), we simply choose a time step $t=t_{m}$, and let $\bm{v}=(\bm{u}-\bm{u}_{\mathcal{N}})(t,\cdot)$ on $(t_{m},t_{m+1})$. Note that the discretized approximation’s evaluation at $t$, which may not be the temporal grids $t_{m}$ or $t_{m+1}$ are naturally defined using the basis in (2.8). Using the error representation in (E.9), we have the $(\text{I})$ and $(\text{II})$ terms in $\left\langle R(\bm{u}_{\mathcal{N}}),\bm{u}-\bm{u}_{\mathcal{N}}\right\rangle$ are

$$ $\displaystyle(\text{I})+(\text{II})=$ $\displaystyle\;\bigl{(}\partial_{t}(\bm{u}-\bm{u}_{\mathcal{N}}),\bm{u}-\bm{u} _{\mathcal{N}}\bigr{)}+\bigl{(}\nu\nabla(\bm{u}-\bm{u}_{\mathcal{N}}),\nabla( \bm{u}-\bm{u}_{\mathcal{N}})\bigr{)}$ (E.13) $\displaystyle=$ $\displaystyle\;\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\|\bm{u}-\bm{u}_{ \mathcal{N}}\|_{L^{2}}^{2}+\nu\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|_{L^{2}}^ {2}$ $$

For the $(\text{III})$ and $(\text{IV})$ terms we have

$$ $(\text{III})+(\text{IV})=\overbrace{\bigl{(}(\bm{u}\cdot\nabla)(\bm{u}-\bm{u}_ {\mathcal{N}}),\bm{u}-\bm{u}_{\mathcal{N}}\bigr{)}}^{=0\text{ by Lemma \ref{ lemma:skew-symmetry}}}+\overbrace{\left(\left((\bm{u}-\bm{u}_{\mathcal{N}}) \cdot\nabla\right){\bm{u}_{\mathcal{N}}},\bm{u}-\bm{u}_{\mathcal{N}}\right)}^{ =:(*)}.$ (E.14) $$

Note by Lemma E.2, $(*)$ is also:

$$ $(*)=\left(\left((\bm{u}-\bm{u}_{\mathcal{N}})\cdot\nabla\right){\bm{u}},\bm{u} -\bm{u}_{\mathcal{N}}\right),$ (E.15) $$

since $c(\bm{z},\bm{z},\bm{z})=0$ for $\bm{z}=\bm{u}-\bm{u}_{\mathcal{N}}$.
Consequently, by a vector calculus identity, we have

$$ $\int_{\Omega}(\bm{z}\cdot\nabla)\bm{u}\cdot\bm{z}\,\mathrm{d}\bm{x}=\int_{ \Omega}\nabla\bm{u}:\left(\bm{z}\!\otimes\!\bm{z}\right)\mathrm{d}\bm{x}.$ $$

Thus, we reach the following error equation:

$$ $\left\langle R(\bm{u}_{\mathcal{N}}),\bm{u}-\bm{u}_{\mathcal{N}}\right\rangle= \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{L^{2 }}^{2}+\nu\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|_{L^{2}}^{2}+\int_{\Omega} \nabla\bm{u}:\left((\bm{u}-\bm{u}_{\mathcal{N}})\!\otimes\!(\bm{u}-\bm{u}_{ \mathcal{N}})\right)\mathrm{d}\bm{x}.$ (E.16) $$

Now, by Assumption $($2.1.2$)$, Poincaré inequality from Lemma E.3,

$$ $\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\|\bm{u}-\bm{u}_{ \mathcal{N}}\|^{2}+\nu\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}\leq$ $\displaystyle\;\frac{\mathrm{d}}{\mathrm{d}t}\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{ 2}+(\alpha-\beta)\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}$ $\displaystyle\leq$ $\displaystyle\;\frac{\mathrm{d}}{\mathrm{d}t}\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{ 2}+\alpha\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}-\beta\|\bm{u}-\bm{u}_{ \mathcal{N}}\|$ $\displaystyle\leq$ $\displaystyle\;B(\bm{u}-\bm{u}_{\mathcal{N}},\bm{u}-\bm{u}_{\mathcal{N}};\bm{u }).$ $$

By (E.13), (E.14), and (E.15),

$$ $\displaystyle B(\bm{u}-\bm{u}_{\mathcal{N}},\bm{u}-\bm{u}_{\mathcal{N}};\bm{u})=$ $\displaystyle\;(\text{I})+(\text{II})+(\text{III})+(\text{IV})=\langle R(\bm{u }_{\mathcal{N}}),\bm{u}-\bm{u}_{\mathcal{N}}\rangle$ $\displaystyle\leq$ $\displaystyle\;\|R(\bm{u}_{\mathcal{N}})\|_{\nu^{-1/2},-1}|\bm{u}-\bm{u}_{ \mathcal{N}}|_{\nu^{1/2},1}$ $\displaystyle\leq$ $\displaystyle\;\frac{1}{2}\|R(\bm{u}_{\mathcal{N}})\|_{\nu^{-1/2},-1}^{2}+ \frac{1}{2}|\bm{u}-\bm{u}_{\mathcal{N}}|_{\nu^{1/2},1}^{2}.$ $$

As a result, we have

$$ $\frac{\mathrm{d}}{\mathrm{d}t}\|\bm{u}-\bm{u}_{\mathcal{N}}\|^{2}+\nu\|\nabla( \bm{u}-\bm{u}_{\mathcal{N}})\|^{2}\leq\|R(\bm{u}_{\mathcal{N}})\|_{\nu^{-1/2}, -1}^{2}.$ (E.17) $$

Integrating from $t_{m}$ to any $t\in(t_{m},t_{m+1}]$ yields

$$ $\|(\bm{u}-\bm{u}_{\mathcal{N}})(t,\cdot)\|^{2}-\|(\bm{u}-\bm{u}_{\mathcal{N}}) (t_{m},\cdot)\|^{2}+\nu\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}_{L^{2}(t_{m },t;\mathcal{H})}\leq\int_{t_{m}}^{t}\|R(\bm{u}_{\mathcal{N}})(t,\cdot)\|_{\nu ^{-1/2},-1}^{2}\,\mathrm{d}t.$ $$

Since $t\in(t_{m},t_{m+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 to prove the reliability is 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,

$$ $\text{norm of the error}\simeq(\text{I})+(\text{II})=\langle R(\bm{u}_{ \mathcal{N}}),\bm{u}-\bm{u}_{\mathcal{N}}\rangle-(*)\not\leq\langle R(\bm{u}_{ \mathcal{N}}),\bm{u}-\bm{u}_{\mathcal{N}}\rangle$ $$

We note that Assumption $($2.1.2$)$ is equivalent to putting a threshold on $\|\bm{u}(t,\cdot)\|_{L^{\infty}}$, which is commonly assumed in the analysis of numerical methods for convection-dominated problems.
A stronger alternative assumption than Assumption $($2.1.2$)$ would be imposing a stronger constraint on $\gamma$ in Assumption $($2.1.1$)$
Using the “sufficient close” assumption (E.10) on the nonlinear term, we have

$$ $|(*)|\leq\|\left((\bm{u}-\bm{u}_{\mathcal{N}})\cdot\nabla\right){\bm{u}_{ \mathcal{N}}}\|_{\mathcal{V}^{\prime}}\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{ \mathcal{V}}\leq\gamma\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2}.$ $$

Now if $\bar{\nu}:=\nu-\gamma\geq 0$, one would reach a similar estimate as (E.17) by replacing $\nu$ with $\bar{\nu}$.
If no extra assumption on $\gamma$ is imposed.

Another way is to consider a Sobolev embedding directly after applying the Hölder inequality, for (E.15) we have

$$ $|(*)|\leq\int_{\Omega}|\nabla\bm{u}||\bm{u}-\bm{u}_{\mathcal{N}}|^{2}\leq\| \nabla\bm{u}\|\,\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{L^{4}}^{2}\leq\rho\|\nabla\bm {u}\|\,\|\nabla(\bm{u}-\bm{u}_{\mathcal{N}})\|^{2},$ $$

where $\rho$ is the constant in the following Sobolev embedding

$$ $\|\bm{u}-\bm{u}_{\mathcal{N}}\|_{L^{4}}\leq\rho\|\nabla(\bm{u}-\bm{u}_{ \mathcal{N}})\|.$ $$

Now if $\tilde{\nu}:=\nu-\rho\|\nabla\bm{u}\|\geq 0$, one would reach a similar estimate as (E.17) by replacing $\nu$ with $\tilde{\nu}$.

###### Lemma E.4 (Contuinity and embedding results for the convection term) .

###### Proof.

###### Lemma E.5 (Energy stability of NSE in a bounded domain) .

###### Proof.

###### Lemma E.6 (Fréchet derivative of the convection term) .

###### Proof.

###### Theorem 2.1 Part I (Functional-type a posterior error estimate is efficient (rigorous version)) .

###### Proof.

###### Theorem 2.1 Part II (Functional-type a posterior error estimate is reliable (rigorous version)) .

###### Proof.

###### Remark E.7 (Necessity of the “sufficient close” and the regularity conditions) .

### E.2 Proof of Theorem 2.2

###### Theorem 2.2 (Equivalence of functional norm and negative Sobolev norm) .

Consider the Gelfand triple $\mathcal{V}\Subset\mathcal{H}\Subset\mathcal{V}^{\prime}$, then

$$ $\|f\|_{\mathcal{V}^{\prime}}=|f|_{-1}\;\text{ for }f\in\mathcal{H}/\mathbb{R}.$ (E.18) $$

###### Proof.

This proof leverages the spectral basis for $\mathcal{H}$ without the extra technicality of Schwartz space involved if one ought to assume that $f\in\mathcal{V}^{\prime}$. Recall (2.8), and define an infinite version

$$ $\operatorname{\mathcal{S}}_{\infty}:=\operatorname{span}\left\{e^{i\bm{k}\cdot \bm{x}}:\bm{k}\in 2\pi\mathbb{Z}^{2}\right\}.$ (E.19) $$

In what follows, we 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\in\mathcal{H}/\mathbb{R}$ with its Fourier series in $\operatorname{\mathcal{S}}_{\infty}$,

$$ $f(\bm{x})=\sum_{\bm{k}\in 2\pi\mathbb{Z}^{2}}\hat{f}(\bm{k})e^{i\bm{k}\cdot\bm {x}}.$ $$

Note since $f$ lies in the quotient space, $\hat{f}(0,0)=0$. Then, the duality pairing can be identified as an $\bm{L}^{2}$-inner product, as well as the test function $v\in\mathcal{H}$,

$$ $\displaystyle\langle f,v\rangle_{\mathcal{V}^{\prime},\mathcal{V}}=$ $\displaystyle\;(f,v)=\left(\sum_{\bm{k}\in 2\pi\mathbb{Z}^{2}}\hat{f}(\bm{k})e ^{i\bm{k}\cdot\bm{x}},\sum_{\bm{m}\in 2\pi\mathbb{Z}^{2}}\overline{\hat{v}(\bm {m})e^{i\bm{m}\cdot\bm{x}}}\right)$ $\displaystyle=$ $\displaystyle\;\sum_{\bm{k}\in 2\pi\mathbb{Z}_{n}^{2}\backslash\{0\}}\hat{f}( \bm{k})\overline{\hat{v}(\bm{k})}$ $\displaystyle\leq$ $\displaystyle\left(\sum_{\bm{k}\in 2\pi\mathbb{Z}_{n}^{2}\backslash\{0\}}|\bm{ k}|^{-2}|\hat{f}(\bm{k})|^{2}\right)^{1/2}\left(\sum_{\bm{k}\in 2\pi\mathbb{Z} _{n}^{2}\backslash\{0\}}|\bm{k}|^{2}|\hat{v}(\bm{k})|^{2}\right)^{1/2}$ $\displaystyle\leq$ $\displaystyle\;|f|_{-1}|v|_{1}.$ $$

As a result, by the definition in (2.14) and the estimate above

$$ $\|f\|_{\mathcal{V}^{\prime}}=\sup_{v\in\mathcal{V}/\mathbb{R},|v|_{\mathcal{V} }=1}|\langle f,v\rangle|=\sup_{v\in\mathcal{V}/\mathbb{R}}\frac{(f,v)}{|v|_{1} }\leq|f|_{-1}.$ $$

To show the other direction, let $v=(-\Delta)^{-1}f\in H^{1}(\mathbb{T}^{2})$ by the well-posedness of $-\Delta v=f$ for $f\in L^{2}(\mathbb{T}^{2})$. We have

$$ $\displaystyle\langle f,(-\Delta)^{-1}f\rangle$ $\displaystyle=\bigl{(}-\Delta v,v\bigr{)}=\bigl{(}f,(-\Delta)^{-1}f\bigr{)}$ $\displaystyle=\left(\sum_{\bm{k}\in 2\pi\mathbb{Z}^{2}}\hat{f}(\bm{k})e^{i\bm{ k}\cdot\bm{x}},\sum_{\bm{m}\in 2\pi\mathbb{Z}^{2}}(-\Delta_{\bm{x}})^{-1} \overline{\hat{f}(\bm{m})e^{i\bm{m}\cdot\bm{x}}}\right)$ $\displaystyle=\left(\sum_{\bm{k}\in 2\pi\mathbb{Z}^{2}}\hat{f}(\bm{k})e^{i\bm{ k}\cdot\bm{x}},\sum_{\bm{m}\in 2\pi\mathbb{Z}^{2}}\frac{1}{|\bm{m}|^{2}} \overline{\hat{f}(\bm{m})e^{i\bm{m}\cdot\bm{x}}}\right)$ $\displaystyle=\sum_{\bm{k}\in 2\pi\mathbb{Z}_{n}^{2}\backslash\{0\}}|\bm{k}|^{ -2}|\hat{f}(\bm{k})|^{2}=|f|_{-1}^{2}.$ $$

Note that by the construction of $v$ above, $|v|_{1}=|f|_{-1}$, thus,

$$ $\sup_{v\in\mathcal{V}}\frac{(f,v)}{|v|_{1}}\geq\frac{\bigl{(}f,(-\Delta)^{-1}f \bigr{)}}{|v|_{1}}=|f|_{-1},$ $$

which proves the theorem.

∎

###### Theorem 2.2 (Equivalence of functional norm and negative Sobolev norm) .

###### Proof.

### E.3 Convergence result of fine-tuning

In the following theorem, Theorem E.9, we show that a sufficient condition for the optimizer to converge is to get in a neighborhood of the true solution, thus corroborating the necessity of training.
Note in both Theorems 2.1 and E.9, the error term of $\bm{u}^{(m)}_{\mathcal{N}}-\bm{u}(t_{m},\cdot)$ is present. We have to acknowledge that the “initial value” for the predicted trajectory, which is the last snapshot in the input trajectory, may have errors. This is the major motivation that we opt to use $\bm{u}^{(\ell)}_{\mathcal{N}}$, which is the input trajectory’s last snapshot in the skip-connection in $\widetilde{Q}$

$$ $\bm{u}^{(m+1)}_{\mathcal{N}}\leftarrow\bm{u}^{(\ell)}_{\mathcal{N}}+\widetilde {Q}_{\theta}(\bm{v}_{\mathcal{N}}^{(m+1)}),\;\text{ for }m=\ell,\dots,\ell+n_{ t}-1.$ $$

In some sense, SFNO learns the *derivative* $\partial_{t}\bm{u}$’s arbitrary-lengthed integral. If 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 $\widetilde{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 (2.8) under the functional norm. Let
$\theta^{*}=\arg\min_{\theta}\|\mathcal{R}(\bm{u}_{\mathcal{N}}(\theta))(t_{m+1
},\cdot)\|_{-1}$, then $\widetilde{Q}_{\theta^{*}}(\bm{u}_{\mathcal{N}})=\arg\min_{\bm{v}\in
\operatorname{\mathcal{S}}}\|\bm{u}(t_{m+1},\cdot)-(\bm{u}(t_{m},\cdot)+\bm{v}
)\|_{*}$.

Figure: Algorithm 2 An error-guarantee fine-tuning strategy.

###### Lemma E.8 (Local strict convexity for the fine-tuning loss) .

Define $\|\cdot\|_{*,\delta}$ to be a (dual) graph norm on $L^{\infty}\big{(}\mathcal{T}_{\delta};\bm{L}^{2}(\mathbb{T}^{2})\big{)}\cap L^
{2}\big{(}\mathcal{T}_{\delta};\bm{H}^{1}(\mathbb{T}^{2})\big{)}$ , where $\mathcal{T}_{\delta}:=[t-\delta,t+\delta]$

$$ $\displaystyle\|\bm{v}\|_{*,\delta}=\left\{\fint_{\mathcal{T}_{\delta}}\|\bm{v} (t,\cdot)\|^{2}\mathrm{d}t+\fint_{\mathcal{T}_{\delta}}\left\|\left(\partial_{ t}\bm{v}+\bm{v}\cdot\nabla\bm{v}\right)(t,\cdot)\right\|_{\mathcal{V}^{\prime} }^{2}\mathrm{d}t\right\}^{1/2}$ $$

For $\bm{u}\in L^{\infty}\big{(}\mathcal{T}_{\delta};\bm{L}^{2}(\mathbb{T}^{2})\big
{)}\cap L^{2}\big{(}\mathcal{T}_{\delta};\bm{H}^{1}(\mathbb{T}^{2})\big{)}$ the weak solution to (2.2) on $\mathcal{T}_{\delta}$ that is sufficiently smooth, there exists $\delta,\epsilon\in\mathbb{R}^{+}$, such that on

$$ $\mathcal{B}(\bm{u};\epsilon):=\left\{\bm{v}\in L^{\infty}\big{(}\mathcal{T}_{ \delta};\bm{L}^{2}(\mathbb{T}^{2})\big{)}\cap L^{2}\big{(}\mathcal{T}_{\delta} ;\bm{H}^{1}(\mathbb{T}^{2})\big{)}:\|\bm{u}-\bm{v}\|_{*,\delta}\leq\epsilon \right\},$ $$

the functional

$$ $J(\bm{v}(t,\cdot)):=\frac{1}{2}\|R(\bm{v})\|_{\mathcal{V}^{\prime}}^{2}\;\text { where }\;R(\bm{v}):=\bm{f}-\partial_{t}\bm{v}-(\bm{v}\cdot\nabla){\bm{v}}+ \nu\Delta{\bm{v}}$ $$

is strictly convex.

###### Proof.

First, by Theorem 2.2,

$$ $J(\bm{v})=\frac{1}{2}\left\langle R(\bm{v}),(-\Delta)^{-1}R(\bm{v})\right \rangle_{\mathcal{V}^{\prime},\mathcal{V}}$ $$

Then,

$$ $\displaystyle DJ(\bm{v};\bm{\xi}):=$ $\displaystyle\;\lim_{\tau\to 0}{\frac{J(\bm{v}+\tau\bm{\xi})-J(\bm{v})}{\tau}} =\left.{\frac{\mathrm{d}}{\mathrm{d}\tau}}J(\bm{v}+\tau\bm{\xi})\right|_{\tau=0}$ $\displaystyle=$ $\displaystyle\;\langle DR(\bm{v})\bm{\xi},(-\Delta)^{-1}R(\bm{v})\rangle_{ \mathcal{V}^{\prime},\mathcal{V}}.$ $$

where

$$ $DR(\bm{v})\bm{\xi}=\partial_{t}\bm{\xi}+\frac{1}{2}\bigl{(}\bigl{(}\bm{\xi} \cdot\nabla\bigr{)}\bm{v}+(\bm{v}\cdot\nabla)\bm{\xi}\bigr{)}-\nu\Delta\bm{\xi}$ $$

is the Fréchet derivative $DR(\bm{v}):\mathcal{V}\to\mathcal{V}^{\prime}$ by Lemma E.6. The Hessian is then

$$ $\operatorname{Hess}J(\bm{v};\bm{\xi},\bm{\zeta})=\left\langle DR(\bm{v})\bm{ \xi},(-\Delta)^{-1}DR(\bm{v})\bm{\zeta}\right\rangle+\left\langle\bm{\zeta} \cdot D^{2}R(\bm{v})\bm{\xi},(-\Delta)^{-1}R(\bm{v})\right\rangle.$ $$

Now, since $R(\bm{u})=0$ on $\mathcal{T}_{\delta}$, we have

$$ $\operatorname{Hess}J(\bm{u};\bm{\xi},\bm{\xi})=\left\langle DR(\bm{v})\bm{\xi} ,(-\Delta)^{-1}DR(\bm{v})\bm{\xi}\right\rangle.$ $$

If we assume that $DR(\bm{u})\bm{\xi}\in\mathcal{V}^{\prime}$ has its functional norm bounded above and below in $\mathcal{B}(\bm{u};\epsilon)$, one has for any $\bm{\xi}$

$$ $\|\bm{\xi}\|_{\mathcal{V}}^{2}\lesssim\|DR(\bm{u})\bm{\xi}\|_{\mathcal{V}^{ \prime}}^{2}\leq\operatorname{Hess}J(\bm{u};\bm{\xi},\bm{\xi}).$ $$

Simply choosing $\epsilon$ small enough such that $\bm{v}$ is sufficiently close to $\bm{u}$ in graph norm associated with the PDE to make the coercivity above still true for $\operatorname{Hess}J(\bm{v};\bm{\xi},\bm{\xi})$ yields the desired local convexity.
∎

###### Theorem E.9 (Guaranteed convergence of the fine-tuning) .

In addition to the same assumptions with Theorem 2.1, suppose Lemma E.8 holds for $\epsilon\in(0,1)$,
and a given $\bm{u}_{\mathcal{N}}$ can be embedded in $\mathcal{B}(\bm{u};\epsilon^{\prime})$ for a $0<\epsilon^{\prime}\leq\epsilon$. Denote $\bm{u}^{(k)}_{\mathcal{N},j}$ the evaluation in Line 11 of Algorithm 2, and $j$ the iteration of optimizer in Line 12, then the fine-tuning using the new loss function (2.16) produces a sequence
$\{\bm{u}^{(k)}_{\mathcal{N},j}\}_{j=1}^{\infty}\subset\mathcal{B}\big{(}\bm{u}
;\epsilon^{\prime}\big{)}$
Furthermore, suppose that the optimizer in Line 12 of Algorithm 2 has a learning rate converging to 0 0. Then, then the fine-tuning using the new loss function (2.16)
produces a sequence of evaluations
converging to the best possible approximation $\bm{u}^{(m+1)}_{\mathcal{N},\infty}\in\operatorname{\mathcal{S}}$ of $\bm{u}(t_{m+1},\cdot)$ starting from $\bm{u}^{(m)}_{\mathcal{N}}$, in the sense that for $m=\ell,\dots,\ell+n_{t}-1$

$$ $\|\bm{u}^{(m+1)}_{\mathcal{N},\infty}-\bm{u}(t_{m+1},\cdot)\|_{\mathcal{V}} \leq\|\bm{u}^{(l)}_{\mathcal{N}}-\bm{u}(t_{m},\cdot)\|_{\mathcal{V}}+c_{1}n_{t }n^{-2}|\bm{u}(t_{m+1},\cdot)|_{2}+c_{2}n_{t}(\Delta t)^{\alpha-1}.$ (E.20) $$

###### Proof.

For simplicity, we denote $\bm{u}_{0}:=\bm{u}(t_{l},\cdot)\in\mathcal{V}$, and $\bm{u}_{\mathcal{N}}:=\bm{u}^{(m)}_{\mathcal{N}}\in\operatorname{\mathcal{S}}|
_{t=t_{m}}=:\operatorname{\mathcal{S}}$. First, by line 12 fine-tuning algorithm, if one solves the optimization in the functional norm exactly, we have

$$ $\theta^{*}=\arg\min_{\theta}\delta t\left\|\overline{R}\bigl{(}\mathcal{G}_{ \alpha}\bm{u}_{\mathcal{N}}(\theta)\bigr{)}\right\|_{\operatorname{\mathcal{S} }^{\prime}}^{2}\;\text{ where }\;\bm{u}_{\mathcal{N}}(\theta):=\bm{u}_{ \mathcal{N}}+\tilde{Q}_{\theta}(\bm{v}_{\mathcal{N}})$ (E.21) $$

where $\overline{R}(\cdot)$ is the residual functional computed using the time derivative term from an extra-fine-step solver’s result

$$ $\overline{R}(\bm{v}):=\bm{f}-D_{t}\bm{v}-(\bm{v}\cdot\nabla){\bm{v}}+\nu\Delta {\bm{v}}$ $$

Unlike representing the derivative using output from the neural operator, one of the keys of our algorithm is that the error for $\partial_{t}\bm{u}-D_{t}\bm{u}_{\mathcal{N}}(\theta)$ can be explicitly estimated using the framework to develop estimates for truncation error 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 $\bm{v}_{\mathcal{N}}$, $\tilde{Q}_{\theta}(\bm{v}_{\mathcal{N}})\in\operatorname{\mathcal{S}}$, as a result, solving (E.21) is equivalent to solve the following: denote $\bm{u}_{\delta}:=\mathcal{G}_{\alpha}(\bm{u}_{\mathcal{N}}+\delta\bm{u})$ for $\delta\bm{u}\in\operatorname{\mathcal{S}}$

$$ $\min_{\delta\bm{u}\in\operatorname{\mathcal{S}}}(\delta t)^{1/2}\left\|D_{t}( \bm{u}_{\delta})+(\bm{u}_{\delta}\cdot\nabla){\bm{u}_{\delta}}-\nu\Delta{\bm{u }_{\delta}}-\bm{f}\right\|_{\operatorname{\mathcal{S}}^{\prime}}.$ $$

Replacing $D_{t}(\cdot)$ by $\partial_{t}(\cdot)$ we have a truncation error term, whose error is of order
$(\Delta t)^{\alpha-1}$ due to taking the time derivative:

$$ $\displaystyle\left\|D_{t}(\bm{u}_{\delta})+(\bm{u}_{\delta}\cdot\nabla){\bm{u} _{\delta}}-\nu\Delta{\bm{u}_{\delta}}-\bm{f}\right\|_{\operatorname{\mathcal{S }}^{\prime}}$ $\displaystyle\leq$ $\displaystyle\;\left\|\partial_{t}\bm{u}_{\delta}+(\bm{u}_{\delta}\cdot\nabla) {\bm{u}_{\delta}}-\nu\Delta{\bm{u}_{\delta}}-\bm{f}\right\|_{\operatorname{ \mathcal{S}}^{\prime}}+\left\|D_{t}(\bm{u}_{\delta})-\partial_{t}\bm{u}_{ \delta}\right\|_{\operatorname{\mathcal{S}}^{\prime}}$ $$

We focus on estimating the first term above,

$$ $\min_{\delta\bm{u}\in\operatorname{\mathcal{S}}}(\delta t)^{1/2}\left\| \partial_{t}\bm{u}_{\delta}+(\bm{u}_{\delta}\cdot\nabla){\bm{u}_{\delta}}-\nu \Delta{\bm{u}_{\delta}}-\bm{f}\right\|_{\operatorname{\mathcal{S}}^{\prime}}$ (E.22) $$

By the fact that $\|\cdot\|_{\operatorname{\mathcal{S}}^{\prime}}$ inherit the scaling law and the triangle inequality from $\|\cdot\|_{\mathcal{V}^{\prime}}$, it is convex as well in this neighborhood. As a result,
any gradient-based optimizer with a converging step size shall converge to the minimum, achieved at $\bm{u}^{(m)}_{\mathcal{N},\infty}$, with a linear convergence rate. Moreover, $\bm{u}^{(m)}_{\mathcal{N},\infty}:=\bm{u}_{\mathcal{N}}^{(m)}+(\delta\bm{u})^{*}$ is the (nonlinear) Galerkin projection of $\bm{u}(t_{m},\cdot)\in\mathcal{V}$.

∎

###### Lemma E.8 (Local strict convexity for the fine-tuning loss) .

###### Proof.

###### Theorem E.9 (Guaranteed convergence of the fine-tuning) .

###### Proof.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
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?

2. 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?

3. 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?

4. 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?

5. 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?

6. 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?


7. How different are the training and testing initial conditions?
soundness: 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.

2.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.
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.
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.

2. 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.

3. 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.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Can ST-FNO be applied to other equations and still perform effectively?

In line 65, the phrase “Similar to the traditional time marching solvers” is confusing. Did the authors intend to mean "Unlike traditional time marching solvers"?

In lines 217-218, what is the definition of $q$?

In lines 1678-1679, there is a '??'. Please correct this.
soundness: 3
strengths: - The concepts of the spatiotemporal adaptation technique and the hybrid operator learning paradigm appear novel.
- The method is solidly supported by theoretical evidence.
- Overall, the paper is well-written and clear.
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.
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:

  - 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?

  - 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.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
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?

Given the model’s complexity, how interpretable are the results?

How 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?
soundness: 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.
The 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.
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.
weaknesses: As with many neural network-based approaches, the interpretability of the results may be challenging.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: Just some questions:
1. 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?
2. 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?
3. 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?
4. 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?
5. 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?
6. 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?
soundness: 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.
2. 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.
3. 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.
4. Strong empirical results: The method demonstrates significant improvements in accuracy and efficiency on challenging turbulent flow benchmarks.
5. Open-source implementation: The authors provide code to reproduce their results, enhancing reproducibility and potential impact.
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:
1. A new spatiotemporal adaptation (ST-FNO) of Fourier Neural Operators to enable learning maps between Bochner spaces, allowing arbitrary-length temporal predictions.
2. A novel training-fine-tuning paradigm that combines limited epochs of end-to-end training with targeted fine-tuning of a spectral convolution layer.
3. 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.
4. Empirical demonstration of significant improvements in both computational efficiency and accuracy compared to existing methods.
Overall 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.
weaknesses: 1. Limited scope: The method is primarily designed for rectangular domains and uniform grids, which may limit applicability to more complex geometries.
2. Dependency on numerical solvers: The fine-tuning process relies on traditional numerical solvers for computing extra field variables, which may introduce computational overhead.
3. 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.
4. 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.

Meta-review: The work proposes a spatiotemporal adaption to the FNO, combining end to end training with fine-tuning in a negative Sobolev norm. The work is theoretically well-founded and numerical experiments convincingly show the advantages of the method.

Decision: Accept (Poster)