Learning to correct spectral methods for simulating turbulent flows
Authors: Gideon Dresdner
Dmitrii Kochkov
Peter Christian Norgaard
Leonardo Zepeda-Nunez
Jamie Smith
Michael Brenner
Stephan Hoyer
Venue: Accepted by TMLR ()

Abstract: Despite their ubiquity throughout science and engineering, only a handful of partial differential equations (PDEs) have analytical, or closed-form solutions. This motivates a vast amount of classical work on numerical simulation of PDEs and more recently, a whirlwind of research into data-driven techniques leveraging machine learning (ML). A recent line of work indicates that a hybrid of classical numerical techniques and machine learning can offer significant improvements over either approach alone. In this work, we show that the choice of the numerical scheme is crucial when incorporating physics-based priors. We build upon Fourier-based spectral methods, which are known to be more efficient than other numerical schemes for simulating PDEs with smooth and periodic solutions. Specifically, we develop ML-augmented spectral solvers for three common PDEs of fluid dynamics. Our models are more accurate (2-4x) than standard spectral solvers at the same resolution but have longer overall runtimes (~2x), due to the additional runtime cost of the neural network component. We also demonstrate a handful of key design principles for combining machine learning and numerical methods for solving PDEs.


Published in Transactions on Machine Learning Research (04/2023)
Learning to correct spectral methods for simulating turbulent flows
Gideon Dresdner
gideond@gmail.o com
Google Research and ETH Zurich, Department for Computer Science
Dmitrii Kochkov
dkochkov@@google..cm
Google Research
Peter Norgaard Google Research
pnorgaardagoogle.com
Leonardo Zepeda-Núñez Google Research
beepedamumezz@google..oom
Jamie A. Smith Google Research
jamiens@google.com
Michael P. Brenner Google Research
mbrennerr@google.com
Stephan Hoyer
shoyer@gooole.ccom
Google Research
Reviewed OpenReview: https: openreview. net/ porm?? id= Jn
Abstract
Despite their ubiquity throughout science and engineering, only handful of partial differ- ential equations (PDEs) have analytical, or closed-form solutions. This motivates vast amount classical work on numerical simulation of PDEs and more recently, a whirlwind of research into data-driven techniques leveraging machine learning (MT). recent line of work indicates that hybrid of classical numerical techniques and machine learning can offer significant improvements over either approach alone. In this work, show that the choice of the numerical scheme is crucial when incorporating physics-based priors. We build upon Fourier-based spectral methods, which are known to be more efficient than other numerical schemes for simulating PDEs with smooth and periodic solutions. Specifically, we develop ML-augmented spectral solvers for three common PDEs fluid dynamics. Our models are more accurate (2 4x) than standard spectral solvers at the same resolution but have longer overall runtimes (~~ 2x), due the additional runtime cost of the neural network component. We also demonstrate handful of key design principles for combining machine learning and numerical methods for solving PDEs.
1 Introduction
The numerical simulation of nonlincar partial differential equations (PDEs) permcates science and engineer- ing, from the prediction of weather and climate to the design of engineering systems. Unfortunately, solving PDEs on the fine grids required for high-fidelity simulations is often infeasible due to its prohibitive compu- tational cost. This leads to inevitable trado-offss between runtime and accuracy. The status quo is solve PDEs on grids that are coarse enough to be computationally feasible but are often too coarse resolve all phenomena of interest. One classical approach is to derive coarse-grained surrogate PDEs such as Reynold's

---

Published in Transactions on Machine Learning Research (04/2023)
Averaged Navier Stokes (RANS) and Large Eddy Simulation (LES) (Pope, 2000), which in principle can be accurately solved on coarse grids. This family of approaches has enjoyed widespread success but is tedious to perform, PDE--ppecffi, and suffers from inherent accuracy limitations (Durbin, 2018; Pope, 2004).
Machine learning (ML) has the potential to overcome many of these limitations by inferring coarse-grained models from high-resolution ground-truth simulation data. Turbulent fluid flow is an application domain that has already reaped some of these benefits. Pure ML methods have achieved impressive results, in terms of accuracy, on diverse set of fluid flow problems (Li et al., 2021; Stachenfeld et al., 2022). Going beyond accuracy, hybrid methods have combined classical numerical simulation with ML to improve stability and generalize to new physical systems, e.g., out of sample distributions with different forcing setups (Kochkov et al., 2021; List et al., 2022).
However, hybrid methods have been limited to low-order fiitte-difference and finite volume methods, with the exception of one recent study (Frezat ct al., 2022). Beyond finite differences and finite volumes, there a broad field of established numerical methods for solving PDEs. In this paper, we focus on spectral meth which are used throughout computational physics (Trefethen, 2000; Burns et al., 2020) and constitute the core of the state-of-the-art weather forecasting system (Roberts et al., 2018). Spectral methods, when applicable, are often preferred over other numerical methods because they can be more accurate for equations with smooth solutions. In fact, their accuracy rivals that of the recent progress made by ML. This begs the question: Can we improve spectral solvers of turbulent fluid flows using learned corrections of course grained simulations?
Our contributions are as follows:
We propose a hybrid physics machine learning method that provides sub-grid corrections to classical spectral methods.
2. We explore two toy 1D problems: the unstable Burgers' equation and the Kuramoto-Sivashinsky the hybrid model is able to make the largest improvements.
3. We compare spectral, finite-volume, MI.-only, and our hybrid model on 2D forced turbulence task. Our hybrid models provide some improvement on the accuracy of spectral-only methods, which themselves perform remarkably well compared to recently proposed ML methods. Furthermore, our results show that a key modeling choice for both hybrid and pure ML. models is the use of velocity- rather than vorticity-based state reprosentations.
1.1 Related work
The study of turbulent fluid dynamics vast. We refer to Pope (2000) for a thorough introduction. Classical approaches lend to derive mathematical approximations to the governing equations in an a priori manner. Recently, there has been an explosion of work in data-driven methods the interface of computational fluid dynamics (CFD) and machine learning (ML). We loosely organize this recent work into three main categories:
Purely Learned Surrogates fully replace numerical schemes from CFD with purely learned surrogate models. A number of different architectures have been explored, including multi-scale convolutional neural networks (Ronneberger et al., 2015; Wang et al., 2020), graph neural networks (Sanchez SSanchez-Gonzalez et al., 2020), and Encoder-Process- Decoder architectures (Stachenfeld et al., 2022).
Operator Learning seeks to learn the differential operator directly by mimicking the analytical properties of its class, such as pseudo-differentall or Fourier integral operators but without explicit physics-informet components. These methods often leverage the Fourier transform (Li et al., 2021; Tran et al., 2021) and the off-diagonal low-rank structure of the associated Green's function (Fan et al., 2019; Li et al., 2020).
Hybrid Physics-ML is an emerging set of approaches which aim to combine classical numerical methods with contemporary data-driven deep lcarning techniques. These approaches use high-resolution simulation data to learn corrections low-resolution numerical schemes with the goal of combining the best of both worlds - the simplicity of PDE-based governing equations and the expressive power of neural networks.
2

---

Published in Transactions on Machine Learning Research (04/2023)
Our work falls into this category, along with growing body of work (Mishra, 2018; Bar-Sinai et al, 2019; Kochkov et al., 2021; Bruno et al., 2021; List et al., 2022; Frezat et al., 2022; Huang et al., 2022).
Some of these works stand out as being closely related to ours but with important differences. Huang et al. (2022) also develops a hybrid approach but focuses exchusively on ODEs, not PDEs, and it is not clear whether their techniques for stable simulations with large time-step size would also apply PDFs Tn their recent work, Frezat et al. (2022) augmented spectral solvers for closure models exclusively focused on climate modeling. We also learn corrections to spectral methods but rather than focusing on climate models, tackle more general PDEs, provide key ML architecture choices, and include comparisons to classical numerical methods on multiple grid resolutions.
Unlike our work, San & Maulik (2018) and Subel et al. (2021) take fundamentally different approaches to ours. San & Maulik (2018) improves classical approaches for closure modeling by using a neural network with a single hidden layer to learn fino-tuned local corrections instead of the usual global correction methods. Their work requires specialized knowledge of classical methods (POD, closure models, etc.) which is orthogonal to our larger goal of automating these approaches using data-driven techniques. Subel et al. (2021) focus a data augmentation scheme to overcome stationary shocks that arise in Burgers' equation. Their proposed shifting mechanism could be used in conjunction with our work, but represents a different research direction.
2 Spectral methods for fluids
Spectral methods are a powerful method for finding high-accuracy solutions to PDEs and are often the method of choice for solving smooth PDEs with simple boundary conditions and geometries. There is an extensive literature the theoretical and practical underpinnings of spectral methods (Gottlieb & Orszag, 1977; Gottlieb & Hesthaven, 2001; Canuto et al., 2007; Kopriva, 2009), particularly for methods based on Fourier spectral collocation (TTrefethen, 2000; Boyd, 2001). Below, we provide a succinct introduction.
2.1 Partial differential equations for turbulent fluids
Let u Rd x + + Rdd be time-varying vector field, for dimensions d and d'. We study PDEs of the form
0 = Du + N(u).
(1)
plus initial and boundary conditions. Here D is lincar partial differential operator, and N is non-lincar term. Equations of this form dictate the temporal evolution of u driven by its variation in space. In practice, PDEs are solved by discretizing in space and time, which converts the continuous PDE into a set of update rules for vectors of cocfficients to approximate the state u in some discrote basis, c.g., on a grid. For time-dependent PDEs, temporal resolution must be scaled proportionally to spatial resolution to maintain an accurate solution. Thus, runtime for PDE solvers is O(n+1),, where is the number of spatial dimensions and n is the number of discretization points per dimension.
For simulations of fluids, the differential operator is typically either diffusive, D= 8. or hyper-diffusive D=
8z. D fluid velocity to become uniform due to internal friction, and convection is the tendency of the fluid to be transported by its own inertia. Turbulent flows are characterized by a convective term that is much stronger than diffusion.
And, the non-linearity is convective term, N= (u²)) = u Diffusion is the tendency of the
For most turbulent flows of interest, closely approximating the exact PDE (known as Direct Numerical Simulation) is computationally intractable because it requires prohibitively high grid resolution. Instead, coarse-grained approximations of the PDE are solved, known as "Large Eddy Simulation" (LES). LES aug- ments Equation (1) by adding a correction term determined by a closure model to account for averaged effects over fine spatial length scales. In practice, this term is often omitted due to the difficulty of deriving appropriate closure formulas ("implicit LES") and the PDE is simply simulated with the finest computa- tionally fcasible grid resolution. In our case, correction terms are an opportunity for machine learning. If we can accurately model PDEs on coarser grids with suitable correction terms, we may be able to significantly reduce the computational cost of large-scale simulations.
3

---

Published in Transactions on Machine Learning Research (04/2023)
32x32
64x64
128x128 28
1024x1024
Average Time Until Currelation < 0.95
10-
(Converged)
FVM
Spectral
5+
4-
FVM
ML. Spl. Op.
EPD
128
256
5i2
1024
Resolution
Figure 1: Comparing finite volume (FVM) to spectral convergence for 2D turbulence. Left: Vorticity fields evolved up to time 9.0 with different grid sizes and numerical methods, starting from an identical initial condition. Qualitatively, it is clear that at resolution 64x64, the spectral method has already captured most of features of the high-resolution state. The FVM looks sharp, but clearly differs. At sufficiently high resolution, the methods converge to the same solution. Right: Here we compare each method with high- resolution 2048x2048 ground truth. We plot the time until the first dip below 0.95 of the correlation with the ground truth. Note that the initial conditions for the finite volume and the spectral method are sampled from the same distribution, but are not identical. At resolution 256, the spectral method has converged within the accuracy limitations of single precision, so we omit higher resolutions. For the EPD baseline and our ML Split Operator (ML Spl. Op.) models, we show performance across five different neural network parameter initializations.
2.2 The appeal of the Fourier basis for modeling PDEs
Let us further assume that u(x, t) [0, 2m] x R + IR in Equation (1) is 2x-periodic, square-integrable for all times i, and for simplicity, one-dimensional. Consider the Fourier coefficients û' of u(r, t), truncated to lowest K + frequencies (for even K):
2n u(3r, -
n' (ύ' K0'
(k/a)) where ii
ik- dr.
(2)
2T
n' a vector representing the solution at time 1, containing the coefficients of the Fourier basis e for k [-K/2,..... k/2). In general, the integral üt has no analytical solution and so we approximate it using a trapezoidal quadrature on K +11 points. Namely, we approximate ae it by sampling u(x, t) on an equispaced grid. The Fast Fourier Transform (FFT) (Cooley & Tukey, 1965) computes these Fourier coefficients efficiently in log-linear time. Spectral methods for PDFs leverage the fact that differentiation in the Fourier domain can be calculated by element-wise multiplication according to the identity 0,û = ikūg. This, in turn, makes inverting linear differential operators easy since it is simply dement-wise division (Trefethen, 2000).
When time-dependent PDEs include non-linear terms, spectral methods evaluate these terms on grid real-space, which requires forward and inverse FFTs inside cach timc-step. These transforms introduce two sources of error. First, the quadrature rules for पi produce only an approximation to the integral. Second, there will be truncation error when the number of frequencies K + is less than the bandwidth of u(z, t). Remarkably, it is well-known fact of Fourier analysis that both of these errors decay super-algebraically (i.e., faster than any power of 1/K) if u(x, !) is periodic and is in O (Trefethen, 2000).* 1 Thus, relatively few discretization points are needed to represent solutions which are C
Bocause of these favorable convergence properties, spectral methods often outshine their finite difference counterparts. For example, spectral methods are used for large-scale simulations of turbulence on GPU super-computers (Yeung & Ravikumar, 2020). See Figure 1 for simple comparison in the case of 2D turbulence. This motivates to start from spectral methods, rather than finite difference methods, and further improve them using machine learning.
Convergencee exponential for analytic (Tadmor, 1986).
4

---

Published in Transactions on Machine Learning Research (04/2023)
3 Learned split operators for correcting spectral methods
Similarly to classical spectral numerical methods, we solve Equation (1) by representing state in finito Fourier basis as ú' and integrating forward in time. We model the fully known differential operators D and N using a standard spectral method, denoted Physics-Step. The machine learning component, denoted Caarned-Corrections -;),, contains tume-able parameters 0 which are optimized during training to match high-resolution simulations.
Building upon traditional physics- based solvers, we use a simple explicit-Euler time integrator for the cor- rection term (See App. D for brief explanation of explicit methods.) This yields the following update equation:
a Physics-Step(ii) - h- Bhaaned-Corrections 0),
where h E IR is the time step size. Ultimately, our work is aimed at two and three dimensional problems so rather than present an abstract schematic as presented in the equation above (Eq. (3)), we chose to present a more detailed, realistic schema of 2D turbulence (Eq. (6)). The spatial state variable thus not. the generalized as in Equation (3) but rather, the vorticity w as defined in Section 4.5
3.1 Convolutional layers for encoding the physical prior of locality
For small time step h, the solution at time + h only depends locally in space the solution at time (within a dependency cone usually encoded by the Courant-Friedrichs-aare previous work (see Sec. 1.1, Hybrid Physics-ML) we incorporate this assumption into the machine learning component of our model using Convolutional Neural Networks (CNN or ConyNet) which, by design, only learn local features. We now provide a high-level view our modeling choices:
(CFL) condition). Following
Real-space vs. frequency-space. Since the Fourier basis is global, each coefficient úf contains informa- tion from the full spatial domain as shown in Equation (2). Thus, in order to maintain spatial features, which are local, we apply the ConvNet component of our Neural Split Operator model in real-space. This is accomplished efficiently via inverse-FFT and FFT map the signal back and forth between frequency- and real-space.
ConvNet Padding. In this work, we tackle problems with periodic boundary conditions. This makes periodic padding natural choice for the ConvNets.
Neural architecture. For both 1D and 2D problems we used an Encoder-Prooess--eer Decoder (EPD) archi- tecture (Stachenfeld et al., 2022) facilitate fairer comparisons pure neural network baselines. While various forms of optimization is possible, such model parallelism among others, we consider this to be outside the scope of this exploratory work. See Section 4.2 for detailed description of EPD models and Appendix C for further information.
3.2
The split operator method for combining time scales
Due to a variety of considerations - numerical stability, computational feasibility, etc. - each term of a PDE often warrants its own time-advancing method. This motivates split operalor methods (Strang, 1968; McLachlan & Quispel, 2002), popular tool for solving PDEs which combine different time integrators. In this work, we use the split operator approach incorporate the additional terms given by the neural network.
In the usual pattern of spectral methods, Physics-Step itself is split into two components corresponding to the D and N terms of Equation (1). The term of Physics-Step is solved with a Crank-Nicolson method and the N part solved using explicit 4th-order Kange-Kuttar which effectively runs at a time-step of h/4. The
component is solved using vanilla first-order Euler timo-step with step size h, which when compared to incorporating the learned component in the 4th-order Runge-Kutta solver, is more accurate and has 4x faster runtime (see Sec. 4). Alternatively, omitting Physics-Stee from Equation (3) gives
5

---

Published in Transactions on Machine Learning Research (04/2023)
Learned-Correction
Velocity-Solve
IFFT
(v vty 201
NN
FFT
-h
WHt1
+
Physics-Step
L+ Physics-Step
Diagram summarizing our model described in Equation (3) for the 2D Navicr-Stokes equation. The input, vorticity W², is processed by two independent components Learned-Correction and Physics-Slep operating at different time-scales. The output Learned- Correction is weighted by h, as in a basic first-order Euler stepping scheme, and combined with the output of the Physics-Step to give the state the next time step.
Figure 2:
a Neural ODE model (Chen et al., 2018) with first-order time-stepping which is precisely the EPD model described by Stachenfeld et al. (2022). This model serves as a strong bascline shown in Section 4.
3.3 Physics-based solvers for calculating neural-net inputs
A final important consideration is the choice imput representation for the machine learning model. We found that ML models operate better in velocity-space whereas vorticity is the more suitable represen- tation for the numerical solver. We were able to improve accuracy by incorporating data pre-processing step, e.g.; velocity-solve operation, for the inputs of the neural network compo- nent (See App. D for brief explanation of velocity-solve). Our overall model is Cearned-Correetion (F (tar://tttttttafoorooffoo Transform("'))) where NN is implemented a periodic ConvNet. For our 1D test prob- lems, Mtate-Transonne is the identity transformation, but for 2D Navier-Stokes (depicted in Fig. 2) we perform a velocity-solve operation to calculate velocity from vorticity. This turns out to be key modeling choice, as described in Section
physics-based
3.4 Training and evaluation
Data preparation for spectral solvers. For ground truth data, we use fully resolved simulations, which we then downsample to the coarse larget resolution. Choosing the downsampling procedure is key decision for coarse-gre coarse-grained solvers (Frezat et al., 2022). In this work, the fully resolved trajectories are first truncated to the target wavenumber (i.c., an ideal low-pass filter). Then, we apply exponential filter of the form ûg exp(-a]/t akk////aamax htt//max/"" (2) U where üg denotes the filtered field and & is the k-th wavenumber (Canuto al., 2006). We obtained stable trajectories using relatively weak filter with = 6 and 16. The exponential filter smooths discontinuities in the PDE solution, which otherwise manifest themselves globally in real space as ringing artifacts known as Gibbs Phenomena (Canuto et al., 2006).
Filtering also used to correct aliasing errors in spectral methods that arise when evaluating non-lincar terms (Gottlieb & Hesthaven, 2001). Ideally, filtering for aliasing errors would be spatially adaptive (Boyd, 1996). In practice, however, filters for both aliasing and truncation often chosen heuristically. While one might aspire to learn these heuristics from data, attempts at doing this is because insufficiently filtered simulations will often entircly diverge rather than accumulate small errors. Thus, in addition to filtering the downsampled training data, we also used the exponential filter on the outputs of Physics-Step of Equation (3) and consider this to be another component of Physics- Step. Omitting this filtering also resulted global errors which were impossible for the Icarned component correct. Figure summarizes our data generation pipeline, including an explicit filtering step. In Section we present failure mode model without filtering (see Fig. 5).
unsuccessful In part
6

---

Published in Transactions on Machine Learning Research (04/2023)
Initial State
High Res Sim
Downsampled
Filtered
Training Example
Downsampled Initial State
Prediction
Solver
Model Params
Figure 3: Diagram of our training pipeline. Starting with high resolution initial state, we run it forward using a high-resolution spectral solver. Then we downsample by truncating higher frequencies in the Fourier representation. This can cause "ringing effects" for which the standard approach is to apply filtering operator. Finally, we train a solver to mimic this process as closely as possible, as measured in f--loss.
We experimented with larger time-step sizes, c.g. 0.1, 1.0 but found the resulting models to be unstable. This is because the physics-componett of the model quickly becomes unstable at large step sizes which propagates to the overall model.
Training loss. We train models to minimize the squared-error over some number of unrolled prediction steps, which allows our model to account for compounding errors (Um et al., 2020). Let u(z, t) denote our prediction at time t, then our training objective is given by B Ex, u t) u2x t) The scaling constant 3 is chosen so that the loss of predicting "no change" is one, i.e., 3 Em.er |u(z, 0) u(x, t) Relative to the thousands of time-steps over which we hope to simulate accurately, unroll over relatively small number during training (e.g., T = 32 for 2D turbulence) because training over long time windows is less efficient and less stable (Kochkov et al., 2021).
Measuring convergence. We seek to optimize the accuracy of coarso-grained simulations. More specif- ically, we aim to make a coarse resolution simulation as similar possible to high-resolution simulation which has been coarse-grained in post-processing. Following Kochkov et al. (2021), we measure conver- gence to a fully resolved, high-resolution ground-truth at each time in terms of mean absolute error (MAE), correlation, and time until correlation less than 0.95. MAE defined as Eg [u(a,) t) u(xx t)l and correlation is defined Corrfu(-, (), ()) E, (T(,,)) () u(a, ())/([[((2,, ()ll2lll((a l)22 (since flows have mean zero). Finally, we compute the first time step in which the correlation dips below 0.95, i.e. u(,,,,
min(t Corrlu(-,t), u(- < 0.95). While MAE is precise up to floating point precision, it is not readily in- terpretable. Whereas, correlation is less sensitive but more interpretable. For this reason, we prefer to report correlation. This potential redundancy is demonstrated in our experiment on the KS equation (Sec. 4.3), where we included both MAE and correlation.
For the 2D Navier-Stokes equation, we only report correlation because it is sufficient. For the unstable Burgers' equation, the situation was the opposite: We only reported MAE because measuring correlation did not provide additional insight - all models had correlation values close to 1.0, whereas MAE sensitive to the improved performance of our method.
4 Results
4.1 Model equations
We showcase our method using three model equations which capture many of the algorithmic difficulties present in more complex systems: 1D Kuramoto- -Siassii--ivassiskky (KS) equation, 1D unstable Burgers' equation, and 2D Kolmogorov flow, variant of incompressible Navier-Stokes with constant. forcing.
7

---

Published in Transactions on Machine Learning Research (04/2023)
The KS equation has smooth solutions which spectral methods are well-designed to solve. Therefore, it is not surprising that, while our method does improve over spectral-only methods, that improvement is not significant. On the other hand, the unstable Burgers' equation presents a test case in which classical spec- tral methods struggle near discontimuities. Here, our method outperforms spectral-only methods. Finally, with two-dimennional Kolmogorov flow, we demonstrate our method on a more challenging fluid simulation. Kolmogorov flow exhibits multiscale bohavior addition to smooth, chaotic dynamics. Without any modifi- cation, spectral-only methods are already competitive with the latest hybrid methods (Kochkov et al., 2021). Similar to the 1D KS equation, our method is able to provide some improvement in this case.
4.2 Baselines
Our hybrid spectral-MT. method naturally gives rise to two types of baselines: spectral-only and MT-only. On all three equations, compare to spectral-only at various resolutions, e.g., "Spectral 32" refers to the spectral-only method with 32 grid points.
For Kolmogorov flow, we compare to two additional baselines: (1) the Prccode-PProces--
(EPD)
model (Stachenfeld et al., 2022), a statc-of-the-art ML-only model, and (2) a finite volume method (FVM) as implemented by Kochkov et al. (2021) a standard numerical technique which serves as an alternative the spectral method.
The EPD architecture acts in three steps. First, the spatial state is embedded into a high-dimensional space using a leed-forward neural network encoder architecture. Then, either another feed-forward or recurrent architecture is applied to process the embedded state. Finally, the output of the process step is projected to back to the spatial state using another feed-forward decoder module. This method has recently achieved state-of-the-arr on a wide range of one,, two-, and diree-dimeiioonl problems. For consistency, we used an identical EPD model for the learned component of our model (Eq. (3)).
Finally, we included second-order FVM on staggered grid to illustrate the strength of the spectral baselines. This gives context to compare to previous work in hybrid methods (Bar-Sinai et al., 2019; Kochkov et al., 2021) which use this FVM model as baseline.
We avoided extensive comparisons lo classical subgrid modeling, such as Large Eddy Simulation (LES), since LES is itself a muanced class of methods with many tunable parameters. Instead, our physics-only baselines are implicit LES models - a widely-used and well-understood model class which serves as a consistent, parameter-fi bascline. In contrast, explicit sub-grid-scale models such as Smagorinsky models include tunable parameters, with the optimal choice depending on the scenario of interest. Furthermore, explicit LES models typically focus matching the energy spectrum rather than minimizing point-wise errors (List et al., 2022).
The Prramoto-SSvashinske (KS) equation
4.3
(KS) equation is model of unstable flame fronts. In the PDE literature, it is popular model system because it exhibits chaotic dynamics in only a single dimension. The KS equation is defined as
0u= = ---u-- U.u.
(4)
The three terms on the right hand side of this equation correspond to convection, hyper-diffusion and anti- diffusion, which drives the system away from equilibrium.
Because solutions are smooth, spectral methods are able to capture the dynamics of the KS equation ex- tremely well. At a resolution of 64, the simulation is already effectively converged to the limits of single precision arithmetic, ie., it enjoys perfect correlation with a high-resolution ground truth.
Considering how well-suited spectral methods are for modcling this equation, remarkable that our ML- Physics hybrid model able to achieve any improvement at all. Looking qualitatively the left panel of Figure 4, there is clear improvement over the spectral 32 baseline.
8

---

Published in Transactions on Machine Learning Research (04/2023)
Spectral 32
Spectaal 1924
MI Spitt Diwrbbo
Mean Absplute
Currilation
0.5
Spectral
0.0
10
Spectral Spt Op.
-20
40
40
20
40
10
20
30
40
50
66
10
20
30
50
Simulotion
Simcdotion Time
Tinoo
Time
Figure 4: Comparing our model to spectral-only basclines on the KS equation. The spectral method is essentially converged to a 1024 ground truth model at a resolution of 64. Vertical dashed line indicates the first time step in which our model's correlation with the ground truth less than 0.9. Our model still able to some improvement over coarse resolution (Spectral 32) bascline.
Spectral
Spectrl
Solft
Mean Abuuhate Breer
10r
81
10a4
60-
Spectral 54
Spectral
10a4
Spectra 256
Spectnal 130
ML. Spt. 0p. 32
-82-
720
1400
150
1600
1000
Semulation Time
Simalation Tims
Somalattion Time
Simulation Time
Simulation Time
Figure 5: Comparing Vertical dashed line indicates the first time step in which model's correlation with the ground truth is less than 0.9. Not only is our model stable unlike the coarse (Spectral 32) bascline, it also performs on par with a medium-grain (128) baseline, demonstrating a 4x improvement in spatial accuracy. If the Physics-Step component Equation (3) does not include the exponential filter as described in Section 3.4, the Correcttooooeti component diverges completely (Left panel, "No Filter").
model to spectral-only baselines the unstable Burgers' equation (Eq. (5)).
4.4 Unstable Burgers' equation
Burgers' equation is simple one-dinensional non-linear PDE which used toy model of compressible fluid dynamics. The characteristic behavior of Burgers' equation is that it develops shock waves. Practition- ers use this equation to test the accuracy of discretization schemes near discontinuities.
Like turbulent flows, the behavior of Burgers' equation dominated by the convection, but unlike turbulent flows, not chaotic. Prior studies ML applied to Burgers' equation imposed random forcings (Bar-Sinai et al., 2019; Um al., 2020), but due the non-choatic nature of the equation, discretization errors decay rather than compounding over time. Instead, we use an unstable viscous Burgers' equation which slightly amplifies low frequencies in order to make the dynamics chaotic. Sakaguchi (1999) provides the following definition:
( g(x - x' L'))
0,u = uu. u+ + voLu -
(5)
for viscosity v > and domain size L. convection, diffusion and a scale-selective amplification of low-frequency signals. The convolution g(æ - x²) amplifics low froquencics, docaying smoothly zero as the frequency increases but best described in the Fourier space, which defer to Appendix A.
The three terms in this unstable Burgers' equation correspond to
Shock waves challenging to model using the Fourier basis because there are discontinuities in the solution (Sec. 3.4). Due to the mixing of length scales introduced by convection, these errors can quickly propagate to affect the overall dynamics. This often results in instability at low resolutions, evidenced in our results (Fig. 5, resolution 32). While there are classical methods for addressing discontinuities, e.g., adding a hyper- viscosity term, they often modify the underlying dynamics that they are trying to model. The inadequacy of unaltered spectral methods for solving this problem explains why our method is able to achieve approximately 4x improvement in spatial resolution spectral baseline. Results summarized in Figure 5.

---

Published in Transactions on Machine Learning Research (04/2023)
Time Ston
000
TIme Step
900
Time Stop 1000 Time Stop
1100
Simelation
Speetml 204022000
7.6-
Mis opt op w urp..
5.0-
0.9
NR at-tpp op uooo
1P0
-6.0
-26-
ML Openite Spinin
0.2
Spectral
Soectrol
Spectrol 64
Mr. Sab fo
upp from
-7.5
arnes
-10.0
000 1000
1000
2000
2300
2.0
3.0
5.0 7.0
Figure 6: Benchmarking our Neural Split Operator (ML Spl. Op.) to spectral-only baselines and state-of-the art MT-only, EPD model (Stachenfeld et al., 2022). Left: Visualizing model states starting at the time when divergence begins to occur. At time step = 800, EPD, Neural Split Operator model, and the high-resolution spectral-only 1024x1024 correlate well. By time step = 900, both learned models have begun to decorrelate with the ground truth albeit in different ways. Right: Comparing model variants across different random initializations. On the y-axis we measure the time to the first time step with correlation<O..55 with the high-resolution spectral ground truth. Then compared our model with the EPD model across three data representation variations: velocity only (v), vorticity only (w) and velocity-vorticity concatenation ((v,w) well as with Nonlinear Term Correction model (N.L. Crrtn.). Shown here are the best performing models over five different random neural network parameter initializations.
4.5
2D turbulence
We consider 2D turbulence described by the incompressible Navier-Stokes equation with Kolmogorov forc- ing (Kochkov et al., 2021). This equation can be written either in terms of velocity vector field v(r, v) = (V., Vy) or a scalar vorticity field w == V2 0,V. (Boffetta & Ecke, 2012). Here we use a vorticity formulation, which is most convenient for spectral methods and avoids the need to separately enforce the incompressibility condition vvv= 0. The equation is given by
d.w - - -V - aw+
(6)
where the terms correspond to convection, diffusion, linear damping and a constant forcing / (App. B). The velocity vector field can be recovered from the vorticity field by solving a Poisson equation - 2 w for the stream function & and then using the relation v(z, y) - -0, This velocity-solve operation is computed via element-wise multiplication and division in the Fourier basis.
For this example, we compared four types of models: spectral-only, FVM-only, EPD (MI-onnly, and two types of hybrid methods: our split operator method and a nonlinear term correction method. The ML models are all trained to minimize the error of predicted velocities, which are obtained via a velocity-solve for the spectral-only, EPD and hybrid models that use vorticity as their state representation. The nonlinear term correction mothod uses a neural network to correct the inputs to the nonlincar term v Vw and avoids using classical correction techniques, e.g., the so-called three-over-two rule (Orszag, 1971). Specifically, this approach takes inputs (v, Vw) and outputs scalar field which is incorporating into the explicit part of the physics solver, e.g. dي = ((v Ve + maalincar-Correcttions w; 0)) + 2 aw
Results are summarized in Figure 6. Comparing Spectral 128 to FVM 1024 we can see that the spectral- only method, without any machine learning, already comes close to a similar improvement. Thus, a -2x improvement over the spectral-only baseline achieved by both the EPD and our hybrid model is comparable to the -8x improvement over FVM achieved in Kochkov al. (2021) in which they used a hybrid MT.-Physics model.
10

---

Published in Transactions on Machine Learning Research (04/2023)
State representation. We compare three representations for inputs to neural networks in our EPD and split-operator methods: velocity, vorticity, and velocity-vorticiyy concatenation. All models represent internal state with vorticity and learn corrections to vorticity, with the exception of velocity-only EPD models. These models never compute vorticity - they use velocity to represent state as well as Iearned corrections, thus malching previous work by Kochkov et al. (2021) and Stachenfeld et al. (2022). Interestingly, this velocity- only representation has the best performance across the board, particularly for the fully Iearned EPD model, even outperforming velocity-vorticity concatenation (Fig. 6, right panel).
The nonlinear term Vw in Equation (6) makes of both vorticity and velocity representations. This indicates that in order to model the dynamics, the network must learn to solve for velocity - global operation making it challenging for a ConvNet restricted to local convolutions to learn. In contrast, computing vorticity from velocity only requires evaluating derivatives, which can be easily evaluated with local convolutions, e.g., via finite differences. This may explain the relatively-worse performance of Fourier Neural Operators using the velocity-basee based representation of Stachenfeld et al. (2022) versus the vorticity- based representation of Li el. al. (2021).
Full versus nonlinear-term-onl correction. To understand the value of time-splitting for the learned correction, we also performed an experiment where we instead incorporated our MI. model into the nonlinear part of Physics-Step (Eq. (3)). This entails solving the convection term with the same 4th order explicit Runge-Kutta method used in the numerical solver. Since velocity-space performed best with the split operator model, we also used it here. The non-split learned correction has significantly worse performance. Anecdotally, we observed that training with high order Runge-Kutta Kutta methods was significantly less stable. Models were much less robust to small changes in learning rates and dilation rates (App. C) than models trained with first-order Runge-Kutta, We believe this may be due to the increased difficulty of propagating gradients through a network which is effectively four times deeper. Another possibility is that the polynomial approximation of the high-order Runge-Kutee Kutta method introduces larger errors at the beginning of training when the model is less accurale.
Analyzing learning curves. To further compare the merits of different ML approaches, we compare the learning curves for select models in Figure 7. First, we notice that although our best EPD and hybrid models are similarly accurate once trained, the hybrid models require drastically less compute to achieve fixed evaluation metrics. For example, squared error of 2000 at 1024 time-steps requires only 3600 training steps for the median ML split operator velocity model, versus 164 600 training steps for the EPD velocity model, corresponding to 4.4 versus 170 TPU-core hours.
These learning curves also reveal that our validation loss, in this case caleulated over 32 unrolled time-steps, is not necessarily indicative of validation performance over much longer unrolls, Although our models have never seen the exact validation data during training, many of them are still able to "overfit" to the task of predicting short. trajectories, at the cost of generalization to long unrolls. We also observe that different ML architectures overfit to different extents. In particular, pure ML models and models with access to vorticity overfit more than our best hybrid model, the ML split operator with only velocity inputs. The ML split operator with access to both velocity and vorticity is particularly interesting, because it seems to undergo a phase transition at around 30 training steps, with correlation for long unrolls dropping dramatically as the model shifts into the "memorization" regime.
5 Discussion
In this paper, we demonstrated the potential of MI-augmented spectral solvers to improve upon the accuracy of spectral-only methods. We also identified several key physically motivated modeling choices velocity- representation, first-order time stepping to improve sensitivity to hyperparameters, and the removal of global spatial artifacts which improve training for both ML-only and hybrid models. Pure ML models can match the accuracy of hybrid models, but are considerably more expensive to train and show more indications of overfitting.
11

---

Published in Transactions on Machine Learning Research (04/2023)
Training loss on validation dota
100
Squared error at time step 1024
10-1
10
ML Spl. Op. V
10-2
EPD: V
Ope
EPD: (v.
10.2
10.
10P
Training steps
Training stops
Figure 7: Learning curves for selected models, showing evaluation metrics after different numbers of training steps on the held-out validation dataset. The training loss is evaluated over the first 32 unrolled steps, whereas squared error is evaluated at 1024 steps (t = 7.2). The latter better indicates the model performance we care about for predicting long trajectories. Separate lines show learning curves for each of our five randomly initialized training runs.
Traditional spectral methods are a powerful set of approaches for solving equations with smooth, periodic solutions. is yet unclear whether ML-based solvers can achieve meaningful computational speed-ups over classical spectral methods on PDEs of this class. In contrast. to prior work which showed computational speed-ups of up 1-2 orders of magnitude over baseline finite volume (Kochkov et al., 2021) and finite difference methods (List et al., 2022), the roughly 2x decrease in grid resolution for 2D turbulence with the MT. split operator would allow for at most 8x reduction in computational cost. However, our neural network for learned corrections is about 10x slower than Physics-Step, which counteracts this potential gain. Small speed-ups might be obtained by using smaller networks or applying corrections less frequently, but overall there is little potential for accelerating smooth, periodic 2D turbulence beyond traditional spectral solvers.
Hybrid ML-spoctral methods may enjoy more significant improvements on other problems. Examples might include PDEs with less smooth solutions, such as 3D turbulence, where energy decays k-5/3 versus k-3 in 2D (Pope, 2000). Or, global atmospheric models, where spectral methods do not achieve exponential convergence (Williamson, 2008). Cases where the exact governing equations are partially unknown, such as physical parameterizations for climate models (Brenowitz & Bretherton, 2018), also present an opportunity for combining physical simulation with machine learning components.
12

---

Published in Transactions on Machine Learning Research (04/2023)
References
Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P. Brenner. Learning data-driven discrotizations for partial differential equations. Proceedings of the National Academy of Sciences, 116, July 2019. ISSN 0027-8424, 1091-6490. doi: 10.1073/ 21173//pas..1//// 1814488116 3, 8, 9, 17)
Guido Boffetta and Robert E. Ecke. Two-Dimensional Turbulence. Annual Review of Fluid Mechanics, 44, January 2012. ISSN 0066-4189, 1545-4479. doi: 10.1140/anmurev-- 1201111007000000114 (A 10)
John P. Boyd. The Erfc-Log Filter and the Asymptotics of the Euler and Vandeven Sequence Accelorations.
In Min, A V and Scott, Ridway (ed.), Proceedings of the Third International Conference on Speciral and High Order Methods, Houston, 1996. Houston Journal of Mathematics. (^
John P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Books on Mathematics, Dover Publications, Mineola, NY, second edition, 2001. ISBN 0486411834 9780456411835, (*
N. D. Brenowitz and C. S. Bretherton. Prognostic Validation of Neural Network Unified Physics Parame- terization. Geophysical Research Letters, 45, 2018. doi: https://doi.org//
(* 12)
Oscar P. Bruno, Jan S. Hesthaven, and Daniel V. Leibovici. FC-based shock-dynamics solver with neural- network localized artificial-viscosity assignment, 2021. (^
Keaton J Burns, GeolTrey M. Vasil, Jeffrey S Oishi, Daniel Lecoanet, and Benjamin P Brown. Dedalus: A flexible framework for numerical simulations with spectral methods. Phys. Rev. Research, 2, April 2020. 2)
C. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods: Evolution to Compler Ge- ometries and Applications Fluid Dynamics. Scientific Computation. Springer Berlin Heidelberg, 2007. ISBN 9783540307280.. (* 3, 16)
Claudio Canuto, M Yousuff Hussaini, Alfio Quarteroni, and Thomas A Zang. Spectral methods: Fundamentals in single domains. Scientific computation. Springer, Berlin, Germany, 1 edition, April 2006. ( 6)
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural Ordinary Differential Equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Newral Information Processing Systems, volume 31. Curran Associates, Inc., 2018. (^
J. W. Coolcy and J. W. Tukey, An Algorithm for the Machine Calculation of Complex Fourier Scries. Math. Compul., 19, 1965. ISSN 00255718, 10886842. (
Paul Durbin. Some Recent Developments Turbulence Closure Modeling. Annu. Rev. Fluid Mech., 50, 2018. 2*
Yuwei Fan, Jordi Feliu-Fabà, Lin Lin, Lexing Ying, and Leonardo Repeda-Núñez. A multiscale neural network based on hierarchical nested bases. Research in the Mathematical Sciences, 6, Mar. 2019. ISSN 2197-9847. doi: 10. 10.1007/ 000777//00687-01-00155
2)
Hugo Frezat, Julien Le Sommer, Ronan Fablet, Guillaume Balarac, and Redouane Lguensat. A posteriori learning for quasi-geostrophic turbulence parametrization. Xrv, April 2022. 2, 3)
D. Gottlieb and J. S. Hesthaven. Spectral methods for hyperbolic problems. Journal of Compulational and Applied Mathematics, 128, March 2001. ISSN 0877-0427. doi: 10.1066/10377--017000)011111 (* 3, 6)
David Gottlieb and Steven A. Orszag. Numerical Analysis of Spectral Methods. Society for Industrial and Applied Mathematics, 1977. doi: 1011137/1.7811611904455 (*3)
Zhongzhan Huang, Senwei Liang, Hong Zhang, Haizhao Yang, and Liang Lin. Accelerating Numerical Solvers for Large-Scale Simulation of Dynamical System via NeurVec, August 2022. S1Ti:2208.006880 (^3)
13

---

Published in Transactions on Machine Learning Research (04/2023)
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine learning learning-accelerated accelerated computational fluid dynamics. Proc. Nall. Acad. Sei. U. S. A., 118, May 2021. ( 2, 3, 7, 8, 10, 11, 12, 17)
David A. Kopriva. Implementing Spectrol Methods for Partial Differential Equations: Algorithms for Sci- entists and Engineers. Springer Publishing Company, Incorporated, 1st odition, 2009. ISBN 9048122600.
3)
Zongyi Li, Nikola Kovachki, Kamyar Azizzadeneshelii Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Multipole Graph Neural Operator for Parametric Partial Differential Equations. In Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS'20, 2020. 2*
Zongyi Li, Nikola Kovachki, Kamyar Avizadeneshell, Burigede Liu, Kaushik Bhattacharya, Andrew Stu- art, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations. Xrxiv:2010. 0.8800888550008 fcs, math], May 2021. arXiv: 2010.08895. 2* 11)
Björn List, Li-Wei Chen, and Nils Thuerey. Learned Turbulence Modelling with Differentiable Fluid Solvers. arXiv:2202.0600 [physicss, February 2022. arXiv: 2202.06988. ( 2, 3, 8, 12)
Robert I. McLachlan and G. Reinout W. Quispel. Splitting methods. Acto Numerica, 11, 2002. doi: 10.0017/0092299222000055 5*
Siddhartha Mishra. A machine learning framework for data driven acceleration of computations of differential equations. Mathematics in Engineering, 1, 2018. ISSN 2640-3501. doi: 10.3934/Minee 2018.1.118. M.1.../ (*
Steven A. Orszag. On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-
Wavenumber Components. Journal of the Almospheric Sciences, 28, September 1971. ISSN 0022-4928, 1520-0469. doi: 10.1175/153 0eee//152000-4999 (489(991)088-107700111111
(A 10)
Stephen B. Pope.
Turbulent Flows. Cambridge University Press, August 2000. ISBN 9780521598566. (^ 2,
12)
Stephen Pope. Ten questions concerning the large-eddy simulation of turbulent flows New Journal of Physics, 6, March 2004. ISSN 1367-2630. doi: 11188/33772230/6/// 1/0/035. 2*
Christopher D Roberts, Retish Senan, Franco Molteni, Souhail Boussetta, Michael Mayer, and Sarah PE Keeley. Climate model configurations of the ECMWF Integrated Forecasting System (ECMWF-IFS cycle 43r1) for HighResMIP. Geoscientific model development, 11, 2018. 2)
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional Networks for Biomedical Image Segmentation. In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi (eds.), Medical Image Computing and Computer-Assisted Intervention MICCAI 2015, Cham, 2015. Springer International Publishing. ISBN 978-3-319-245444. 2*
Hidotsugu Sakaguchi. Chaotic dynamics of an unstable Burgers equation. Physica D: Nonlinear Phenomena, 129, May 1999. ISSN 0167-2789. doi: 10.1066/50000--775000))00111
(* 9,9, 16)
Omer San and Romit Maulik. Extrome learning machine for reduced order modeling of turbulent geophys- ical flows. Physical Review E, 974):044,3 April 2018. doi: 20.1100//hyyaRee...7..0..
Publisher:
American Physical Society. (*
Alvaro Sanchez- Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying. Jure Leskovee, and Peter W. Battaglia. Tearning to Simulate Complex Physics with Graph Networks. Xix:2202.009400 [physics, stat], September 2020. arXiv: 2002.09405. (*
Kimberly Stachenfeld, Drummond B. Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Efficient Turbulence Simulation. arXiv:2112, 1522:511 [physics], Jamuary 2022. arXiv: 2112.15275. 2* 5, 6,
Learned Coarse Models for
8, 10, 11, 17)
14

---

Published in Transactions on Machine Learning Research (04/2023)
Gilbert Strang. On the Construction and Comparison of Difference Schemes. SIAM Journal on Numerical Analysis, 5, 1968. doi: 10.1137/0705041. (5
Adam Subel, Ashesh Chattopadhyay, Yifei Guan, and Pedram Hassanzadeh. Data-driven subgrid-scale
modeling of forced Burgers turbulence using deep learning with generalization to higher Reynolds numbers
via transfer learning. Physics of Fluids, 3333)003, March 2021. ISSN 1070-6631, 1089-7666. doi:
10.1063/5.004025) arXiv:2012.0666: [physics]. 3*
Eitan Tadmor. The Exponential Accuracy of Fourier and Chebyshev Differencing Methods. SIAM Journal on Numerical Analysis, 23, 1986. doi: 00.133/0723001. 4)
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized Fourier Neural Operators. 1.X.::2111.1802 [cs], November 2021. arXiv: 2111.13802. 2*
Lloyd N. Trefethen. Spectral Melhods in MATLAB. Society for industrial and applied mathematics (STAM), 2000. (^ 2, 3, 4)
Kiwon Um, Raymond, Fei, Philipp Holl, Robert Brand, and Nils Thuerey. Solver-in-the-Love Learning from Differentiable Physics to Interact with Iterative PDE-Solvers. In NewrIPS 2020, June 2020. 7,
Rui Wang, Karthik Kashinath, Mustafa Mustafa, Adrian Albert, and Rose Yu. Towards Physics- informed Deep Learning for Turbulent Flow Prediction. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2020. 2)
David L Williarson.. Equivalent finite volume and Eulerian spectral transform horizontal resolutions estab- lished from aqua-planet simulations. Tellus Ser. A Dyn. Meteorol. Oceanogr., 60, January 2008. (A 12)
P K Yeung and K Ravikumar. Advancing understanding of turbulence through extreme-scale computation: Intermittency and simulations at large problem sizes. Phys. Rev. Fluids, 5, November 2020. (^ 4)
Z. Yin, H.J.H. Clerex, and D.C. Montgomery. An easily implemented task-based parallel scheme for the Fourier pscudospectral solver applicd to 2D Navier-Stoii turbulence. Computers & Fluids, 33(4):509-52) 2004. ISSN 0045-7930. doi: 10.1016/1 htt00///..omfffui....... 200.:/.....000000000..0. (^ 19)
Fisher Yu and Vladlen Koltun. Multi-Scale Context Aggregation by Dilated Convolutions. In Yoshua Bengio and Yann LeCun (eds.), 4th International Conference on Leorning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016. 18)
15


Reviews:
[{"review_id": "QyFPbX9rXj", "reviewer": "Reviewer_xrQa", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nI have major concerns or questions about the motivation, setting, presentation, and significance. The research problem is vague, the method motivation is vague, and much of the model is undefined. I’m also quite confused what this paper actually presents: is this a PDE solver that can be used to solve arbitrary PDE systems? \n\nI find the motivation of the paper to be vague. The introduction claims that coarse direct simulations can be improved by adding correction terms, which makes sense. But then we move to spectral methods, and it’s not clear why would we need corrections there. Isn’t the point of spectral methods that they are more precise, and thus have less need for corrections? What problem are we then solving? To me it seems that the key quantity is K: we only need corrections if the K is small enough. But why do we want to have small K, can’t we handle quite a lot of frequencies with traditional spectral method? The paper needs to make it clearer what is the research problem, and how does the K play into this.\n\nIn 2.1. the corrections are motivated for direct numerical methods, but then corrections are applied to spectral methods. I find this to be confusing. If the paper is about spectral methods, why are we talking about direct methods at all? I don’t see much relevance. Furthermore, I don’t see what are the issues in spectral methods that the correction fixes. Please state explicit what problems of spectral methods eq 3 solves.\n\nI wonder what is the motivation for including the correction step. The paper does not describe this much. I feel that adding the corrector acts as extra capacity or features that can augment eg. the insufficient 32 spectral solver. One could even see the DNN providing neural features that complement the spectral features used. This is obviously useful if the neural features are extracting different information than the frequencies (which they likely do). One could also think of the balance between frequencies and neural weights or features. For instance, if the neural network has 100K weights and we have 32 frequencies, why not instead have eg. 64 frequencies and 50K weights? Shouldn't we look at the total size of the parameterisation, and also study different balances between them? It would be good to discuss this. \n\nThe paper is written with heavy engineering focus, and the method is described in terms of neural network architectures and implementation tricks. This is not sufficient: the method needs to be presented using rigorous math throughout with proper motivations. Currently the method description basically stops at eq 3. I would add mathematical equations for the velocity solve, FFTs, neural netowrks, EPD and also for the SO-stepping.\n\nThe experiments show that the method can correct a spectral 32 model somewhat, but its still unable to make the solutions accurate. I wonder what is then the point of these experiments. Don’t we need accurate solutions? If the ML32 or ML64 is not accurate enough, why not do ML128 or ML256? I’m a bit confused overall what claims the experiments are trying to show.\n\nThe experiments also lack clear runtimes comparisons: I would expect that one goal is to show that eg. ML64 is faster than Spectral2048 (or something like that).\n\nFinally, I wonder what the contributions of this paper are. The listed contributions in introduction are vague: “we explore” “we compare” “some improvement”. What research problem do you solve? It seems that the paper’s contribution lies in providing some steps in making the ML parts of ML-augmented solvers better, even if the ML-augemnted solvers are not yet realistic. Whatever the claims and contributions are, they need to be more explicit.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a spectral PDE solver, which uses neural network corrections for improved efficiency. The method is interesting and potentially useful.\n", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "(See also the S&W)\n\n* Does K vary in fig1? Please make explicit.\n* The fig2 seems out of sync with rest of paper. Suddenly we have omegas instead of u-hats, theta is not present, and we have {velocities, IFFT, v, NN, FFT} emerging as components that have not been defined in the main text properly. Please present the complete method in math.\n* The NSO has not been introduced\n* The velocity stuff has not been defined, and I’m not sure where it comes from\n* The convnet padding is not explained clearly enough\n* The EPD is not explained\n* The convnet motivation could be strenghtened. CNNs are more flexible than differentials, and thus can model (slightly) non-local effects. Does this make sense for the underlying model?\n* The motivation for split operators is unclear. I don’t see why the spectral+correction method would need it. It’s also not clear what is referred to with the “each term of PDE”. Don’t we have just two terms? \n* In 3.2. it is stated that we use runge-kutta, while earlier we used euler. I think the eq 3 needs to describe the different steppings more clearly, or one should introduce a new equation that describes these in detail. Also it seems that “t\" is continuous, but we still use t+1 stepping. The method would then fail for systems that converge very quickly. Maybe “t” is just an index then?\n* Vorticity has not been explained\n* Can you describe how is the new solver used in practise. It seems that the neural correction is custom to the specific system. If the PDE changes, can we still apply the ML-solver? If not, does this then have some or any real-world use? Please explain\n* please include running times to the experiments. \n\n", "broader_impact_concerns": "No issues.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "kt3Lhwwy4v", "reviewer": "Reviewer_43tW", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n* Nice paper (but see writing weaknesses below)\n* Intuitive method\n* Results show promise\n\n**Weaknesses**\n\n* While I appreciate that the authors tried to succinctly provide all details within the 12 page limit, this paper doesn't seem like it is written for the typical ML reader. For example, I am only loosely familiar with several terms that are probably commonplace in the PDE literature:\n   - \"integrating forward in time\"\n   - \"simple explicit Euler time integrator\"\n   - the CFL condition\n   - a \"velocity-solve\" operation\n   - the mathematical details of the various operations that went into Figure 3\n   - etc\nwhich are all core components of the method, but not explained very well in this paper. Without precisely knowing what these mean this method is very hard for the layperson to reproduce. As such I feel this paper is either more suitable for JCP or similar, or has to be substantially rewritten in order to appeal to the broader TMLR audience. \n\n* In the same vein, for an ML audience it may be helpful to situate the performance results against more than just one ML-style benchmark. Right now the comparisons are with Stachenfeld et al 2022 (Fig 6), but this only paints a partial picture. How about other techniques (e.g. all the recent works on Fourier Neural Operators and their derivatives, which also learn representations in the spectral domain)? And if these are not appropriate comparison points, why not?. Especially, the 1D Burgers equation is a very common test case in most ML-based methods.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors study the problem of solving nonlinear PDEs using machine learning techniques. For certain families of PDEs, it is well known that spectral solvers are efficient compared to other numerical schemes; specifically, solving in the spectral (Fourier) domain enables resolving coarse-to-fine features. The contribution of this paper is to develop neurally augmented spectral solvers that learn a \"correction\" term to simulate errors caused at high resolution.\n\nThe authors showcase performance of their method on three classical PDE systems: 1D Kuramoto–Sivashinsky (KS) equation, 1D unstable Burgers’ equation, and 2D Kolmogorov flow.\n\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "+ More mathematical details (particularly on the PDEs parts) would help in terms of reaching out to a bigger section of the ML audience.\n+ A more thorough discussion about competing ML-based methods for PDE solutions, and perhaps more comparisons with such methods.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6FbX6c8FE3", "reviewer": "Reviewer_mxe5", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. From a big-picture perspective - the idea of having hybrid approaches that leverage existing strengths of numerical solvers and data-driven augmentations, to me, represents the pathway to greatest impact for scientific machine learning for partial differential equations. Yes there may be issues with computational costs when using the ML models, but data-driven corrections can be used for better shock-capturing or turbulence closure modeling tasks (as the authors mention in their discussion).\n\n2. The benchmark experiment on the Unstable Burgers' equation is timely and investigates the appropriateness of their proposed approach for potentially shocked solutions - something that has not been given the attention it deserves. Figure 5 tells me that their ML correction learns an effective hyperviscosity that prevents the explosion of their spectral baseline.\n\n3. A comparison across numerical discretization techniques is welcome - this is also something that isn't appropriately considered during most SciML studies for PDEs. While the conclusions from this study are quite well-known in the CFD community (that when applicable spectral methods outperform FVM), the experiment performed by the authors reveals some valuable quantitative insight.\n\nWeaknesses:\n1. While TMLR acceptance criteria are technical correctness and clarity of presentation rather than significance or impact, I am assuming there is a certain threshold of novelty which the authors must demonstrate before their submission is accepted. In that regard - I am curious how their work differs from existing literature that has been around for a while:\n- https://arxiv.org/pdf/2208.03680.pdf (where an NN learns a correction to an arbitrary numerical solve)\n- https://journals.aps.org/pre/abstract/10.1103/PhysRevE.97.042322 (where a shallow NN learns a correction to a reduced-order numerical solve)\nIn both these papers, ML corrections were seen to dramatically accelerate the solve of the numerical solver (by virtue of being able to use larger timesteps). The first example also has some interesting error analysis for bounding the total error due to the addition of the additive correction. \n\n2. Following on from point 1 - why not assess corrections that allow for larger timesteps during roll-outs? Wouldn't that make these methods much more attractive and also address the disappointing computational benefits? This wasn't clear in the paper - but I suspect the authors used the same time-step for roll-outs as for the training data. The main benefits of the ML additive correction would be realized if timesteps were significantly coarser (as demonstrated in the above references). In the context of their Burgers experiment - the authors can look at https://arxiv.org/pdf/2012.06664.pdf where a similar correction idea was used in conjunction with a physics-solve to enable much larger timesteps.\n\n3. The authors say - \"To understand the value of time-splitting for the learned correction, we also performed an experiment where we instead incorporated our ML model into the nonlinear part of Physics-Step\" - this is not clear. Did the authors use ML to correct the nonlinear term before it was sent into the time-integration? If so the behavior they observe tracks with past studies of the so-called structural closure modeling paradigm in turbulence - which has ultimately remained impractical because of similar (in)stability issues.\n\n4. Finally, this is not so much a weakness as a comment - the main consideration for using spectral versus FVM stems from issues other than merely accuracy (for example scaling on distributed computing hardware) - so while the proposed method and claims are sound for the problem identified (in the authors' words \"equations with smooth, periodic solutions\"), things would change dramatically for anything remotely practical. For example - their ConvNet architecture would require some degree of model-parallelism the moment a problem is distributed on multiple ranks. In that case, local corrections to a FVM scheme would be much more computationally tractable. \n\n5. The authors mention compute costs specific to the one ConvNet architecture they have used (for the NS problem). They have also mentioned that smaller architectures may make a difference about their claims but then double down by saying that \"there is little potential for accelerating smooth, periodic 2D turbulence beyond traditional spectral solvers\" - did they try using smaller architectures with fewer parameters? They note that much less compute was required to train this identical architecture in comparison to the EPD study. I was of the impression that if the heavy lifting is done by the physics-solve, you could get away with more efficient correction mechanisms.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a hybrid classical and data-driven technique for accurate time-integration of partial differential equations (PDE). Their approach builds a data-driven correction to a standard numerical integrator (here given by a spectral solver) and the authors observe gains in accuracy albeit at the cost of added computational expense. Their approach (the split operator) utilizes a baseline convolutional neural network (from a previous study) that learns local spatial correlations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Requested changes: \nI am unsure about requesting changes to this study. I would like to see their model perform for larger timesteps as compared to their baseline ground truth data and then provide some assessments for speed up - but as mentioned previously - we have seen previous studies demonstrate this. If the editor believes that the current work (despite my observations of somewhat insufficient novelty) is valuable to the TMLR audience - I have no reservations supporting its acceptance. \n", "broader_impact_concerns": "No concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers are learning towards accept: the main concern is 'incremental nature' of the approach (this is a valid claim), on the other hand, the method is new and functional. The paper is well-written and clear to follow. 
It would be interesting to see generalization and scaling from 2D turbulence to other classes of problems.
Decision: Accept (poster)