Recently, Deep Learning (DL) pipelines have proven particularly effective in addressing problems formulated by Partial Differential Equations (PDEs). In this paper, we study the problem of learning Neural Operators (NOs) between infinite-dimensional function spaces (Boullé & Townsend, 2023;Kovachki et al., 2023;Azizzadenesheli et al., 2024;de Hoop et al., 2022;2023). Specifically, the problem in consideration is for the Navier-Stokes Equations (NSE) in the turbulent regime (Re = O(10 3 ) to O(10 4 )). In computation, the difficulties of solving NSE in this regime roots from its "stiffness" attributed to the nonlinear convection with a nearly singular viscous diffusion, as well as the spatiotemporal nature of being highly transient. For this problem, we propose to synergize operator learning with classical numerical PDE methods, complementing one's drawbacks and limitations with the other's strengths.
Compromises and drawbacks of traditional numerical methods. To solve NSE, traditional time stepping schemes include Adams-Bashforth/-Moulton or Runge-Kutta (RK) families, e.g., see Canuto et al. (2007, Appendix D), Karniadakis et al. (1991). If one opts to use an explicit scheme, or there exists a certain portion of the forcing terms (e.g., the convection term in NSE) computed via explicit schemes, then extremely small time step sizes must be used: ∆t ∼ O ((∆x) α ), α ≥ 1 (Rannacher, 2000, Chapter 4). This necessity is often referred to as the "stiffness" of the PDE, and the threshold or constraint of the time step is called the Courant-Friedrichs-Lewy (CFL) condition, e.g., see Johnston & Liu (2004); Wang (2012); He & Sun (2007). CFL poses a sufficient condition on the step size for "stability", and this requires the step size usually much smaller than the requirement for "consistency". Here, this temporal consistency usually refers to a first-order optimal local truncation error, e.g., how the original Butcher tableau is derived for RK (Butcher, 1965). The CFL puts a threshold for all explicit time-marching schemes on how fast the local errors in different frequencies can propagate, and thus prevent error accumulation "pollutes" the approximation globally to ensure stability. This means that, for any time-marching schemes, finer mesh (better spatial consistency) requires the time steps to be much smaller to prevent errors from traveling to neighboring nodes and elements, which greatly increases the computational cost.
Spatiotemporal operator learning. Among the end-to-end operator learning for NSE, the common approach is the so-called autoregressive "roll-out". During rolling-out, several snapshots of solutions are concatenated as channels in the input tensor to the neural operator, which outputs an approximated solution at the subsequent time step. This procedure can be repeated recurrently until reaching the model's stability capacity. In contrast to the traditional time marching solvers whose step sizes are restricted by the CFL condition, roll-out can withstand a much bigger time step size. NOs used in rollout approach include the original FNO2d (Li et al., 2020), and those in Li et al. (2022); Brandstetter et al. (2023); Gupta & Brandstetter (2023); Fonseca et al. (2023); Zheng et al. (2024). However, for an end-to-end operator learner, the roll-out approach faces the same dilemma of super-resolution capacity in the temporal dimension as a traditional solver: finer time steps cost prohibitively more, while larger time steps does not guarantee stability. To solve this problem, we turn to the prevalent functional analytic framework for studying the convergence and stability of solution trajectories for NSE: the theory of Bochner spaces, e.g., Aubin-Lions lemma (Lions & Magenes, 2012, Chapter 3), (Temam, 1995, Chapter 2), (Evans, 2022, Chapter 7). Inspired by this theoretical insight, we propose a SpatioTemporal adaption for all Fourier Neural Operator (ST-FNO) variants. ST-FNO now can perform arbitrary zero-shot super-resolution not only spatial dimensions, but also allows the temporal dimension to vary for the first time. For a prototype spatiotemporal operator learner FNO3d in Li et al. (2020), it learns a map between a fixed number of spatial "snapshots" (product spaces). The newly proposed ST-FNO directly learns a map between Bochner spaces L 2 (T 0 ; V) to L 2 (T ; V) on non-overlapping time intervals T 0 and T . This makes the model become "trajectory-to-trajectory", where the initial trajectory is the input of an NO to obtain the output evaluation as an approximation to the subsequent trajectory of the solution. Here V denotes the spatial Hilbert space in which snapshots of the solution at a specific time reside, and for a detailed notation list, we refer the readers to Appendix A.
Limits and the lack of accuracy for neural operators. The NO approach has the potential to bypass various difficulties and compromises of traditional schemes of numerical PDEs. However, in all end-to-end operator learning benchmarks of NSE, even the state-of-the-art NOs still fall short in achieving high-accuracy solutions. For example, an end-to-end roll-out approach suffers from error accumulation experimentally, e.g., see Li et al. (2022, Figure 9). To our best knowledge, in terms of the relative difference with the ground truth under the Bochner norm, no NO-only-based operator learning approach can break the barrier of an even 2-digit accuracy. Moreover, to our best knowledge, NO-only approaches have no theoretical stability estimate, such as the error propagation operator is a contraction. Recently, a remarkable advancement called PDE-Refiner (Lippe et al., 2023) learns an extra error correction NO under the Denoising Diffusion framework (Ho et al., 2020). For a single instance, it can get a stable long roll-out, and achieve for the first time O(10 -8 ) relative difference with the ground truth after a single time marching step, and O(10 -6 ) in the Bochner norm. Nevertheless, for all models above and their learning frameworks, the optimization is to minimize the difference between the outputs from the NO, namely u N , and the ground truth u S , generated by a traditional numerical PDE solver. The difference with the true solution u under a certain norm is not directly optimized, as the analytical expression of u is unknown most of the time in real-life applications. For difficult PDEs such as the NSE, in a linear time-stepping scheme, the ground truth u S may on its own only have 3-digit accuracy in the Bochner norm. Note that this may already Main contributions. The main contributions of this work are summarized as follows. For the difference between the common roll-out approach versus direct spatiotemporal learning between Bochner spaces, we refer the readers to Figure 1.
• Spatiotemporal Fourier Neural Operator. We design the first spatiotemporal adaption technique for all FNO variants (ST-FNO) to enable them to learn maps between Bochner spaces. • New hybrid operator learning paradigm. We propose to train and evaluate ST-FNOs using a new strategy for spatiotemporal predictions, which has significantly improved over the existing methods in speed and accuracy. Only a few epochs of training combined with a spectral-refining finetuner yield highly accurate approximations, extra field variables are obtained by auto-differentiable multistep solvers (such as RK4), the input of which are neural predictions. The auto-differentiability of the solver makes it possible for the fine-tuning to optimize the parameters in FNO. • Functional a posteriori error estimation. Leveraging the spectral structure of the new ST layer in ST-FNO, Parseval identity, and Gelfand's triple, the fine-tuning utilizes a new loss, a posteriori
where u N := u N {θ} nt m=1
Figure 1: Schematic differences between approaches. 4th-order Runge-Kutta (RK4): small time steps bounded by the CFL condition, de-aliasing filter needed; autoregressive NO rolling-out: using previous evaluation as input repetitively, large time steps, no stability guarantees; Spatiotemporal FNO (ST-FNO) with hybrid fine-tuning: large time steps, yielding arbitrary-length temporal prediction in a single evaluation, parallel-in-time optimization for spectral fine-tuning.
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 knowledge of 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 (Kochkov et al., 2021;Dresdner et al., 2023) with enhanced functionality for tensor operations such as the facilitation of fine-tuning for the latent fields, publicly available at https://github.com/scaomath/torch-cfd, with scripts to replicate the experiments as well as the data generation. The data are available at https://huggingface.co/datasets/scaomath/ navier-stokes-dataset .
Both drawbacks and advantages of traditional numerical methods and NO-based methods in Section 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 to modify a generic FNO-based neural architecture to obtain an operator learner between Bochner spaces.
For 2D NSE, in the the velocity-pressure (V-P) formulation (2.2), the velocity field u(t, x) : [0, T ] × Ω → R 2 is seen as an element in the Bochner space L p ([0, T ], V) where V is a spatial Sobolev space in which each snapshot at t of the solution u(t, •) resides. As in Li et al. (2020), we also consider the vorticity-streamfunction (V-S) formulation (2.1) with periodic boundary condition (PBC). In [0, T ] × Ω, the scalar-valued vorticity ω := ∇×u, and streamfunction ψ satisfies u = rot ψ. These two formulations read (Vorticity-Streamfunction)
For all analyses in line with the Hilbertian framework, V = H 1 (T 2 ) for vorticity and H 1 (T 2 ) for velocity, where T 2 denotes the unit torus, i.e., Ω = (0, 1) 2 with a component-wise PBC. For a fixed forcing function in V ′ , the initial condition is either ω(0, x) = ω 0 (x) or u(0,
Here ω 0 is drawn from a compactly supported probability measure µ, in which the compactness corresponds to certain power/enstrophy spectrum decay law to produce isotropic turbulence (McWilliams, 1984), and we refer the reader to Appendix C for details. Then, we can consider the following map
where t ℓ+nt ≤ T . The input and the output in the training data for G µ are snapshots obtained from a solution trajectory with the same initial condition ω 0 that is drawn from certain compactly supported probability measure µ. The operator learning task for NSE using a prototype spatiotemporal operator learner, e.g., FNO3d (Li et al., 2020), is to learn this G between a fixed number of Cartesian products of spatial Sobolev spaces. Specifically,
In the task for this prototype "snapshots learner", the number of snapshots ℓ ∈ Z + is fixed. As such, X µ and Y represent spaces of two non-overlapping temporal segments of solution snapshots. Based on these snapshots' discretized approximations from data, these snapshots learners learn an operator G θ : X → Y, where the number of parameters in G θ is independent of the spatial discretization size, yet does depend on the number of snapshots ℓ. This dependence in the setting of this task makes it not "trajectory-to-trajectory".
In this study, the learning aims to recast (2.3) to a trajectory-to-trajectory operator learning problem, conforming to the Hilbertian formulation of NSE using Bochner spaces. Using the fact that the weak solutions to (2.1) and ( 2.2) at a given time interval T is in L 2 (T ; V), the operator to be learned is
In what follows, we shall present how to modify any common FNO variant to become a Bochner space operator learner that can learn maps from arbitrary-sized discretization in R ℓ×n×n×d to R nt×n×n×d (d = 1 in V-S; d = 2 in V-P). This operator can be trained by a standard end-to-end supervised learning pipeline using lower-resolution data, and for inference, the newly proposed spatiotemporal trajectory-to-trajectory learner, ℓ, n t , n can all be of variable sizes. In contrast, the snapshot learner FNO3d predicts the subsequent n t snapshots with the first ℓ snapshots as input, with both n t and ℓ fixed, and only the spatial resolution n can be much higher than the input-output pairs used in the training. Hereafter we omit the d dimension if no ambiguity arises.
For any FNO variant, such as FNO3d (Li et al., 2020) or Factorized FNO (Tran et al., 2023), with the following meta-architecture:
where P is a lifting operator, Q is a projection operator that does a pointwise reduction in the channel dimension, σ j can be a pointwise universal approximator or simply chosen as a nonlinearity. K j := K ϕj is the parametrized spectral convolution that uses a spatiotemporal 3D FFT. During training, the operator to be learned is restricted to finite-dimensional subspaces X ⊃ X n → Y n ⊂ X , in the sense that the continuous spatial Sobolev space V in the product spaces is replaced by a finite-dimensional subspace V ⊃ S ≃ R n×n with continuous embeddings {a S ∈ R ℓ×n×n } ↩→ X , and {u S ∈ R nt×n×n } ↩→ Y. The positional encodings p := (t i , x j ) ℓ i=1 ∈ R 3×ℓ×n×n represents the underlying spatiotemporal grid, and is concatenated to a S before feeding to P .
For spatiotemporal problems, all FNO variants novelly exploit a convenient architectural advantage of operator-valued NNs: the input temporal dimension ℓ is treated as the channel dimension of an image, thus enabling channel mixing as a learnable temporal extrapolation. However, this neat trick coincidentally makes the lifting operator the biggest hurdle for G θ to become a trajectory-to-trajectory operator learner between Bochner spaces, since the channel dimension of P must be hard-codedfoot_0 and thus depends on the input pair's time steps.
In what follows, we use FNO3d as an example to present three new modifications to FNO3d to become a trajectory-to-trajectory learner ST-FNO3d, that is, to act as an operator that maps an arbitrary-time-step input to an arbitrary-time-step output. We also note that this modification applies to any FNO variant that predicts spatial and temporal dimensions at the same time. These changes are so simple that the trajectory-to-trajectory modification can serve as drop-in replacements for their snapshot learner counterparts, when the temporal input dimensions are fixed in a task. For the schematic difference, please refer to Figure 2, Figure 3, and Appendix D.2 for more comments. : layer normalization to replace a hard-coded global normalization. Combined with channel mixing, the first spectral convolution layer serves as a time-depth-wise separable (global) convolution, after which the time dimension is shrank to a fixed "latent" time dimension through iFFT's resampling. : the spatiotemporal spectral convolution layer as the final layer is fine-tuned after the training phase.
Modifications to 3D FFT. The first change is in the spatiotemporal FFT from an analytical point of view. Using the spectral convolution K(•) in FNO3d as an example, with a slight abuse of notation, denote the latent dimension (width) by d v , let the vector-valued latent representation with d v channels be continuously embedded into Π nt i=1 V by v S ↩→ v in each channel, and denote Λ := [t ℓ+1 , . . . , t ℓ+nt ], then K(•) is "semi-discrete" in a sense as follows
where t takes values only in a discrete set Λ, and x ∈ Ω. Here W ∈ R (dv+1)×dv is a pointwiseapplied affine linear operator, κ ∈ C (Λ × Ω) × (Λ × Ω); R dv×dv , and m denotes an approximation to the Lebesgue measure on Ω. Evolving K into a spatiotemporal spectral convolution acting on Bochner space is straightforward. We change (2.5) slightly by adopting of an atom-like measure δ(•) in the temporal dimension, which then generalizes to the variable-length temporal discretization for any (t, x) ∈ (a, b) × Ω, i.e., one can obtain arbitrarily many snapshots on the interval of interest
(2.6)
During training using the discretized data, each hidden layer is a map from R dv×dt×n×n → R dv×dt×n×n , where d t is a "latent" dimension of time steps that is chosen ≤ n t . The layerwise propagation mechanics remains the same:
, where F and F -1 denote the spatiotemporal Fourier transform and its inverse, respectively. The global spatiotemporal interaction characterized by the kernel is truncated in terms of modes in the frequency domain at
R ϕ is parametrized as a C τmax×kmax×kmax×dv×dv -tensor. Here (τ, k) denotes the coordinate in the spatiotemporal Fourier domain. Note that the integral represented by matrix product with Fκ(•, (s, y)) can then be viewed as residing in the continuous space as the Fourier basis (3.1) can be evaluated at arbitrary point.
Modifications to the lifting operator P . The other major hurdle for FNO3d to become trajectory-totrajectory is that a global normalization is applied with a hard-coded temporal dimension using training data. 2 In ST-FNO3d's modification, a periodic padding along the temporal dimension is applied to the input, and then a time-depth separable convolution layer I with variable time steps is used to map an arbitrary number of snapshots to a fixed number of channels d v with a fixed latent time steps d t . The spatiotemporal positional encodings p S := W p p, with W p a learnable random projection with periodically padded p as input. This makes p S match the latent fields' dimensions such that they can be concatenated. The latent fields, concatenated with p S , are then normalized by a learnable layer normalization layer Ln(•), instead of a global fixed pointwise normalizer for the raw data in FNO3d. Finally, in ST-FNO3d, the lifting operator becomes P := Ln p S + I(•) : R ℓ×n×n → R dv×dt×n×n .
Modifications to the projection operator Q. In FNO3d, the original out projection operator Q maps a tensor in R nt×dv×n×n to another in R nt×n×n . In the spatiotemporal-adapted FNO Q : R dv×dt×n×n → R nt×n×n , which maps the last latent representation to match the dimension of the output u S . In Q, the key for an arbitrary-length inference is that we compose another spectral convolution K S after channel reduction. K S can be thought of as a post-processing layer (also as the only layer to be fine-tuned). It takes advantage of the FFT/iFFT's natural super-resolution capacity, especially in the temporal dimension, by zero-padding the latent temporal step dimension (d t ) to given arbitrary output time steps. For the necessity of this padding, we refer the reader to Figure 7 for an illustrative example. For the V-P formulation, to impose the divergence-free condition for, S is implemented as an optional non-learnable Helmholtz decomposition layer in Q.
Built upon a reasonably accurate approximation, the fine-tuning of ST-FNO is proposed. It is able to yield accuracy on par with traditional numerical methods on the same time horizon, and computational resources used are comparable to the evaluation of NOs. Taking advantage of the efficient ST-FNO zero-shot arbitrary-length temporal inference, this new approach does not need thousands of marching steps like traditional methods. Meanwhile, it borrows the wisdom from traditional Galerkin methods to improve the accuracy (consistency) of the NO approach, liberating the scheme from trade-offs that the traditional methods must make to ensure stability and convergence.
We shall present the proposed fine-tuning using the V-P formulation in subsequent subsections. Denote an ST-FNO evaluation at a specific
where k := 2π(k j ) j=1,2 and τ m := 2πm/(Tt l ). Then, a temporal continuous approximation u N := u N (t, •) can be naturally defined by allowing t vary continuously on T := [t ℓ+1t p , T + t p ] thanks to the spectral basis of S, where t p is the temporal periodic padding in Section 2.2. Define residual functional R(u
2) is violated by the finite-dimensional approximation u N , not in a localized/pointwise fashion, but rather in a global way by representing the error based on its inner product against arbitrary v. At a specific time t, the functional norm of R(u N )(t, •) ∈ V ′ defined as follows is then an excellent measure of the error to be a candidate for a loss function in view of Theorem 3.1:
(3.3) Theorem 3.1 (A posteriori error bound for any fine-tuned approximations, informal version). Let the weak solution to (2.2) be u ∈ L 2 (T ; V), and ∂ t u ∈ L 2 (T ; V ′ ). For u that is sufficiently regular, the dual norm of the residual is efficient to estimate the error for any u N :
Published as a conference paper at ICLR 2025 Moreover, if u and u N are "sufficiently close", then it is reliable to serve as an error measure:
(3.5) where T m := (t m , t m+1 ], and the constants depend on the regularity of the true solution u.
The functional norm (3.3) is "global" because it does not have a natural ℓ 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 Peetre (1975) and Gel'fand & Shilov (2016, Chapter 3)), we can define the pairing between V and V ′ as follows without getting into the intricate natures of the tempered distribution:
The spatial Sobolev space H s (T 2 ) can be alternatively identified using norm ∥ • ∥ s and seminorm | • | s (e.g., see Ruzhansky & Turunen (2009, Def. 3.2.2)) as follows for any
Moreover, we have the subsequent lemma in our specific case.
Spatially, we realize a regularized negative Sobolev norm by the Fast Fourier Transform (FFT):
, where Z 2 n := (Z/nZ) 2 , and α ≥ 0.
(3.8)
With this norm at hand, (3.3) and the Bochner norm in (3.4) of the residual are realized to serve as the loss function in the fine-tuning to "refine" the approximation in the spectral domain
(3.9)
Parseval identity. In the context of using an optimization algorithm to train an NN as a function approximator, it is known (e.g., Siegel et al. (2023)) that whether the integral in the loss function is accurately computed affects whether the NN can achieve the scientific computing level of accuracy. For example, on a uniform grid, the accuracy of the mesh-weighted spatial ℓ 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 (Trefethen & Weideman, 2014, 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, in Li et al. (2024d), 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 Li et al. (2022). To our best knowledge, the H -1 -functional norm has not been applied in either function learning or operator learning problems. In fact, the relation of the Gelfand triple is so simple and elegant, leading to Theorem 3.2, that the functional norm is nothing but an inverse frequency weight 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 Lippe et al. (2023) 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 8 and Figure 6. 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.
With the trajectory-to-trajectory learner and the error estimators, we propose a new training-finetuning paradigm. Another important motivation of this new paradigm is that, experimentally, all FNO variants capture the statistical properties of the 2D turbulence (Benzi et al., 1987;Boffetta & Ecke, 2012) quickly in training. Even after a single epoch, the evaluations of ST-FNO converge to a reasonably small neighborhood of the ground truth in terms of the frequency signature of the energy cascade of the flow (see e.g., Figure 8, and Appendix D.3 for more details). Then, 95% of the FLOPs spent in training only contribute to marginal improvements. These observations motivate 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 (Lippe et al., 2023) 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 (3.2), one needs to compute ∂ t u N . Instead of a common approach of using NO to represent the extra field variables (Wen et al., 2022;Brandstetter et al., 2023), we opt to apply a traditional numerical solver B γ (•) with an extremely fine time step, e.g., O((∆t) γ ) for γ ≥ 2, to compute these extra field variables (Line 8 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 , . . . , t ℓ } aiming to predict the trajectories at {t ℓ+1 , . . . , t ℓ+nt }, the new paradigm to train and fine-tune ST-FNO is outlined in Algorithm 1.
Algorithm 1 The new parallel-in-time spectral refining fine-tune strategy in small data regime.
time stepping scheme B γ (•); optimizer D(θ, ∇ θ (•)); training dataset: solution trajectories at [t 1 , . . . , t ℓ ′ ] as input and at [t ℓ ′ +1 , . . . , t ℓ ′ +n ′ t ] as output. 1: Train the ST-FNO model until the energy signature matches the energy cascade. 2: Freeze Θ in G Θ of ST-FNO G θ,Θ , keep Q θ trainable. 3: Cast all nn.Module involved and tensors to torch.float64 and torch.complex128 hereafter. Input: Evaluation dataset: solution trajectories at [t 1 , . . . , t ℓ ] as input, output time step n t . 4:
Extract the latent fields v (m+1) N output of G Θ at t m+1 and hold them fixed. 6: By construction of ST-FNO:
March one step with (∆t) γ using B γ and approximate the ∂ t u N as follows:
) for all evaluation time steps.
10:
Apply the optimizer to update parameters in Q: θ ← D(θ, ∇ θ (η 2 )).
Forward pass only through Q to update u
) for all m. Output: A sequence of velocity profiles at corresponding time steps {u
Interpretations of the theoretical results. Theorem 3.1 establishes that the functional norm of the residual R(u N ), without accessing u, is a good representation of the error uu N in Bochner norms. Estimate (3.4) indicates that reducing the a posteriori error estimator is a necessary condition for reducing the true error. While estimate (3.5) is more delicate in that it is only reliable when u N gets "close" to u. Theorem 3.2 lays the foundation to accurately evaluate this functional norm via FFT. Nevertheless, (3.5) serves as a guideline to design the "refining" procedure (lines 11 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.
In this illustrative example, we examine the 2D Taylor-Green vortex model (Taylor & Green, 1937), whose analytical solution is known and frequently employed as a benchmark for evaluating traditional numerical schemes for NSE (Gottlieb & Orszag, 1977). 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 FNO3d used in this example is scaled down to 1 layer.
This meta-example contains a series of examples of isotropic turbulence (McWilliams, 1984), featured in various work such as Kochkov et al. (2021); Li et al. (2020). The energy and the enstrophy spectra satisfy the energy law of turbulence first proposed by A. N. Kolmogorov (Kolmogorov, 1941). The data are generated using a second-order time-stepping scheme that is proven in theory to preserve the inverse cascade of the energy/enstrophy spectra (Wang, 2012;Gottlieb et al., 2012). We consider two cases:
(I) NSE benchmark in Li et al. (2020), ν = 10 -3 , ω 0 ∼ N (0, (-∆ + τ 2 I) -α/2 ), the energy density in wavenumber k
, no drag, fixed force with a low wavenumber. (II) The famous decaying turbulence discovered in McWilliams (1984), the initial power spectrum Tables 2 and 3 report the results for example (I) and example (II), respectively. In Appendix D.2,Table 5 and Table 6 report the architectural changes and the computational efficiency comparison with roll-out and traditional solver.
We designed a new pipeline to train and fine-tune a new spatiotemporal modification of FNO to get close to the true solution (not the ground truth generated by the numerical solver) of NSE in the turbulent regime. The evaluation errors in benchmark problems are up to 10 5 times better than the non-fine-tuned FNO variants trained by a simple end-to-end pipeline. Due the exploitation of the intricate connections with traditional spectral methods, e.g., the optimally learned parameters correspond to a Fourier-Galerkin projection (Kovachki et al., 2021) using the Fourier basis, only FNO-based neural operators benefit from this advancement. How to generalize this new pipeline to other types of operator learners in Appendix B will be worthy of exploration.
Table 1: Results for Taylor-Green vortex, the relative errorsa at the final time step, ε := ωtrue -ωN .
After fine-tuning
Table 3: Evaluation metrics of the McWilliams isotropic turbulence example. All models are trained using on 64×64 mesh and evaluated on 256×256 mesh. ε := u S -u N or ω S -ω N . S stands for the Fourier approximation space (3.1) on a 256 × 256 fine grid.
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 10.
After fine-tuning
A NOTATIONS T 2 := S 1 × S 1 where S 1 is homeomorphic to [0, 1) with the point 1 being 0.
the Bochner spaces containing {u :
Bochner space containing {u : T → H ess sup t∈T ∥u(t, •)∥ V < +∞} u the true weak solution to the NSE u(t, •) ∈ V for any t u (l) S the ground truth data generated at t l by the numerical solver in S ⊂ V u (m) N neural operator evaluations at t m that can be embedded in
both are viewed as column vectors B DETAILED LITERATURE REVIEW AND MOTIVATIONS Interplay of deep learning and PDEs. PDE solvers are function learners to represent PDE solutions using neural networks (Han et al., 2018; Raissi et al., 2019; Chen & Koohy, 2024). PDE discovery encompasses all the techniques dedicated to identifying and optimizing PDE coefficients from data (Brunton et al., 2016; Champion et al., 2019). Recently, hybrid solver approach has been explored in Greenfeld et al. (2019); Kochkov et al. (2021); Huang et al. (2023); Taghibakhshi et al. (2023; 2021); Hu & Jin (2024); Huang et al. (2022). Reinforcement learning has been applied in the field of mesh generation to build more efficient solving pipelines (Yang et al., 2023;Gillette et al., 2024). For neural operators, remarkable outcomes are achieved in diverse applications, such as weather forecasting (Keisler, 2022) and turbulent fluids simulations (Shu et al., 2023;Li et al., 2023a;Lienen et al., 2023), 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 (Zheng et al., 2023); and the ability to provide mesh-free and resolution-independent solvers in cases of irregular domains (Hao et al., 2023;Li et al., 2024c).
Neural operators. Looking at neural operator architectures, we can identify two different 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) (Lu et al., 2021;Goswami et al., 2023;Wang et al., 2021), wherein aggregation occurs via linear combination. Similarly, Fourier Neural Operator (FNO) (Li et al., 2020), and some of its variants (Gupta et al., 2021;Rahman et al., 2023;Tran et al., 2023;Li et al., 2024c;Wen et al., 2022), aggregate different latent representations through convolution with a learnable kernel in the frequency domain. Additional examples with a linear aggregation with a U-Net meta-architecture can be found in Raonic et al. (2024); He et al. (2024). Bartolucci et al. (2024) further explores the error analyses for the linear aggregations through the lens of the frame(let) theory. The second 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 (Vaswani et al., 2017), which builds an instance-dependent kernel.
In the context of operator learning applications, Transformer-based neural operatos have been studied in Cao ( 2021 Numerical methods for NSE. The NSE falls into the category of a stiff PDE (system) ∂ t u = Lu + N (u) + f , where f is the external forcing, L and N (•) are a linear and a nonlinear operator, respectively. Trefethen noted back in Kassam & Trefethen (2005) on the difficulty to design a time-stepping scheme as N (•) and L have to be treated differently. There are a long history of numerical methods for NSE we draw inspiration from. Petrov-Galerkin methods have been developed for NSE in Boffi et al. (2013); Girault & Raviart (2012). Nonlinear Galerkin method (Marion & Temam, 1990) inspires us to prove Theorem E.9. Chorin (1968) uses a clever trick to impose the divergence free condition without constructing a divergence-conforming finite element subspace. Shen (1994) designed various Galerkin methods in the space of orthogonal polynomials. Pioneered by Chorin, Shen, and E, projection methods are among the most popular schemes to solve NSE (see e.g., Weinan & Liu (1995) for a summary), also serves as an inspiration to add the Helmholtz layer. Bernardi et al. (1992) proposed a mixed discretization for the vorticity-streamfunction formulation.
Consistency-stability trade-offs. In view of the Lax equivalence principle ("consistency" + "stability" =⇒ "convergence" (Lax, 2002, 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 (Lele, 1992), yet the lack of stability is more severe and needs higher-order temporal smoothing. On the other hand, implicit schemes can be unconditionally stable for stiff or highly transient NSE with relatively large time steps 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. One famous example of the consistencystability trade-off is the 3/2-rule (also known as 2/3-dealiasing dealiasing filter) for the nonlinear convective term (Orszag, 1971b;Patterson & Orszag, 1971;Hou, 2009;Gottlieb & Orszag, 1977) for pseudo-spectral methods (Orszag, 1971a;1972). The highest 1/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. Numerical results suggest that for pseudo-spectral spatial discretization with no higher-order Fourier smoothing temporally, the dealiasing filter is indispensable (Tadmor, 1987), as the time marching may experience numerical instability (Kreiss & Oliger, 1979;Goodman et al., 1994) without it. This is due to the nonlinear interaction in the convective term, which is caused by the amplification of high-frequency "aliasing" errors when the underlying solution lacks sufficient smoothness. In this study, the hybrid approach we adopted combines the strengths of NOs and traditional numerical solvers. There is 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.
Why and why not functional-type a posteriori error estimation? For the sake of presentation, we make some handy assumptions: (1) A : V → V ′ denotes a linear operator on the Gelfand triple; (2) u S is obtained through Galerkin methods (Evans, 2022), which is a projection onto a finite-dimensional approximation subspace of S ⊂ V; and (3) an evaluation of u N can be continuously embedded into S. Then, we have the following enlightening Pythagorean-type identity for the true solution u
Through the Riesz representation theorem, u S is solved through (u-u S , v) a := ⟨A(u)-A(u S ), v⟩ = 0, for any v ∈ S. This bilinear form (•, •) a induces a (semi)norm ∥ • ∥ a inheriting the topology from V. Given this orthogonality, minimizing the difference between u N and u S becomes fruitless if ∥uu S ∥ a is relatively big in the first place, and unnecessary computational resources may have been spent to get closer to u S . Rather, (B.1) indicates that it is more efficient if one can design a method to reduce the error of ∥uu N ∥ a directly, while circumventing the fact that u is not accessible.
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 (Veeser & Verfürth, 2012) or inverse inequalities (Carstensen & Funken, 1999;Veeser & Verfürth, 2009), see also Verfürth (2013, §1.6.2). Functional-type a posteriori error estimation (Repin, 2008) consider the error as a functional, which is equivalent to error to bilinear form-associated norm as follows
where the weak form of the PDE is
The common approach is using the help from extra "flux" or "stress" dual variables (Repin, 2008, §6.4), see also Repin (2000) for how to get an accurate representation of the error functional in L 2 . For example, for the Stokes flow, which can be viewed as a steady-state limit of viscous fluid, such that ∂ t u = 0 and no nonlinear convection in (2.2), Re ∼ O(1), Repin (2008, §6.2) estimates error of an H(div)-L 2 mixed discretization as follows
where C P is the Poincaré constant of the compact embedding and (σ, 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 (σ, 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 ∆u 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 ∇ • σ where σ ∈ H(div), and is constructed to be closer to the true solution's gradient ∇u than ∇u 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 (Girault & Raviart, 2012, Chapter III §1 Sec 1.2). Recently, Fanaskov et al. (2024) 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 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.
Error correction for Bayesian inverse problems. The a posteriori error-correction approaches through sampling for the output of surrogate NNs in Bayesian inverse problems (Yan & Zhou, 2020). More recently, in Cao et al. (2023), the error equation in an inverse problem is approximated by solving for a Galerkin projection of a linearized error equation, while in this paper, we leverage the Gelfand triple directly to compute a nonlinear Galerkin projection directly by minimizing the spectral norm.
In the context of function learning problems, Du et al. (2024) applied an L 2 -spectral loss for Physics-informed Neural Networks (PINN) (Raissi et al., 2019) and found superior results over mesh-weighted spatial losses. Du et al. (2024) realized the PDE residual as a functional on L 2 and exploited the Parseval identity on this space, which computes sup v∈L 2 ⟨R, v⟩/∥v∥ L 2 . We note that this is not the natural space to measure the PDE residual according to the Galerkin formulation of NSE. Using the Poisson equation -∆u = f as a simpler example, if u ∈ H 1 , then the residual f + ∆u ∈ H -1 . -∆ is a bounded operator from H 1 to H -1 which guarantees the stability estimate
The dual space of H 1 (T 2 ) is the natural space to consider the residual functional. The reason is that the solution to NSE belongs to H 1 (T 2 ) at a given snapshot. The dual space of H 1 (T 2 ) happens to be H -1 using the Gelfand triple. In this paper, we compute sup v∈H 1 ⟨R, v⟩/∥v∥ H 1 . However, this H -1 -norm is a non-localizable functional norm in the spatial domain, and traditional wisdom (FEM or FVM) has to circumvent the direct evaluation of it by localizing the residual, which introduced inaccurate compromises in the residual-based error estimation.
There is a vast amount of literature for the a posteriori error estimation for the viscous flow (1/Re > 0) problems under the Hilbertian framework for traditional numerical methods. Our methods draw inspiration from these pioneers and try to address the drawbacks. The most popular a posteriori error estimator for stationary Stokes problem is from Verfürth (1989), and it is of residual-type by computing a mesh-weighted L 2 -norm elementwise, and the singular distribution ∆u S is represented by the magnitude of flux jumps across the facets in this mesh. In Bank & Welfert (1991), 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 Oden et al. (1994). 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 (Um et al., 2020); interpolation to achieve spatial super-resolution using solutions produced by traditional numerical solvers marching on coarser grids (Kochkov et al., 2021). Sun et al. (2023) uses orthogonal polynomial features (HiPPO from Gu et al. (2020)) to achieve temporal super-resolution in addition to the spatial one. McGreivy & Hakim (2023) explores how to incorporate energy-preserving schemes such as upwinding into the ML solvers through traditional numerical methods. Dresdner et al. (2023) novelly explores the possibility of using neural networks to correct the result of the traditional time-marching solvers. Qi & Sun (2024) identifies important frequency bands in FNO parametrization and improves FNO's rolling-out performance in NSE benchmarks.
Here we derive the vorticity-streamfunction formulation for the convenience of the readers. Consider the standard NSE in 3D, u's z-component is 0, and u = u(t, x, y) has only planar dependence
Taking ∇×(•) 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 ω = ∇ × u Next, one can have ω • ∇u = 0, thus the equation becomes
We need to introduce a streamfunction ψ to recover the velocity u = ∇ × (0, 0, ψ). The main simplification comes from the assumption of dependence on xand y-variable only. Thus, in 2D this becomes u = rot ψ (a π/2 rotation of the 2D gradient ∇ψ). Therefore,
and the vector vorticity can be equivalently represented by a scalar vorticity ω ω = ∇×u = (0, 0, curl(u 1 , u 2 )) = (0, 0, ∂ x u 2 -∂ y u 1 ) =: (0, 0, ω)
The equation becomes a nonlinear system for ω and ψ:
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π) 2 such that the inflow/outflow occurs on the "boundary" (see Figure 4). The exact solution is given by: u(t, x, y) = e -2κ 2 νt sin(κx) cos(κy) cos(κx) sin(κy) or e -2κ 2 νtcos(κx) sin(κy) sin(κx) cos(κy),
where ν = 1/Re is chosen as 10 -3 . For a sample trajectory with κ = 1, please refer to Figure 4. 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 κ ranging from 1 to 11, among which κ = 11 is chosen as the test trajectory.
The example featured in the original FNO paper (Li et al., 2020) can be viewed as a special case of the isotropic turbulence with a strongly regularized initial condition, no drag, and a force with low wavenumber. This is our example (I). The example (II) has weakly-regularized initial condition, small drag, and no external force. The initial condition spectra are chosen such that, after the warmup time, the spectra are the comparable with the Kolmogorov flow featured in Kochkov et al. (2021). The similarity is in the sense of the energy cascade formed in the Fourier domain (see Figure 8 and 9). The training data used in both (I) and (II) are generated using a pseudo-spectral method with a de-aliasing filter using an RK4 marching for the convection term and implicit Euler for the diffusion term. The Reynolds number used in this example is 1000. The time steps for both are ∆t = 10 -3 on a 256 2 grid, then downsampled to a 64 2 grid. The warmup time is 4.5. The δt for the ST-FNO training and evaluations are δt = 55∆t. The testing data are generated on a 512 2 grid with a 5 × 10 -4 time step then downsampled to 256 2 . There are 1024 trajectories in the training data. The input of model has 10 snapshots from t = 4.5 to t = 4.5 + 9 δt. The output ground truth has the 10 subsequent snapshots from t = 4. 0 20 40 60 0 10 20 30 -20 -10 0 (a) 0 20 40 60 0 10 20 30 -20 -10 0 (b) Figure 6: Pointwise values of residual of ST-FNO in the frequency domain in ln scale of Example (I). (a): the residual of ST-FNO before correction, time derivative computed using the ground truth; (b): the residual of ST-FNO after fine-tuning for 10 ADAM iterations where the time derivative is computed using an extra-fine-step numerical solver.
The training uses torch.optim.OneCycleLR (Smith & Topin, 2019) 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 • lr max . The lr max = 10 -3 . Despite using the negative Sobolev norm is quite Figure 7: A simple illustration of necessity of padding using a 1D non-periodic data. FNO3d cannot handle temporal non-periodic data. The free temporal super-resolution by FFT/iFFT has Runge-like phenomena near the end points without padding. Padding in the ST-FNO3d modification resolves this problem. Since data themselves are intrinsically lower-dimensional, we opt not to pad the latent representations, but rather apply the padding only in the last spectral convolution layer to get the output.
efficient in fine-tuning, using it in training not be the most efficient in minimizing the L 2 -norm due to the optimization being non-convex, we observe that all norms are equivalent "bad" due to nonlinearity of NOs. The result demonstrated is obtained from fixing the random number generator seed to be 1127825. All models are trained on a single RTX A4500 or RTX A6000, the inference is done with a power limiter on the GPU to ensure fairness. The codes to replicate the experiments are open-source and publicly available at https://github.com/scaomath/torch-cfd.
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 W ⊂ 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 (Girault & Raviart, 2012, Chapter III 1 Section 1.1). Inspired by the postprocessing to eliminate ∇p (Ku et al., 1987, eq. ( 15)) together with the neural Clifford layers (Brandstetter et al., 2023), we add a Helmholtz layer S after each application of σ j • (W j + K j ). S performs a discrete Helmholtz decomposition (Girault & Raviart, 2012, Chapter 1 §1 Section 3.1) in the frequency domain to project the latent fields to be solenoidal.
Difference with Fourier Neural Operator 3d In view of (2.5), the FFTs in FNO3d transforms are continuous integrals in the spatial dimensions yet a discrete sum in the temporal dimension due to the channel treatment. In the original FNO3d architecture, despite that no explicit restriction imposed on the output time steps, the d out cannot be trivially changed as the data are prepared by applying a pointwise Gaussian normalizer that depends on d out . Denote the lifting operator (channel expansion) in FNO3d by P : {a h } ⊕ {p h } → R dv×dout×n×n . Before the application of P in FNO3d, a h is artificially repeated d out times, then P is in a space of linear operators ≃ R (din+3)×dv . This observation makes it independent of the spatial discretization size n × n by dependent on d in . In FNO3d, the channel reduction is a pointwise nonlinear universal approximator, yet in ST-FNO3d, this is linear, which eventually facilitate the convexity of the fine-tuning optimization problem.
In implementations, we choose Iter max = 100 and the ADAM optimizer with learning rate 10 -1 in Algorithm 1, i.e., an ADAM is simply run for 100 iterations to parameters that count only a fraction of a spectral layer since the number parameters accounts for a single channel. To train or to stop? Investigating the frequency domain signatures In Stuart's paper on the convergence of linear operator learning (de Hoop et al., 2023), the evaluation error scales with number of samples, in the few-data regime, the overfitting comes fast. Part of the reason is that the operator learning problems are in the "small" data regime.
Denote ω(t, x) := ∇ × u, and ω := Fω where F is applied only in space, then we can define enstrophy spectra as follows:
In Example (II), the enstrophy spectrum decays as O(k -3 ), which is slightly faster than kinetic energy. This is usually referred to as the direct cascade, as opposed to the inverse cascade O(k -5/3 ) discovered by Kolmogorov. For ST-FNO, it learns the frequency signature of the data after a single epoch. See Figure 8 and Figure 9.
In this subsection, we perform the following additional experiments regarding the possibility of fine-tuning ST-FNO3d for "out-of-distribution" data. All ST-FNO3d models are trained for 15 epochs (5 more than the model in Table 3) with the data generated using the McWilliams's isotropic turbulence data from Example (II) with ν = 10 -3 . The initial energy distribution satisfies:
The parameter k 0 in Example (II) roughly means that the initial energy is concentrated at |k| ≈ k 0 .
In Table 3 k 0 is chosen as 4 according to McWilliams (1984). For all examples, the same "burn-in" time is used (t = 4.5), the training output time steps are 10 and evaluation output time steps are 40. Please refer to Figure 10 for sample trajectories and Figure 11 for the energy cascade after the "burn-in" time period. (Ex1) The initial energy is concentrated at k 0 = 2. A sample trajectory can be found in Figure 10(a). (Ex2) Same with above, k 0 = 8. A sample trajectory can be found in Figure 10(c).
(Ex3) k 0 = 4, ν = 1/Re = 1/5000. In order to resolve the small viscosity, the mesh size becomes 512 × 512.
(Ex4) k 0 = 8, ν = 1/Re = 1/10000, mesh size is 1024 × 1024 for the same reason as above.
The evaluation and fine-tuning result can be found in Table 7. We found that the new paradigm works reasonably well for first three example, yet fails to converge for the last example. For an example plot of the prediction comparison please refer to Figure 12. The correlation graph between the ST-FNO3d's predictions and the ground truth can be found in Figure 13, in which the FNO(2+1)d is adapted from the original FNO2d code. This roll-out-based model predicts the next snapshot based on 10 given snapshots as input, and the predicted snapshot is concatenated to the current input as the latest snapshot to form the input for the next prediction.
Table 7: Evaluation metrics of the McWilliams isotropic turbulence example for "out-of-distribution" data. For details please refer to Section D.4 ε := ω Sω N . r := R(ω N ). Y := H -1 (T 2 ). All error norms are evaluated at the final time step. The training is under the L 2 -norm and the fine-tuning is under the H -1 -norm. Symbol "-" means the errors are the same as the entry on the row above.
In-distribution evaluation Out-of-distribution evaluation After fine-tuning
Assumption E.1 (Assumptions for Theorems 3.1, 3.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 V ⋐ H ⋐ V ′ hold, where V = H 1 (T 2 ), and
The time interval 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 u 0 ∈ V, and f ∈ L 2 (T ; V ′ ). Lemma E.2 (Skew-symmetry of the trilinear term). For z, u, v ∈ H 1 (T 2 ), where T 2 denotes Ω := (0, 1) 2 equipped with component-wise PBC. Denote
Proof. We note that this result is common in NSE literature, see e.g., Temam (1995, Part I Sec 2.3), for homogeneous boundary or whole space. In this case, it suffices to verify that the boundary term (a) (b) (c) Figure 11: Energy spectrum densities datasets for additional experiments for Example (II). (a): k 0 = 2; (b): k 0 = 4. (c): k 0 = 8.
vanishes. Here we still provide a short argument for the convenience of readers. First by the product rule, and a vector calculus identity of ∇ • (v • u), see e.g., Balanis (2012, Appendix II.3.2).
By Gauss divergence theorem, we then obtain:
with n is outer normal with respect to ∂Ω, and n i denotes that with respect to edge e i . Now on e 1 := {x = 0} × {y ∈ (0, 1)}, n 1 = (-1, 0) ⊤ , w| e1 = w| e2 by PBC where e 2 := {x = 1} × {y ∈ (0, 1)} and w ∈ {z, u, v}. The integrals on e 1 and e 2 cancel with one another since n 2 = -n 1 . Furthermore, if ∇ • z = 0, integrating (E.2) on Ω yields the desired result. Lemma E.3 (Poincaré inequality). V/R ⋐ H/R is compact and the following Poincaré inequality holds with constant 1: for any v ∈ V/R ∥v∥ ≤ |v| 1 (E.3)
Proof. A common textbook proof exploits the equivalence of sequential compactness and compactness, once the norm topology is we prove this in a very intuitive way under our special setting: by
E.1 PROOF OF THEOREM 3.1
To prove Theorem 3.1, we need Lemmas E.4, E.5, and E.6, all of which can be found in classical textbook on analysis and approximations of NSE, such as Temam (Temam, 1995), or Girault & Raviart (Girault & Raviart, 2012) for homogeneous Dirichlet boundary condition for the velocity 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 H 1 -norm on V = H 1 (T 2 ) for α > 0 on Ω is defined as
1/2 , and |v| α,1 := ∥α∇v∥ L 2 .
For f ∈ V ′ , define a weighted dual norm as follows:
Lemma E.4 (Contuinity and embedding results for the convection term). For trilinear convection term
Proof. The first result is obtained by the Hölder inequality, and the second holds with suite Sobolev embedding results (Temam, 1995, Part I, §2.3, Lemma 2.1).
Lemma E.5 (Energy stability of NSE in a bounded domain).
∥u∥ 2 L ∞ (T ;H) ≤ ∥u 0 ∥ 2 H + ∥f ∥ 2 L 2 (T ;V ′ ) .
(E.6)
Proof. (E.17) is a classical result, see e.g., Temam (1995, Section 3.1).
Lemma E.6 (Fréchet derivative of the convection term). Given v, u ∈ V, the linearization of the difference of the convection term reads Proof. This result is normally used without proof in linearizing NSE, see e.g., Girault & Raviart (2012, Chapter 4 eq. (6.5.2)). Define F (u, G) u ⊤ G, then expanding F at v and G = ∇u yields the result.
Theorem Part I (Functional-type a posterior estimate is efficient (rigorous version)).
Consider the Gelfand triple V ⋐ ⋐ V ′ , and the weak solution u ∈ L 2 (T ; V) ∩ L ∞ (T ; H) to (2.2) be sufficiently regular, then the dual norm of the residual is efficient to estimate the error:
The constant depend on the regularity of the true solution u.
Proof. The proof of the lower bound (efficiency) (E.8) follows the skeleton laid out in Verfürth (2003, Lemma 4.1) for diffusion and temporal derivative terms, while not needing to extend approximation using an affine linear extension operator as Verfürth (2003) does. The treatment of the nonlinear convection term follows Fischer (2015, Lemma 7) for the steady-state NSE (no temporal derivatives), barring the technicality of the divergence-free condition. Consider at t ∈ T , for any test function v ∈ H 1 (T 2 ) and ∇ • v = 0, integrating v against R(u N ) in (3.2) yields
Integrating by parts for the diffusion term, using the same argument for PBC as in Lemma E. .
(E.9) not be the temporal grids t m or t m+1 are naturally defined using the basis in (3.1). Using the error representation in (E.9), we have the (I) and (II) terms in ⟨R(u N ), uu N ⟩ are
As a result, we have
Integrating from t m to any t ∈ (t m , t m+1 ] yields
Since t ∈ (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, norm of the error ≃ (I) + (II) = ⟨R(u N ), uu N ⟩ -( * ) ̸ ≤ ⟨R(u N ), uu N ⟩
We note that Assumption (2.1.2) is equivalent to putting a threshold on ∥u(t, •)∥ L ∞ , 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 γ in Assumption (2.1.1) Using the "sufficient close" assumption (E.10) on the nonlinear term, we have
Now if ν := νγ ≥ 0, one would reach a similar estimate as (E.17) by replacing ν with ν. If no extra assumption on γ is imposed.
Another way is to consider a Sobolev embedding directly after applying the Hölder inequality, for (E.15) we have
where ρ is the constant in the following Sobolev embedding ∥uu N ∥ L 4 ≤ ρ∥∇(uu N )∥.
Now if ν := ν -ρ∥∇u∥ ≥ 0, one would reach a similar estimate as (E.17) by replacing ν with ν. 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 ∈ H/R with its Fourier series in S ∞ ,
Note since f lies in the quotient space, f (0, 0) = 0. Then, the duality pairing can be identified as an L 2 -inner product, as well as the test function v ∈ H, ⟨f, v⟩
As a result, by the definition in (3.7) and the estimate above
Line 97 in the original FNO3d source code fourier_3d.py (link provided for an unaltered fork of the master branch commit de514f2 in the original FNO GitHub repository).
Line 205 and 209 in the original FNO3d source code fourier_3d.py (link provided for an unaltered fork of the master branch commit de514f2 in the original FNO GitHub repository).