Neural Ideal Large Eddy Simulation: Modeling Turbulence with Neural Stochastic Differential Equations
Authors: Anudhyan Boral
Zhong Yi Wan
Leonardo Zepeda-Nunez
James Lottes
Qing Wang
Yi-Fan Chen
John Roberts Anderson
Fei Sha
Venue: NeurIPS 2023 poster (2023)
Keywords: partial differential equations
physics
turbulence
stochastic differential equations
physical simulation
neural differential equations
TL;DR: We synergize ideal large eddy simulation and neural stochastic differential equations for a probabilistic approach to turbulence modeling.

Abstract: We introduce a data-driven learning framework that assimilates two powerful ideas: ideal large eddy simulation (LES) from turbulence closure modeling and neural stochastic differential equations (SDE) for stochastic modeling. The ideal LES models the LES flow by treating each full-order trajectory as a random realization of the underlying dynamics, as such, the effect of small-scales is marginalized to obtain the deterministic evolution of the LES state. However, ideal LES is analytically intractable. In our work, we use a latent neural SDE to model the evolution of the stochastic process and an encoder-decoder pair for transforming between the latent space and the desired ideal flow field. This stands in sharp contrast to other types of neural parameterization of closure models where each trajectory is treated as a deterministic realization of the dynamics. We show the effectiveness of our approach (niLES – neural ideal LES) on two challenging chaotic dynamical systems: Kolmogorov flow at a Reynolds number of 20,000 and flow past a cylinder at Reynolds number 500. Compared to competing methods, our method can handle non-uniform geometries using unstructured meshes seamlessly. In particular, niLES leads to trajectories with more accurate statistics and enhances stability, particularly for long-horizon rollouts. (Source codes and datasets will be made publicly available.)

Neural Ideal Large Eddy Simulation: Modeling Turbulence with Neural Stochastic Differential Equations
Anudhyan Boral
Zhong Yi Wan Google Research Mountain View, CA 94043, USA
Google Research Mountain View. CA 94043, USA anudhyanggoogle.cmm com
wanzy0google. com
Leonardo Zepeda-Núñez
James Lottes
Google Research Mountain View, CA 94043, USA
Google Research Mountain View, CA 94043, USA
1zepedannuezz@google.. com
jlottes@google.com com
Qing Wang Google Research Mountain View, CA 94043, USA
Yi-fan Chen
Google Research Mountain View, CA 94043, USA
wqingdgoogle.com com
y aff.cen@googgee.coom com
John Roberts Anderson
Fei Sha
Google Research Mountain View. CA 94043, USA j janders@google.co com
Google Research
Mountain View, CA 94043, USA fsha@google.com com
Abstract
We introduce data-driven learning framework that assimilates two powerful ideas: ideal large eddy simulation (LES) from turbulence closure modeling and neural stochastic differential equations (SDE) for stochastic modeling. The ideal LES models the LES flow by treating each full-order trajectory random realization the underlying dynamics, as such, the effect of small-scales marginalized obtain the deterministic evolution of the LES state. However, ideal LES is analytically intractable. In our work, we use latent neural SDE to model the evolution of the stochastic process and encoder-decoder pair for transforming between the latent space and the desired ideal flow field. This stands in sharp contrast other types of neural parameterization of closure models where each trajectory is treated as deterministic realization of the dynamics. We show the effectiveness of our approach (niLES neural ideal LES) two challenging chaotic dynamical systems: Kolmogorov flow at Reynolds number of 20,000 and flow past cylinder Reynolds number 500. Compared competing methods, our method can handle non- non-unifon uniform geometries using unstructured meshes seamlessly. In particular, niLES leads trajectories with more accurate statistics and enhances stability, particularly for long- ang-horizon rollouts. (Source codes and datasets will be made publicly available.)
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Time
t=
1 +
wit s.t. wit vit
mean
DNS Space
LES Space
filtera
filter
x02
oft, Att)
Figure I: A cartoon of ideal LES. The LES field is lifted into the space real turbulent fields by applying the multi-valued inverse of the filtering operator. Then the turbulent DNS fields are evolved continuously according to the N-S equations. Finally, at the end of the LES time step, the mean the DNS fields is filtered to obtain the new LES field. Although the ideal LES is the ideal LES evolution, is analytically intractable. DNS image from [54].
Introduction
Multiscale physical systems are ubiquitous and play major roles in science and engineering [5, 7, 16, 25, 57, 74, 38]. The main difficulty of simulating such systems is the need to numerically resolve strongly interacting scales that are usually order of magnitude apart. One prime example of such problems turbulent flows, in which a fluid flow becomes chaotic under the influence of its own inertia. As such, high-fidelity simulations of such flows would require solving non-linear partial differential equations (PDEs) at very fine discretization, which often prohibitive for downstream applications due to the high computational cost. Nonetheless, such direct numerical simulations (DNS) are regarded as gold-standard [62, 43].
For many applications where the primary interest lies in the large (spatial) scale features of the flows, solving coarse-grained PDEs a favorable choice. However, due to the so-called back-scattering effect, the energy and the dynamics of the small-scales can have a significant influence the behavior the large ones [51]. Therefore, coarsening the discretization scheme alone results in highly biased (and often incorrect) characterization of the large-scale dynamics. To address this issue, current approaches incorporate the interaction between the resolved and unresolved scales by employing statistical models based on the physical properties of fluids. These mathematical models are commonly known as closure models. Closure models widely used approaches such as GyymoldsAvvergeee Navier-Stoke Stokes (RANS) and Large Eddy Simulation (LES) [68] have achieved considerable success. But they are difficult to derive for complex configurations of geometry and boundary conditions, and inherently limited in terms of accuracy [26, 67, 41].
An increasing amount of recent works have shown that machine learning (ML) has the potential to overcome many of these limitations by constructing data-driven closure models [55, 77, 47, 56). Specifically, those ML-based models correct the aforementioned bias by comparing their resulting trajectories to coarsened DNS trajectories as ground-truth Despite their empirical success and advantage over classical (analytical) ones, ML-based closure models suffer from several deficiencies, particularly when extrapolating beyond the training regime [9, 71, 79]. They are unstable when rolling the dynamics out to long horizons, and debugging such black box models is challenging. This raises the question: how do incorporate the inductive bias of the desired LES field to make good architectural choice of the learning model?
In this paper, we propose a new ML-based framework for designing closure models. We synergize two powerful ideas: ideal LES [50] and generative modeling via neural stochastic differential equations (NSDE) [80, 52, 44, 45]. The architecture of our closure model closely mirrors the desiderata of ideal

---

LES for marginalizing small-scale effects of inherent stochasticity in turbulent flows. To the best of our knowledge, our work is the first apply probabilistic generative models for closure modeling.
Ideal LES seeks a flow field of large-scale features such that the flow is optimally consistent with an ensemble of DNS trajectories that are filtered to preserve large-scale features [50]. The optimality, in terms of minimum squared error, leads to the conditional expectation of the filtered DNS flows in the ensemble. Ideal LES stems from the observation that turbulent flows, while in principle deterministic, are of stochastic nature as small perturbations can build up exponentially fast and after a characteristic Lyapunov time scale the perturbations erases the signal from the initial condition. Thus, while deterministically capturing the trajectory of the large-scales from single filtered DNS trajectory can be infeasible due to the chaotic behavior and loss of information through discretization, it may be possible to predict the statistics reliably over many possible realizations of DNS trajectories (with the same large-scale features) by marginalizing out the random effect. Fig illustrates the main idea and Section 2 provides detailed description. Unfortunately, ideal LES is not analytically tractable because both the distribution of the ensemble as well and the desired LES field are unknown.
To tackle this challenge, our main idea is to leverage i) the expressiveness of neural stochastic differential equations (NSDE) [52, 83] to model the unknown distribution of the ensemble as well as its time evolution, and ii) the power in Transformer-basee encoder-decoders to map between the LES field to-be-learned and a latent space where the ensemble of the DNS flows presumably resides. The resulting closure model, which we term as neural ideal LES (niLES), incorporates the inductive bias of ideal LES by estimating the target conditional expectations via Monte Carlo estimation of the statistics from the NSDEE-larmed distributions, which are themselves estimated from DNS data.
We investigate the effectiveness of niLES in modeling the Kolmogorov flow, which is classical example of turbulent flow. We compare our framework to other methods that use neural networks to learn closure models, assuming that each trajectory a deterministic realization. We demonstrate that niLES leads to more accurate trajectories and statistics of the dynamics, which remain stable even when rolled out to long horizons. This demonstrates the benefits of learning samples statistics (i.e., conditional expectations) from ensembles.
2 Background
We provide a succinct introduction to key concepts necessary to develop our methodology: closure modeling for turbulent flows via large eddy simulation (LES) and neural stochastic different equations (NSDE). We exemplify these concepts with the Navier-Stokes (N-S) equations.
Navier-Stokes and direct numerical simulation (DNS) We consider the N-S equations for in- compressible fluids without external forcing. In dimensionless form, the N-S equations are:
O.u (u V)u Vp vV with V
where = u(x, t) and p p(z, t) are the velocity and pressure of fluid at spatial point x in the domain Ω C Rd at time f; v is the kinematic viscosity, reciprocal of the Reynolds number Re, which characterizes the degree of turbulence of the flow. We may eliminate pressure p on the hight-hand-sider e.g. by taking the divergence of the momentum equation and using the fact that velocity divergence-free. Hence, we rewrite the N-S equations compactly as
8,u = R (1;
We impose boundary conditions on 0 and initial conditions u(z, 0) uo(x), x Ω in N-S equations throughout the manuscript. We sometimes omit the implicit y on the right hand side. To solve numerically, DNS will discretize the equation on fine grid G. such that all the scales are adequately resolved. It is important to note that as y becomes small, the inertial effects becomes dominant, thus requiring a refinement of the grid (and time- step due the Courant-Friedricee Lewy (CFL) condition [21]). This rapidly increases the computational cost [62].
LES and closure modeling For many applications where the primary interests are large-scalo features, the LES methodology balances computational cost and accuracy [76, 72]. It uses coarse grid G which has much fewer degrees-of-freedom than the fine grid G. To represent the dynamics with
respect to G. we define filtering operator (.) || + IR G which is commonly implemented
3

---

using low- (in spatial frequencies) filters [12]. Applying the filter to Eq. (2), we have
0,ū = R [1; v) + Relouure (u, u)
(3)
where ū is the LES field, Res has the same form as R with i7 as input, u is the DNS field in Eq. (2), and Re (ī7, u) R x Rl + R V (T गे Tæ) is the closure term. It represents the collective effect of unresolved subgrid scales, which are smaller than the resolved scales in g.
However, as unknown to the LES solver, the closure term needs to be approximated by functions of u. How to model such terms has been the subject of a large amount of literature (see Section 5). Traditionally, those models are mathematical ones; deriving and analyzing them is highly challenging for complex cases and entails understanding the physics of the fluids.
Learning-based closure modeling One emerging trend is to leverage machine learning (ML) tools to learn a data-driven closure model [47, 77] to parameterize the closure term,
Relosure (ū, u) R MM 0)
(4)
where 0 is the parameter of the learning model (often, a neural network). With a DNS field ground-truth, the goal is to adjust the 0 such that the approximating LES field
NS (i: v) + M(ū; 9)
(5)
matches the filtered DNS field u. This is often achieved through the empirical risk minimization framework in ML:
0" arg min Σ ||ū u.
where indexes the trajectories in the training dataset, each being DNS field from simulation with a different condition.
Despite their success, learning- based models have also their own drawbacks. Among them, how to choose the learning architecture for parameterizing is more of an art than a science. Our work aims to shed light on this question by advocating designing the architecture to incorporate the inductive bias of designed LES fields. In what follows, we describe ideal LES, which motivates our work. To give preview, our probabilistic ML framing matches very well the formulation of ideal LES in extracting statistics from turbulent flows of inherent stochastic nature.
Ideal LES It has long been observed that while chaotic systems, such as turbulent flows, can be in principle deterministic, they are stochastic in nature due to the fast growth (of errors) with respect to even small perturbation. This has led to many seminal works that treat the effect of small scales stochastically [68]. Thus, instead of viewing each DNS field as deterministic and distinctive realization, one should consider an ensemble of DNS fields. Furthermore, since the filtering operator (.) is fundamentally lossy [75], filtering multiple DNS fields could result in the same LES state. Ideal LES identifies evolution of LES field such that the dynamics is consistent with the dynamics of its corresponding (many) DNS fields. Formally, let the initial distribution T,, (u) over the (unfiltered) turbulent fields be fixed but unknown. By evolving forward the initial distribution according to Eq. (2), we obtain the stochastic process Te (u).
The evolution of the ideal LES field v is obtained from the time derivatives of the set of unfiltered turbulent fields whose large scale features are the same as v [50]:
dv Øt
Ou Ot
Er
Fig illustrates the conceptual framing of the ideal LES. We can gain additional intuition by observing that the field also attains the minimum mean squared error, matching its velocity P0/ôt to that of the filtered field Òu/Ot.
It is difficult to obtain the set {u = w) which is required to compute the velocity field (and infer w). Thus, despite its conceptual appeal, ideal LES analytic intractable. We will show how to derive data-driven closure model inspired by the ideal LES, using the tool of NSDE described below.

---

NSDE NSDE extends the classical stochastic differential equations by using neural-network pa- rameterized drift and diffusion terms [45, 80, 52]. has been widely used data-driven model for stochastic dynamical systems. Concretely, let time 1 E [0, 1], Ze the latent state and X the observed variable. NSDE defines the following generative process of data via a latent Markov process:
Zo ~ Po),), p()) ~ dZ: ho Zt. t) + go (z, t) e dW,, X2 p(,,
(8)
where Po(-) is the distribution for the initial state. W, is the Wiener process and denotes the Stratonovich stochastic integral. The Markov process (Zr)oe provides probabilistic prior for the dynamics, to be inferred from observations {Xe Note that the observation model p(X, Ze) only depends on the state at time t. h(-) and g(-) are the drift and diffusion terms, expressed by two neural networks with parameters 0 and Φ.
Given observation data x = { learning 0 and é is achieved by maximizing the Evidence Lower Bound (ELBO), under variational posterior distribution which is also an SDE [52]
(Z1[z)) ~ diz = (Ze, t, x x)dt + go go(Zi, t) o dW,
(9)
Note that the variational posterior has the same diffusion the prior. This is required to ensure finite ELBO, which given by
{ log( log(msl2)ddt
2 (
ho 2, () (o( go go t)
2
L C=Eg Eg
dt
(10)
Note that both the original SDE parameters 0 and Ø and the variational parameter 1 are jointly optimized maximize L. For detailed exposition of the subject, please refer to [45]. In the next section, we will describe how NSDE is used to characterize the stochastic turbulent flow fields.
3 Methodology
We propose neural SDE based closure model that implements ideal LES. The generative latent dynamical process in neural SDE provides a natural setting for modeling the unknown distribution of the DNS flow ensemble that crucial for ideal LES reproduce long-term statisties.
Setup We are given a set of filtered DNS trajectories {ū:}. Each V. is a sequence of "snapshots" indexed by time and spans a temporal interval Ta over the domain Ω. We use VARIABLE(() to denote the time / snapshot of the VARIABLE (such as ū, (1) or v(()) Those trajectories are treated as "ground- truth""tt"" and we would like to derive data-driven closure model M(u; Q) in the form of Eq. (3) to evolve the LES state
Oee= = R (v) + M(e),
(11)
where we have dropped the model's parameters ④ for notation simplicity. Our goal is to identify
the optimal M(r) such that the trajectories of Eq. (11) have the same long-term statistics as Eq. (2). Following ideal LES, we would render Eq. (11) equivalent to Eq. (7) by implementing the ideal M(v)
as
M(v(t)) = E.. ā.u गे v(t)] RNN (n(t)). Let Pe (u; denote the density of nH (u) restricted to the set {ulū e()):
(12)
(w,)(())) x 8(ā ((t))T((u)
(13)
We can thus rewrite the closure model
M(e(t)) / D.u - R. (u(t)) u(t) du / f(u; f(u:u(t)) du (14)
where fu) denotes the fluctuation exerted on the large scales v the true DNS trajectory was u. This shows that the ideal M(o(t) should compute the mean effect of the small-scale fluctuations fu) by integrating over all possible DNS trajectories u that are consistent with the large scales of u(t). However, just as ideal LES is not analytically tractable, so is this ideal closure model. Specifically, while the term R (c()) can be easily computed using a numerical solver on the coarse grid, the remaining terms are not easily computed. In particular, d.u would require a DNS solver, thus defcating the purpose of seeking a closure model. An approximation to Eq. (14) is needed.
5

---

Time
LES State Neural SDE Trajectories
Latent Space LES Space
ENC
DEC
Mean
Figure 2: Schematic of our modeling approach motivated from ideal LES (cf. Figure for structural correspondence). The evolution in the low-dimensional latent space follows trajectories of data- driven Neural SDE, mirroring the fine temporal resolution of the DNS trajectories. The final states of the latent state at the next time step are decoded into the LES space and averaged over obtain the mean correction due to small-scale fluctuations.
3.1
Main idea
Consider a trajectory segment from to to f where we are told only ū(ta) and ū(t1), how can we discover all valid DNS trajectories between these two times compute the closure model without incurring the cost of actually computing DNS? Our main idea is to use probabilistic generative model to generate representation of those trajectories learn from a corpus of DNS trajectories but the corpus is still too costly to acquire. Instead, leverage the observation that we do not need the full details the DNS trajectories order to approximate (v(()) well - the expectation itself is low-pass filtering operation. Thus, proceed by constructing (parameterized) stochastic dynamic process to emulate the unknown trajectories and collect the necessary statistics. As long the construeted process is differentiable, we can optimize it end to end by matching the resulting dynamics of using the closure model to the ground-truth.
fraction of the cost. One can certainly
We sketch the main components below. The stochastic process instantiated in the latent state space of neural SDE with an encoder whose output defines the initial state distribution controlled by the desired LES state at time to- The desired statistics, i.e., the mean effect small-scales, is computed in Monte Carlo estimation via a parameterized nonlinear mapping called decoder. The resulting correction by the closure model then compared to the desired LES state at time (, driving learning signals to optimize all parameters. See Fig. 2 for an illustration with details in Appendix D.
Encoder The encoder defines the distribution of the initial latent state variable in the NSDE, denoted by Zn E LR4, Concretely, RR b RL+L(L+1)2 maps from v(to) in the LES space the mean and the covariance matrix of a multidimensional Gaussian in the latent space: ((v(t)) = (/4z0) E This distribution used sample K initial latent states in RL: z N( E,,) 5 K).
Latent space evolution The latent stochastic process evolves according to time [0. 1]. This is
an important distinction from the time for the LES field. Since we are interested in extracting statistics from DNS with finer temporal resolution, 7 represents a faster stee-step evolving process whose beginning and end map to the physical time to and 41. (In our implementation, At = 41 to. the time step of the coarse solver, is greater than Ar, the time step for the NSDE.
The latent variable Z, evolves according to the Neural SDE:
dZ,, h2, d) go r) dW,
(15)
where W, the Wiener process on the interval [0, 1]. We obtain trajectories ( (2 }rejo,1 sampled from the NSDE, and in particular, we obtain an ensemble {z(")) 12-1-
Decoder The decoder D RL + from latent state back into the LES space. So can compute the Monte-Carlo approximation of
x maps each the K ensemble members {

---

Algorithm 1 Compute M(v((o))
Input: LES state (()), Closure model parameters O
1. Encode uto)) to a distribution on the latent state: (H2u 22o) = &a(u(to)).
2. Sample K initial latent states in RL; Z NN 22o) (1 S i S K).
3. Sample trajectories - (T E (0, 1D by solving Eq. (15) with initial conditions 2
4. Decode the final latent states _ to the LES space: x(() = De (2 (1 S < K).
5. Take the empirical mean of the samples: v' ( E x(0) Output: v' /tt - to).
M(v(to))(tz to) (cf. Eq. (14) for the definition)
() d'd ((=
vo f fdw (twlwww umm 2
D(Z 1=1
K
(16)
where w is the spatio-temporal lifted version of u, i.e., the space of DNS trajectories in space and time, and fw absorbs both in the closure Eq. (14) and the implicit conditioning in Eq. (13), and pw is a prior for the trajectories. This notational change allows us to approximate the integral directly by a Monte- Carlo approximation on the lifted variables, i.e., the distribution of trajectories (as shown in Fig. 2) which are modeled using the NSDE in Eq. (15). See Algo. 1 for the calculation steps.
We stress that while motivate design via mimicking DNS trajectories, Z is low-dimensional and is not replicating the real DNS field u. However, it is possible that with enough training data, Z might discover the low-dimensional solution manifold in u.
Training objective The NSDE is differentiable. Thus, with the data closure model Eq. (16), apply end-to-end learning match the dynamics with the closure model to ground-truths Concretely, let vto) be ū(t)) and we numerically integrate Eq. (11)
utt1) v(to) + J." R (u(t)dd + ((u((t)))(t - to)
(17)
This gives rise the likelihood the (observed) data in NSDE. Specifically, following the VAE setup of [52], we have
- logg((((11)) = L= ((1)) (20*) ( - 21
(18)
where is scalar posterior scale term. The training objective L(v, u) L (v, u) + KI. NSDE includes a KI. divergence term associated with the neural SDE. See Appendix D for more details on the additional term.
To endow more stability to the training, following [8]] the training loss incorporates S time steps (S > 1) of rollouts of the LES state:
(3)) Σ L(olto kAI), Tlto kAt))
(19)
For a dataset with multiple trajectories, we just sum the loss for each one of them and optimize through stochastic gradient descent.
3.2 Implementation details
We describe details in the Appendix G for the Transformer-bssee encoder and decoder and Ap- pendix D for implementing NSDE with a Transformer parameterizing the drift and diffusion terms.
Experimental results
We showcase the advantage of our approach using two instances of chaotic Navier-Stokes flow: 2D Kolmogorov flow, [47] at Reynolds number of 20,000 and flow past circular cylinder at Reynolds number of 500.
7

---

4.1 Setup
Dataset generation The reference data for Kolmogorov flow consists of an ensemble of trajectories generated by randomly perturbing an initial condition. The flow past cylinder dataset is generated as single long, chaotic trajectory, split up into 14 equal segments with an initial subsegment thrown away to avoid correlation between the segments. Fach trajectory is generated by solving the NS equations directly using a high-order spectral finite element discretization in space and time-stepping is done via 3rd order backward differentiation formula (BDF3) [23], subject to the appropriate (CFL) condition [21]. These DNS calculations were performed on a mesh with 2304 elements for Kolmogorov flow and 292 elements for the cylinder flow with each element having a polynomial order and 13 respectively. The DNS trajectories were then sampled to a coarse grid with order and 5 respectively for the two datasets. For more details see Appendix B.
Benchmark methods For all LES methods, including niLES, we use a 10x larger time step than the DNS simulation. We compare our method against a classical implicit LES at polynomial order 4 using a high order spectral element solver using the filtering procedure of [28]. As reported in [14], the implicit LES for high order spectral methods is comparable to Smagorinsky subgrid-scale eddy-viscosity models. We also train an encoder-decoder deterministic NN-based closure model using similar transformer backend our niLES model. This follows the procedure from prior works [56, 81, 47, 24] in training deterministic NN-based closure models. Additionally, we train deterministic poodde-processor-eeceeee network with neural ODE processor [19] with strictly more parameters than the niLES model.
Metrics For measuring the short-term accuracy we used root-mean- root-mean-square--erro (RMSE) after unrolling the LES forward. Due to chaotic divergence, the measure becomes meaningless after about 1000 steps. For assessing the long-tery tergttemm&e ability to capture the statistics of DNS we use turbulent kinetic energy (TKE) spectra. The TKE spectra is obtained by taking the 2D Fourier transform of the velocity field, after interpolating it to uniform grid. Then we obtain the kinetic energy the Fourier domain and we integrate it in concentric annuli along the 2D wavenumber k.
4.2 Main Results
Long-term turbulent statisties We summarize in Figs. 3 and 5. In the short-term regime, both niLES and deterministic methods achieve higher accuracy than the implicit LES, even when unrolling for hundreds of steps. Beyond this time frame, the chaotic divergence leads to exponential accumula- tion of pointwise errors, until the LES states are fully decorrelated with the DNS trajectory. At this point, we must resort to statistical measures to gauge the quality of the rollout.
When considering reconstruction errors from very short-term rollouts, such as when the model unrolled only times, identical to the training setup, the deterministic approaches and niLES yield similar reconstruction errors. Beyond this regime, and going up to 1000s of rollout steps, the probabilistic approach of niLES provides superior generalization ability even when measuring a deterministic reconstruction loss.
Furthermore, niLES captures long-term turbulent statisties significantly better than the other two approaches, particularly in the high wavenumber regime. The deterministic NN is not stable for long term rollouts due to energy buildup in the small scales, which eventually leads the simulation to blow up.
Inference costs The cost of the inference for our method well the baselines are summarized in Table 1. niLES uses four SDE samples for both training and inference, and each SDE sample resolved using 16 uniformly-spaced time steps corresponding to a single LES time step. The inference cost scales linearly with both the number of samples and the temporal resolution. However, niLES achieves much lower inference cost than the DNS solver while having similarly accurate turbulent statistics for the resolved field. The deterministic NN model has slightly lower inference cost than niLES, since can forego the SDE solves. While the implicit LES is the fastest method and is long-term stable, it cannot capture the statistics especially in the high wavenumber regime.
Limitations A drawback of our current approach stems from using transformers for the encoder- decoder phase of our model, which might result in poor scaling with increasing number of mesh
8

---

3.1
InpiiE
016
ImpEce LES
227
nILES
mILES
0.00
200
110
Inpliii 1P%
0102
331
Draemr
Deterres
PEPS
300
300
1500
0.00
1000
Figure 3: Root mean squared errо (RMSE) over the first 1000 steps (first two columns) and the turbulent kinetic energy (TKE) spectrum ((k) averaged over the first 2500 steps (right two columns) of two independent test trajectories unseen during training or validation. niLES has an improved ability capture the long term statisties accurately compared to both implicit LES and deterministic NN. The energy buildup in the small scales (large wavenumber) in the deterministic NN model eventually leads unstable trajectories.
04
65
XX
Figure 4: Comparison between rollout predictions after 800 LES steps on a held-out trajectory. Velocities in the æ (top row) and y (bottom row) directions respectively. Snapshots of filtered DNS (reference) (a), niLES (b), implicit LES (c) and deterministic NN models (d). The niLES captures several finer scale features of the flow consistent with the reference filtered DNS trajectory. The implicit LES has an overall smoothing effect and some turbulent structures are not captured. The deterministic NN LES shows artifacts which indicate instability.
elements. Alternative architectures which still handle interacting mesh elements should be explored. Additionally, the expressibility of the latent space in which the solution manifold needs to be embedded can affect the performance of our algorithm, and requires further study.
Related work
The relevant literature on closure modeling is extensive it is one of the most classical topics in computational methods for science and engineering [68]. In recent years, there has also been an explosion of machine learning methods for turbulence modeling [9, 25] and multi-scaled systems in general. We loosely divide the related works into four categories, placing particular emphasis on the treatment of effects caused by unresolved (typically small-scaled) variables.
Classical turbulence methods primarily relies on phenomenological arguments to derive an eddy viscosity term [48], which is added to the physical viscosity and accounts for the dissipation of energy from large to small scales. The term may be static [4], time-dependent [76, 31] or multi-scale [39, 40].
Data-driven surrogates often do not model the closure in explicit way. However, by learning the dynamics directly from data finite resolution, the effects of unresolved variables and scales are expected to be captured implicitly and embedded in the machine learning models. A variety of architectures have been explored, including ones based on multi-scaled convolutional neural networks [70, 82, 77], transformers [1 1], graph neural networks [73, 49] and operator learning [65].

---

Deleminng::
WeJratDDE
0.10
0.24
0.30
0.00
43
Figure 5: Average RMSE over Figure 6: Comparison between rollout predictions of the 10 testcases on the cylinder wake dataset. Shaded regions depict the
horizontal component of velocity in the cylinder wake spread (minimum and maximum) dataset after 100 LES steps on a held-out testcase. From among the testcases.
top bottom: snapshots of (a) filtered DNS (reference), (b) niLES, and (c) ILES, The niLES captures correct shapes of the vortices when compared ILES. The latter has an excessive smoothing effect and the vortices appear more rounded.
Table 1: Inference times in wall clock seconds evolve the state for second of simulation time for the Kolmogorov dataset
Method
DNS
Implicit LES
Deterministic NN
niLES (Ours)
Inference time [s]
15600
8.85
11.4
23.8
Hybrid physics-ML contains a rich set of recent methods to combine classical numerical schemes and deep learning models [61, 6, 47, 56, 24, 81, 58, 35]. The former is expected provide reasonable baseline, while the latter specializes in capturing the interactions between modeled and unmodeled variables that accurately represent high-resolution data. This yields cost-effective. low-resolution methods that achieve comparable accuracy to more expensive simulations.
Probabilistic turbulence modeling seeks to represent small-scaled turbulence stochastic pro- cesses [1, 32, 33, 20, 36). Compared to their deterministic counterparts, these models are better equipped to model the backscattering of energy from small large scales (i.e. opposite to the direction energy flow eddy viscosity), which rare in occurrence but often associated with events of great practical interest and importance (e.g. rogue waves, extreme wildfires, etc.).
Our proposed method inspired by methods the last two categories. The closure model we seek includes the coarse solver in the loop while using probabilistic generative model to emulate the stochastic process underlying the turbulent flow fields. This gives rise to long-term statistics that are accurate the inexpensive neural SDE provides good surrogate for the ensemble of the flows.
6 Conclusion
Due to chaotic divergence, is infeasible to predict the state the system with pointwise accuracy. Fortunately, in many systems of interest, the presence an attractor allows much lower-dimensional model to capture the essential statistics of the system. However, the time evolution of the attractor is known, which makes building effective models challenging.
In this work we have argued that taking probabilistic viewpoint is useful when modeling chaotic systems over long time windows. Our work has shown modeling the evolution with a neural SDE is beneficial in preserving long-term statistics and this line of thinking is likely fruitful.
10

---

References
[1] Ronald Adrian. On the role of conditional averages in turbulence theory. Tarbulence in Liquids, pages 323-332. 1977.
[2] Ronald Adrian. Stochastic estimation of the structure of turbulent fields. Springer, 1996.
[3] Ronald Adrian, BG Jones, MK Chung, Yassin Hassan, CK Nithianandan, and AT-C Tung. Approximation of turbulent conditional averages by stochastic estimation. Physics of Fluids A: Fluid Dynamics, 1(6)9922- 998, 1989.
[4] Nadine Aubry, Philip Holmes, John I. Lumley, and Emily Stone. The dynamics of coherent structures in the wall region of turbulent boundary layer. Journal of fluid Mechanics, 192:115 122111-177, 1988.
[5] Gary S Ayton, Will G Noid, and Gregory A Voth. Multiscale modeling of biomolecular systems: in serial and in parallel. Current opinion in structural biology. 17(2):192-198) 2007.
[6] Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116631):15344-55 -5349999444499 2019.
[7] Peter Bauer. Alan Thorpe, and Gilbert Brunet. The quiet revolution of numerical weather prediction. Nature, 525(7567):47-55, 2015.
[8] Andrea Beek, David Flad, and Claus-Dieter Munz. Deep neural networks for data-driven les closure models. Journal of Computational Physics, 398:108910, 2019.
Andrea Beck and Marius Kurz. perspective machine learning methods turbulence modeling. CAMM-Minnilungen, 44(1 1421):220100002,, 2021.
[10] G Berkooz. An observation on probability density equations, or, when do simulations reproduce statistics? Nonlinearity, 7(2):313, 1994.
[11] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Pangu-weather: A 3d high-resolution model for fast and accurate glohal weather forecast. arXiv preprint arXiv:2211.02556, 2022.
[12] Hugh M Blackburn and S Schmidt. Spectral element filtering techniques for large eddy simulation with dynamic estimation. Journal of Computational Physics, 186(2):610-62,, 2003.
[13] Guido Boffetta and Robert E. Ecke. Two-dimensional turbulence. Annual review of fluid mechanics, 44:427-451, 2012.
[14] Christoph Bosshard, Michel O Deville, Abdelouahab Dehbi, Emmanuel Leriche, et al. Udns or les, that is the question. Open Journal of Fluid Dynamics, 5(04):339, 2015.
[15] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclau- rin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne et al. Jax: composable transformations of python+ numpy programs. 2018.
[16] Bradshaw. Turbulence modeling with application to turbomachinery. Progress in Aerospace Sciences. 32(6):575-624, 1996.
[17] Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. arXiv preprint arXi::2202.03370 2022.
[18] Ashesh Chattopadhyay, Jaideep Pathak, Ebrahim Nabizadeh, Wahid Bhinji. and Pedram Hassanzadeh. Long-term stability and generalization of obscrvationally-coom constraidd stochastic data-driven models for geophysical turbulence. Environmental Data Science, 2:e1, 2023.
[19] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[20] Carlo Cintolesi and Etienne Mémin. Stochastic modelling turbulent flows for numerical simulations. Fluids, 5(3):108, 2020.
[21] R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathematical physics. IBM Journal of Research and Development, 11(2):215-25 234, 1967.
15

---

[22] Daan Crommelin and Eric Vanden-Eijuden, Subgrid-scale parameterization with conditional markov chains. Journal of the Atmospheric Sciencex, 65(8)226612265, 2008.
[23] M. O. Deville, P. F. Fischer, and E. H. Mund. High-Order Methods for Incompressible Fluid Flow. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2002.
[24] Gideon Dresdner, Dmitrii Kochkov, Peter Norgaard, Leonardo Zepeda-Núñez, Jamie Smith, Michael P Brenner, and Stephan Hoyer. Learning to correct spectral methods for simulating turbulent flows. arXiv preprint arXii:2207.005555 2022.
[25] Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao. Turbulence modeling in the age of data. Anmual review of fluid mechanics, 51:357-377, 2019.
[26] Paul A Durbin. Some recent developments in turbulence closure modeling, Annu. Rev. Fluid Mech., 50, 2018.
[27] Haoqi Fan, Bo Xiong, Karttikeya Mangalam, Yanghao Li, Zhicheng Yan, Jitendra Malik, and Christoph Feichtenhofer. Multiscale vision transformers. In Proceedings of the IEFF//VF International Conference on Computer Vision, pages 6824-6835, 2021.
[28] Paul Fischer and Julia Mullen. Filter-based stabilization of spectral element methods. Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 332(3):265-270, 2001.
[29] Kai Fukami, Koji Fukagata, and Kunihiko Taira. Super-resolution reconstruction of turbulent flows with machine learning. Journal of Fluid Mechanics, 870:106-1220 120, 2019.
[30] Kai Fukami, Koji Fukagata, and Kunihiko Taira. Machine-lcarming-aase based spatio-temporal super resolution reconstruction of turbulent flows. Journal of Fluid Mechanics, 909:A9, 2021.
[31] Massimo Germano, Ugo Piomelli, Parviz MOin, and William H Cahot. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids A: Fluid Dynamics, 3(7):1760-1765, 1991.
[32] Ian Grooms and Andrew J Majda. Efficient stochastic superpaarametterizationn for geophysical turbulence. Proceedings of the National Academy of Sciencex, 110(12):44444466, 2013.
[33] Ian Grooms and Andrew Majda. Stochastic superparamettriiattoon quasigeostrophic turbulence. Journal of Computational Physics, 271:78-98, 2014.
[34] Yifci Guan, Ashesh Chattopadhyayy Adam Subel, and Pedram Hassanzadeh. Stable a posteriori les of 2d turbulence using convolutional neural networks: Backscattering analysis and generalization to higher re via transfer learning. Journal of Computational Physics, 458: 1090, 2022.
[35] Yifei Guan, Adam Subel, Ashesh Chattopadhyay. and Pedram Hassanzadeh. Learning subgrid-scale closures in the small-data regime for stable and accurate les. Physica D: Nonlinearr Phenom- ena, 443:133568, 2023.
[36] Arthur Guillaumin and Laure Zanna. Stochastic-decy learning parameterization ocean momentum forcing. Journal of Advances in Modeling Earth Systems, 13(9): 2021.
[37] Irina Higgins, Loie Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner. beta-vae: Learning basic visual concepts with constrained variational framework. In International conference on learning representutions, 2017.
[38] Thomas JR Hughes, Gonzalo R Feijóo, Luca Mazzei, and Jean- BapptBtt Quincy. The variational multi- scale method a paradigm for computational mechanies. Computer methods in applied mechanics and engineering, 166(1-2):3-2 24, 1998.
[39] Thomas JR Hughes, Luca Mazzei, and Kenneth E Jansen. Large eddy simulation and the variational multiscale method. Computing and visualization in science, 3:47-59, 2000.
[40] Thomas JR Hughes, Assad Oberai, and I Luca Mazzei, Large eddy simulation of turbulent channel flows by the variational multiscale method. Physics of flidd, 1746:198 2001.
[41] Javier Jimenez and Robert D Moser. Large-eddy simulations: where are we and what can we expect? AIAA journal, 6166126 2000.
[42] Myeongseok Kang, Youngmin Jeon, and Donghyun You. netual- mixed subgrid-scne scale model for turbulent flow. Journal of Fluid Mechanics, 962:A38, 2023.
16

---

[43] Petr Karnakov, Sergey Litvinov, and Petros Koumoutsakos. Computing loaming flows across scales: From breaking waves microfluidics. Science Advances, 86))eabm05500 2022.
[44] Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.024335 2022.
[45] Patrick Kidger, James Foster, Xuechen Li, and Terry J Lyons. Neural sdes International Conference on Machine Learning, pages 5453-5463. PMLR, 2021.
gans. In [46] Patrick Kidger, James Foster, Xuechen Chen Li, and Terry Lyons. Efficient and accurate gradients for neural sdes, Advances Neural Information Processing Systems, 34:18747-18761, 2021.
[47] Dmitrii Kochkov. Jamie A Smith, Ayya Alieva, Qing Wang, Michael Brenner, and Stephan Hoyer. Machine learning-accelrrated computational fluid dynamics. Proceedingy of the National Academy of Sciences, L18(21))e2101784115 2021.
[48] Andrei Nikolaevich Kolmogorov. The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers. Proceedings of the Royal Society London. Series A: Mathematical and Physical Sciences, (199))9)),-,3, 1991.
[49] Remi Lam, Alvaro Sanchez- Matthew Willson, Peter Wirnsberger, Meire Fortunato, Alexander Pritzel, Suman Ravuri, Timo Ewalds, Ferran Alet, Zach Eaton-Rosen, et al. Graphcast: Learning skillful medium-range global weather forecasting. arXiv preprint arXiv:2212. 12794, 2022.
[50] Jacob A Langford and Robert D Moser. Optimal les formulations for isotropic turbulence. Journal of fluid mechanics, 398:321-346, 1999.
[51] CE Leith. Stochastic backscatter in a subgrid-scale model: Plane shear mixing layer. Physics of Fluids A: Fluid Dynamics, 2(3):297-299 1990.
[52] Xuechen Li, Ting-Kam I Leonard Wong, Ricky TQ Chen, and David K Duvenaud. Scalable gradients and variational inference for stochastic differential equations. In Symposium Advances in Approximate Bayesian Inference. pages 1-28. PMLR, 2020.
[53] Yanghao Li, Chao-Yuan Wu, Haoqi Fan, Karttikeya Mangalam, Bo Xiong, Jitendra Malik, and Christoph Feichtenhofer. Mvitv2: Improved multiscale vision transformers for classification and detection. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4804 4814, 2022.
[54] Yi Li, Eric Periman, Minping Wan, Yunke Yang, Charles Meneveau, Randal Burns, Shiyi Chen, Alexander Szalay, and Gregory Eyink. A public turbulence database cluster and applications study lagrangian evolution of velocity increments in turbulence. Journal of Turbulence, 9():N31, 2008
[55] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv:2010.000 08895 [cs, math), May 2021. arXiv: 2010.08895.
physics- based loss functions and optimisation horizons. Journal of Fluid Mechanics, 949:A25, 2022. [56] Bjöm List, Li-Wei Chen, and Nils Thuerey. Learned turbulence modelling with differentiable fluid solvers:
[57] Peter Lynch. The emergence of rumerical weather prediction: Richardson's 's dream. Cambridge University Press, 2006.
[58] Jonathan F MacArt, Justin Sirignano, and Jonathan Freund. Embedded training of neural-network subgrid-scale turbulence models. Physical Review Fluids, 6(5):050502, 2021.
[59] Romit Maulik, Kai Fukami, Nesar Ramachandra, Koji Fukagata, and Kunihiko Taira. Probabilistic neural networks for fluid flow surrogate modeling and data recovery. Physical Review Fluids, 5(10):00440,, 2020.
[60] Romit Maulik, Omer San, Jamey Jacob, and Christopher Crick. Sub-grid scale model classification and blending through deep learning. Journal of Fluid Mechanics, 870:784-812, 2019.
[61] Siddhartha Mishra. A machine learning framework for data driven acceleration of computations of differential equations. arXiv preprint arXiv:1807.09511 2018
[62] Steven A. Orszag. Analytical theories of turbulence. Journal of Fluid Mechanicy, 41(2):363-386 1970.
[63] Shaowu Pan and Karthik Duraisamy. Data-driven discovery of closure models. SIAM Journal on Applied Dynamical Systems, 17(4):2381 1211, 2018.
17

---

[64] Jonghwan Park and Haecheon Choi. Toward neuura-network-basee large eddy simulation: Application to turbulent channel flow. Journal of Fluid Mechanicx, 914:A16, 2021.
[65] Jaideep Pathak, Shashank Subramanian, Peter Harrington. Sanjeev Raja. Ashesh Chattopadlyay. Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li. Kamyar Azizzadenesheli, et al. Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv: 1 /224,, 2022.
[66] Stephen B Pope. Pdf methods for turbulent reactive flows. Progress in energy and combustion science, 11()):1191192 1985.
[67] Stephen B Pope. Ten questions concerning the large-eddy simulation of turbulent flows. New Journal of Physics, 6:35-35, March 2004.
[68] Stephen B Pope and Stephen B Pope. Turbulent flows. Cambridge university press, 2000.
[69] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. In International Conference on Machine Learning, pages 5301-5310. 10. PMLR, 2019.
[70] Olaf Ronneherger, Philipp Fischer, and Thomas Brox, U-net: Convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted InterventionMMICCAl 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234-241. Springer, 2015.
[71] Andrew Ross, Ziwei Li, Pavel Perezhogin, Carlos Femmandez-Grandee and Laure Zanna. Benchmarking of machine learning ocean subgrid parameterizations in an idealized model. Journal of Advances in Modeling Earth Systems, ((1)::2222200033 2023.
[72] Pierre Sagaut. Large eddy simulation for incompressible flows: an introduction. Springer Science & Business Media, 2005.
[73] Alvaro Sanchez- Sanchez-Gonzalez. Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International conference on machine learning, pages 8459-8468. PMLR, 2020.
[74] Tapio Schneider, João Teixeira, Christopher S Bretherton, Florent Brient, Kyle G Pressel, Christoph Schär, and A Pier Siehesma. Climate goals and computing the future of clouds. Nature Climate Change, 7(1):35,, 2017.
[75] Claude Elwood Shannon. Communication in the presence of noise. Proceedings of the Institute of Radio Engineers, 37(1):10-21, 1949.
[76] Joseph Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiment. Monthly weather review, 91(3):99-166 164, 1963.
[77] Kimberly Stachenfeld, Drummond B Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho. Peter Battaglia, and Alvaro Sanchez- Gonzalez. Learned coarse models for efficient turbulence simulation. arXiv preprint arXiv:2112.15273 2021.
[78] Adam Subel. Ashesh Chattopadhyay, Yifei Guan, and Pedram Hassanzadeh. Data-driven subgrid- scale modeling of forced burgers turbulence using deep learning with gemeralization to higher reynolds numbers via transfer learning. Physics of Fluids, 33(3):331 1000311700 100,, 2021.
[79] Salar Taghizadeh, Freddie D Witherden, and Sharath S Girimaji. Turbulence closure modeling with data-driven techniques: physical compatibility and consistency considerations. New Journal of Physics, 2099023 2020.
[80] Belinda Tzen and Maxim Raginsky. Neural stochastic differential equations: Deep latent gaussian models in the diffuuion limit. arXiv preprint /90.:::905.099000 2019.
[81] Kiwon Um, Robert Brand, Yun Raymond Fei, Philipp Holl, and Nils Thuerey. Solver-in-the-loo: Learning from differentiable physics to interact with iterative pde-solvers. Advances Neural Information Processing Systems, 33:61 13:6111-611 16122 2020.
[82] 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, pages 1457-1466, 2020.
[83] Winnie Xu, Ricky TQ Chen, Xuechen Li, and David Duvenaud. Infinitely deep bayesian neural networks with stochastic differential equations. In International Conference on Artificial Intelligence and Statistics. pages 721-738. PMLR, 2022.
18


Reviews:
[{"review_id": "AvqZaHHJ9x", "reviewer": "Reviewer_SnNJ", "summary": "This paper targets learning turbulence closure models for RANS simulations via neural networks. The paper proposes to use a neural SDE on the latent space of a transformer to predict different samples from the distribution of the next state, and then compute an average over these. This process is unrolled and trained for a sequence of multiple steps.", "questions": "Which other, specific applications and scenarios do the authors see for their method?", "limitations": "Limitations are discussed briefly.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Overall, this is a good idea, and I'm not aware of a an NSDE being previously used in this form. Thus, I see the general direction of the paper and the promising approach as strong points.", "weaknesses": "On the other hand, the paper targets a single, two dimensional Kolmogorov flow secnario as the only test case. In addition, only a single deterministic NN is compared to (plus an implicit LES solver). For this single data set, the paper is lacking a stable evaluation: multiple, differentily initialized models evaluated across multiple tests to obtain a stable result are not evaluated. In addition, NeurIPS is a very broad ML venue. Turbulence is definitively an exciting topic here, but even more important would be a broader evaluation, ideally with substantially different secnarios to show that the method has merit beyond turbulence. In its current form, I don't think that the results are sufficient for a NeurIPS paper.\n\nI see two ways to improve this aspect of the submission: either the authors focus their writing on the turbulence secnario (cf. below), and present multiple scenarios in this context, or non-turbulence cases are included to broaden the scope. A more stable evaluation with multiple models and tests should be included either way. This potentially also could help to show the benefits of the method more clearly. Right now the gains in terms of accuracy and the differences in the TKE spectrum seem to be mild. Additional cases could show areas where the approach gives larger improvements.\n\nIn addition, I would also recommend that the authors include additional learned baselines. This is less crucial, but would nonetheless help to put the work into the context of previous methods at NeurIPS, ICLR & co.\n\nI also do want to mention two weak points in the writing. One is that I found the motivation (esp. L44) quite unintuitive: the \"inductive bias of the LES field\" is not very clear, and the summary implies this is \"simply\" a matter of choosing the right architecture. Things get clearer afterwards, but reading the paper front to back, I think this summarizing question is not helping a reader.\n\nThe conclusions are also not a good fit with the rest of the paper: suddently, a transition is is made to generic chaotic systems. The whole previous paper targets a single specific scenario in the form of Kolmogorov turbulence, and presents a single set of results. Hence, this outlook is not supported by the content of the paper. I can understand that the authors have hopes that their method will at some point in the future generalize to other cases, but it should be made clear that this is an outlook. Rather, references to specific works where the authors see potential would be interesting to give here.\n\nOverall, I want to encourage the authors to continue their direction of work. Nonetheless, I find it difficult to directly argue for accepting this paper in its current form.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review and detailed feedback. We appreciate your encouragement and your suggestions on how to improve the quality of the manuscript.\n\n**More scenarios**\n\nWe have included in the uploaded PDF an instance of our proposed methodology applied to cylinder wake at a high Reynolds number of 500, where the system exhibits irregular vortex shedding and chaotic flow. The difficulty of such an example is two-fold [3]: the need of a nonuniform mesh that is usually refined near the cylinder, and the need to capture the boundary layer near the cylinder. We hope this adds diversity to the turbulent flows that this method could be applied to. See Fig 4 in the uploaded PDF for a qualitative rollout and Fig 5 for the accuracy plots with the new dataset. We have also shown the spread (min/max values) among 10 testcases in the dataset.\n\n**Stable Evaluation**\n\nThank you for the suggestion. For both the Kolmogorov datasets, we have added evaluations across 6 different test cases (See Fig 1 and Fig 2 in the attached PDF) We plot both the RMSE and the TKE among the 6 test cases. While our training method unrolls for only 8 steps, the evaluation is unrolled to 100s of steps, well beyond the training horizon.\n\n**Why turbulent flows?**\n\n> I see two ways to improve this aspect of the submission: either the authors focus their writing on the turbulence secnario (cf. below), and present multiple scenarios in this context, or non-turbulence cases are included to broaden the scope\n\nThank you for the comment. We target turbulent flows as a prototypical chaotic system with an impact in real world applications encompassing engineering and science, such as weather and climate.  In fact, many of the difficulties inherent to chaotic systems manifest in turbulent flows. For example, maintaining long-term statistics accurately (such as turbulent kinetic energy) in turbulent flows is especially challenging. Many data-driven approaches become unstable for longer rollouts [1,2]. Your suggestion of investigating the method’s applicability to other flows is appreciated and as mentioned above we have added a new setup for such flows. \n\nWe will modify the text in order to convey this point more clearly.\n\n**Gains in Accuracy and TKE**\n\nWhile the gains in accuracy shown in Fig. 3 (of the manuscript) are seemingly modest, achieving similar gains using traditional computational techniques used in turbulent flows would require the computation of an expensive direct numerical simulation (DNS). For instance, using a high-order DNS that achieves a similar quality is roughly two orders of magnitude slower. \n\nTo the best of our knowledge, similar works on Neural Network-based turbulence closure models do not achieve such gains at such a high Reynolds number (20,000) while being stable for long beyond the training horizon. Furthermore, among spectral element methods, the Implicit LES method is considered to be the state-of-the-art [4]. Even against this approach we have demonstrated that our method produces more accurate simulations while exhibiting long-term stability.\n\n**Additional baselines**\n\nThank you for the suggestion. We are including a Neural ODE model [5] as a learned baseline for deterministic LES.\n\nAs you pointed out, one of the novel aspects of our approach is leveraging the probabilistic formalism of LES to design the algorithmic pipeline of our method. To the best of our knowledge, this has not been proposed before, therefore, we were not able to find other probabilistic LES  approaches to be used as learned baselines. \n\n**Writing**\n\nWe acknowledge your feedback on the clarity of the writing.\nWe will properly nuance the outlook section of the paper, particularly the connection of the results shown in this manuscript and the broader area of chaotic dynamics. We will also shift some of the explanations through the paper to better articulate how the motivation and modeling insights guided our choices in the algorithmic pipeline, in particular the choice of an NSDE in latent space to simulate DNS efficiently as a closure model.\n\n**References**\n\n1. Beck, A., & Kurz, M. (2021). A perspective on machine learning methods in turbulence modeling. GAMM‐Mitteilungen, 44(1), e202100002.\n1. Moser, R. D., Haering, S. W., & Yalla, G. R. (2021). Statistical properties of subgrid-scale turbulence models. Annual Review of Fluid Mechanics, 53, 255-286.\n1. Williamson, C. H. (1995). Vortex dynamics in the wake of a cylinder. In Fluid vortices (pp. 155-234). Dordrecht: Springer Netherlands.\nBosshard, C., Deville, M. O., Dehbi, A., & Leriche, E. (2015). \n1. UDNS or LES, that is the question. Open Journal of Fluid Dynamics, 5(04), 339.\n1. Chen, R. T., Rubanova, Y., Bettencourt, J., & Duvenaud, D. K. (2018). Neural ordinary differential equations. Advances in neural information processing systems, 31.\n"}, {"review_id": "7WA3wiKAAT", "reviewer": "Reviewer_GBgQ", "summary": "This submission proposed a data-driven method to learn a.closure model to simulate the results from DNS. The key part is a latent stochastic process by Neural SDE. And finally compute the Monte-Carlo approximation.", "questions": "1. Since the work treats DNS field as stochastic, why can you compute RMSE in the figure 3?\n2. Why do you design a stochastic process on latent space, instead of directly on original field? ", "limitations": "See questions and weakness.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The model treats the DNS as a stochastic process, instead of a deterministic process as in many previous works.\n\n2. Empirical results indicate it performs well on Kolmogorov flow.", "weaknesses": "1. The number of datasets is small.\n2. It lacks recent deep learning-based methods as baselines.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and feedback.\n\n**RMSE in Fig 3**\n\nYou are correct that we treat DNS fields stochastically. However, the ideal LES field is deterministic, which is approximated by a filtered DNS instance, and is a valid approximation for short time rollouts. \n\nOver longer time horizons, however, chaoticity dominates. In that regime RMSE ceases to be a meaningful metric; here we have provided statistical error metrics such as Turbulent Kinetic Energy (TKE) spectrum.\n \n**Baselines and number of datasets**\n\nWe have added cylinder wake as an additional dataset showcasing the performance of our method. In addition, we added another deterministic NN-based LES as baseline using the Neural-ODE framework. Please see our response to reviewer SnNJ. The conceptual framing of using SDE for probabilistic LES is the first kind to our best knowledge, so we do not have such probabilistic NN-based baselines to directly compare to. We hope our method becomes a baseline for future work.\n\n**Why not latent space on the original field?**\n\nThe main objective of LES modeling is the reduction of the computational cost. In general, the LES formulations tend to be more expensive at the same resolution than the original systems. However, the cost reduction is achieved by using a much coarser grid compared to the original one, which, compounded by the larger time-steps, results in an overall cheaper method. \n\nIn our case, running an SDE at the original resolution would be equivalent to running the DNS, which would ultimately defeat the purpose of LES. \n\nThe main insight/inductive bias we would like to incorporate is that: instead of running multiple DNS trajectories, we assume that there is a low-dimensional latent space where we can efficiently sample multiple short-term trajectories; and whose aggregated statistics remain close to the statistics stemming from directly computed DNS trajectories.\n"}, {"review_id": "T1Z0eo29FU", "reviewer": "Reviewer_a9bB", "summary": "This paper introduces a data-driven method for approximating the closure term in Large Eddy Simulation (LES). The closure term represents the unresolved scaling effect caused by reducing the computational grid through downsampling. In comparison to the traditional physically-informed approach, the learned closure term has the ability to automatically capture relevant features without the need for specific domain expertise or manual design. The model improves temporal resolution by employing latent space evolution with stochastic propagation. The proposed model demonstrates competitive performance compared to the Filtered DNS reference across different wave numbers, and significantly outperforms the baseline model.", "questions": "Q: It would be better if the authors could show the stochastic property in latent space brings diverse and valid results in the LES space.", "limitations": "N/A", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The proposed method outperforms several baselines in terms of the squared error and the accuracy of turbulent kinetic energy spectrum.\n- The latent space evolution implements the procedure to solve the SDE equation, which introduces the domain specific inductive bias into the model.", "weaknesses": "- The model's training solely relies on a single reconstruction loss, which is atypical for stochastic latent variables. The absence of prior regularization with respect to latent space implies that the stochastic component added at each temporal step acts more like a sparsity regularization, as the model is trained to minimize the stochastic effect rather than explicitly capturing the physical dynamics. Apart from prediction accuracy, there is a lack of additional benchmarks to demonstrate the advantages of the stochasticity, such as prediction diversity.\n\n- When comparing the proposed niLES with deterministic Neural Network-based models, it should be noted that niLES introduces additional parameters $h_{\\phi}$ and $g_{\\theta}$ as well as more computational routes, potentially leading to an unfair comparison. The authors should showcase the model's performance while varying the magnitude of the Wiener process in Equation 15 and compare it to the scenario where the noise magnitude is zero.\n\n- The current evaluation protocol does not adequately assess the scalability of the model, particularly with respect to wave numbers. To address this limitation, the authors should exclude a certain range of wave numbers as test cases, enabling a more comprehensive evaluation of the model's scalability.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review. Please find our response below.\n\n**Single reconstruction loss**\n\nWe apologize for the confusion. Eq. (19) indeed contains only the reconstruction loss, which is only a fraction of the training loss we ultimately use for training (see line 218 right after Eq(18)). Eq. (30) in Appendix D shows other loss terms related to the latent space. We will refactor this section and add these additional terms to the main text. \n\n**Prediction diversity**\n\nThis is a great suggestion. We have added a plot in the uploaded PDF (Fig 3) to show how the rollouts in the LES space vary due to the stochastic variation of the latent space trajectories.\n\n**Magnitude of the Wiener process**\n\nWe are adding a new baseline where the magnitude of the Wiener process is zero; see Figs 1,2,5 in the attached PDF. This corresponds to a deterministic NN whose latent space dynamics are parameterized as a neural ODE. We added additional numbers of layers (3x more layers) in the neural ODE drift function to compensate for the fact that the neural SDE contains additional parameters corresponding to prior drift and diffusion. The resulting neural ODE-based architecture has strictly more parameters than our proposed model niLES.\n\n**Scalability of the model**\n\nThank you for your feedback on increasing the comprehensiveness of the evaluation. However, we are not fully sure that we understood what you mean by “exclude a certain range of wave numbers”. Could you please clarify your comment?\n\nTo broaden our evaluation, we have added several more test cases for the evaluation of our model. In particular, for the Kolmogorov dataset we have included the RMSE and TKE spectrum for 6 independent examples. (See Figs 1 and 2 in the attached PDF.)\n\nAdditionally, we have included another scenario – cylinder wake at Reynolds number 500 (some more details are in our response to reviewer SnNJ). This scenario has more than 5x the degrees of freedom (52K vs 9.2K) in the DNS simulation compared to the Kolmogorov flow and exhibits irregular, chaotic vortex shedding. The LES simulation in the cylinder wake case contains 7.7K degrees of freedom compared to 2.3K in the Kolmogorov LES field.\n\nFor the cylinder wake dataset, we have included the average RMSE and the spread (min / max) as obtained by our method and the baseline methods. See Fig 5 in the attached PDF. We hope this expands the evaluation satisfactorily.\n"}, {"review_id": "CMDEM0up89", "reviewer": "Reviewer_hr9G", "summary": "The authors introduce a neural SDE model for LES flow fields. The model structure is motivated by the ideal LES approach and the model learns a data-driven closure term. The learned latent representation captures the variability and fine-scale structure of the fully-resolved DNS solution that is lost by the LES filter.", "questions": "1. Eq (2): Why can the pressure be eliminated?\n1. Line 87: Which boundary conditions do you impose and how do you do so implicitly?\n1. Is there a mix-up of $h_\\theta$, $h_\\psi$, $h_phi$, $g_\\theta$ and $g_\\phi$ happening in Eq (8), (9), (10)?\n1. Eq (14): You could define $f$ more explicitly here.\n1. Eq (16): What do you mean with spatio-temporal lifting? Does w correspond to z? What is $\\mathcal{D}$?\n1. What exactly is the relationship between the SDE solver timestep and the simulation timestep? Are the temporal intervals $\\mathcal{T}_i$ rescaled to $[0, 1]$?", "limitations": "The authors discuss limitations adequately.", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The paper is excellently motivated, presented and embedded into the theory. Stochastic modeling for turbulence as presented in this paper is an important direction to explore and therefore the paper would be valuable contribution to the field.", "weaknesses": "1. The authors seem to use the term \"closure model\" more loosely than usual in the introduction. Out of the citations for data-driven closure models in line 39, only [45] defines a closure model as the term is used in the context of RANS and LES.\n1. Lines 177-182 are slightly misleading as they can be understood to suggest that the model would learn possible DNS realizations as $Z_t$ instead of an opaque latent representation. (Actually, lines 211-213 clarify this but could be moved forward/merged into the first paragraph)\n1. A figure showing the flow (and predictions) after, for example, 200, 400, etc. steps would help the reader appreciate how chaotic the flow is and how much the velocity field evolves over 800 steps.\n1. A figure comparing the DNS and filtered data would also be helpful for the reader to understand the effect of LES filtering better.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the detailed and thoughtful review. We are glad that you found our paper well-motivated and that it provides an important direction to explore, viz. stochastic modeling of turbulence using neural networks.\n\n**Writing**\n\n>The authors seem to use the term \"closure model\" more loosely than usual in the introduction. Out of the citations for data-driven closure models in line 39, only [45] defines a closure model as the term is used in the context of RANS and LES.\n\nThat is a good catch! We’ve indeed used broader strokes in the introduction to give an overview of the extensive literature in closure models and other forms of data-driven approaches for turbulence. We will fine-tune the writing in the introduction to reflect the specific nature of closure modeling for LES and RANS that we undertake in this work.\n\n>Lines 177-182 are slightly misleading as they can be understood to suggest that the model would learn possible DNS realizations as Z_t instead of an opaque latent representation. (Actually, lines 211-213 clarify this but could be moved forward/merged into the first paragraph)\n\nIn line 177-182, we started the explanation by introducing one idea at a time, that of the need to model the DNS, and later on specify we just need a latent representation and not the whole DNS.  You are right that it can create confusion as to what our method is actually doing. Instead of saying that we “generate DNS trajectories” we can be more explicit at the beginning of the first paragraph that we learn a representation of said DNS trajectories. Furthermore we can clarify in these paragraphs that we indeed learn an opaque representation in latent space for these DNS trajectories. \n\n**Figure showing the flow evolution and the DNS snapshots with corresponding filtered versions**\n\nThank you for the excellent suggestions! We will include these figures in the camera-ready version of the paper.\n\n**Elimination of pressure**\n\nThe pressure can be eliminated in the Navier Stokes equations by e.g. taking the divergence of the momentum equation. We use the fact that the divergence of velocity is zero. It leads to a Poisson-like equation for the pressure, which can then be substituted into the momentum equation again. This results in an equation where the velocity is the only unknown. For more details, we provide the reference [1, page 295, section 6.2.2]. We will also add the reference to the corresponding sentence in the camera-ready version of the paper.\n\n\n**Boundary conditions**\n\nThe boundary conditions we impose in the Kolmogorov flow are 2D periodic along both axes. In the newly added cylinder wake dataset, we follow a setup of a static circular cylinder with no slip boundary condition on the cylinder surface, constant Dirichlet boundary condition on the inflow wall, periodic along the vertical axes, and homogeneous Neumann on the outflow wall. (Please see our response to reviewer SnNJ for more details.)\n\nIn line 87, by ‘implicit’ we meant that the boundary conditions are assumed to have been enforced in subsequent equations (Eq. 4, Eq. 5) without mentioning the fact explicitly. We can see how it may cause confusion with implicit numerical methods or linear solves. We apologize for the oversight. We shall rephrase this section to avoid confusion.\n\n**Notational issues in Eq (8), (9), (10) and (14)**\n\nThank you for catching the error. Eq (8) should say $g_\\theta$ instead of $g_\\phi$; and Eq (9) should say $h_\\phi$ instead of $h_\\psi$. We will fix these in the manuscript. Appendix D also reiterates the neural SDE equations with the corrected notation.\n\nThank you for the suggestion on Eq (14), in the text we can spell out $f$ in terms of the small scale fluctuations induced on the LES field from the DNS field.\n\n**Spatio-temporal lifting and SDE solver/simulation timesteps**\n\nWe refer to the DNS trajectory $w$ as being the spatio-temporally lifted version of the latent space trajectory $Z_t$, as the $w$ lives in both spatially higher resolution as well as temporally in the finer timescale. $\\mathcal{D}$ is the decoder (defined on line 203) which sends the latent space trajectory $Z_t$ to the fluctuation on the LES field caused by the DNS trajectory $w$. The output space of $\\mathcal{D}$ has the same dimensionality as the LES field.\n\nThe SDE timesteps are the same as the DNS simulation timesteps, which are an order of magnitude smaller than the LES simulation timesteps. They are not same as $\\mathcal{T}_i$ (or rather the timesteps we use for the filtered DNS training data) since the filtered DNS trajectories are downsampled by ~10x. As the SDE is meant to capture the fluctuations in the fast timescale, in each LES step the SDE evolves forward by ~10 timesteps. We do rescale the timesteps to $[0, 1]$ but that is for simplicity and is an implementation detail.\n\n\n\n**References**\n\n1. Deville, M. O., Fischer, P. F., & Mund, E. H. (2002). High-order methods for incompressible fluid flow (Vol. 9). Cambridge university press.\n"}]

Author rebuttal:
We thank all reviewers for providing thoughtful reviews and constructive feedback. We are encouraged by the positive comments that our method is well-motivated, clearly presented and provides an important direction to explore in the area of data-driven turbulence closure modeling. 

To address the reviewers’ concerns regarding the need for a broader evaluation, we have made the following high level changes:

- **Additional learned baseline.** We have added an additional deterministic NN-based baseline method: this architecture uses the same encoder-decoder but uses a Neural ODE based latent evolution. This corresponds to zeroing out the stochasticity due to the Wiener process in the Neural SDE formulation. Furthermore, to address the concern that the number of parameters in the Neural SDE might be higher because of the prior, posterior drift and diffusion, we have used 3x more layers (12 vs 4) in the drift function of the Neural ODE to compensate for this effect. We have reevaluated the test cases for Kolmogorov flow in Figs 1 and 2 of the attached PDF including the new baseline, and also included four additional testcases from the Kolmogorov flow dataset.
- **Additional turbulence scenario.** We have added an additional dataset – cylinder wake, which is a challenging instance of a chaotic Navier Stokes flow. In addition to its complex geometry, this flow has a high Reynolds number of 500, which is known to exhibit irregular vortex shedding. We have added the evaluation pipeline for our method and the baselines using this new dataset. Our response to reviewer SnNJ contains more details. See Fig 4 in the attached PDF for a qualitative plot and Fig 5 for accuracy plots.

Furthermore, we have addressed (or will address in a final version) the comments regarding typos and clarity issues pointed out by the reviewers. We hope that our explanations are satisfactory, and we are happy to answer any further questions that the reviewers may have.


Meta-review: This paper proposes a new methodology for designing closure models that utilizes a neural SDE incorporating the inductive bias of ideal LES. A Transformer based architecture is employed for the encoder-decoder mapping between the LES field and a latent space. The numerical experiments properly show the effectiveness of the proposed method niLES.  
The idea of the proposed algorithm is interesting and looks promising. The stochastic modeling would be a meaningful direction that should be explored by follow-up researches. In summary, I think this paper deserves publication in NeurIPS2023.
Decision: Accept (poster)