Woudschoten Conference · 2 October 2026

What a coarse grid cannot know

Learned turbulence closures, and how AI agents help us find out

Syver Døving Agdestein · CWI Amsterdam · with Benjamin Sanderse

Supported by NWO (VI.Vidi.193.105, OCENW.XS25.1.149)
Compute on Snellius: SURF (EINF-15798)

Good afternoon. [The programme still lists Benjamin Sanderse with "Everything flows"; say in one line that you are giving this slot instead.] I am Syver, a PhD researcher with Benjamin at CWI. On the left of the background you see a flow in full detail; on the right, what a coarse grid sees of it: one number per cell. This talk is about that dashed line. What can a coarse grid know about the flow it does not resolve, and what can it not know? Three parts on the science, and at the end something different: how AI agents helped us find out. (Timing: 26 numbered slides for about 30 minutes; slide 18 can be skipped. Funding: NWO Vidi "Discretize first, reduce next" and NWO XS "Discovering neural stochastic differential equations to simulate probabilistic turbulence"; computing time from SURF.)

Taylor–Green vortex at Re = 3000 breaking down into turbulence, shown as vortex tubes

The question

How much turbulence can a coarse grid know about?

A simulation on a coarse grid sees the flow only through that grid. What does that view determine, what can it not determine, and what can we do about the rest?

Taylor–Green vortex at Re = 3000 becoming turbulent (vortex structures)

I · Coarse grids

Turbulence has too many scales

Reynolds number Re: how turbulent a flow is

Number of grid points proportional to Re to the power 9/4

grid points to resolve every eddy: direct numerical simulation (DNS)

Large-eddy simulation (LES): a coarse grid, plus a model for what lies beyond the cut.

Grid-point estimate: Pope, Turbulent Flows (2000), ch. 9

I · Coarse grids

Coarse equations are not closed

Velocity slice on the fine grid, full of small eddies

The same velocity slice filtered to a coarse grid

DNS grid, 810³: every eddy

LES grid, 128³: the filtered flow ū

The subfilter stress τ needs the small scales, which the coarse grid does not have.

A closure model must supply τ from ū alone.

How does LES work? Filter the velocity: average it over the size of a grid cell. That gives u-bar, the flow as the coarse grid sees it. Left, a velocity slice on the DNS grid, 810 cubed points; right, the same flow filtered to a 128 cubed LES grid. Filter the Navier–Stokes equations and everything can be written in terms of u-bar, except one term: tau, the subfilter stress, the difference between the filtered product and the product of the filtered velocities. [click] Tau depends on exactly the small scales the coarse grid does not have. [click] So a closure model has to supply it from u-bar alone. That is the closure problem. (Bridge: Roel Verstappen's talk this morning was about this same problem. Our target tau is always the exact stress of the discrete coarse scheme; that detail is on a hidden slide if asked.)

II · Learned closures

Learn the closure from simulation data

input: the resolved velocity gradient Ā, 9 components

output: the stress τ, symmetric: 6 independent

a small network, applied at every grid point

Training data: filtered DNS of forced isotropic turbulence, 3 viscosities × 3 filter widths

So we learn the closure from data. We run fully resolved simulations, filter them, and compute the exact stress. Left: the nine components of the resolved velocity gradient Ā on a slice of the coarse grid. In the middle, a small neural network that is applied independently at every grid point: nine numbers in. [click] Six numbers out: the stress is a symmetric tensor, so the three greyed components below the diagonal are copies of the ones above. [click] In a formula: tau at a point is approximated by a model M-theta of the gradient at that point. Using only the local gradient as input is a modelling choice, a locality choice; classical models such as Smagorinsky make the same choice. The training data cover three viscosities and three filter widths. Which network should we use, and what does it need to know?

II · Learned closures

Symmetry: one flow, many views

Engraving of Archimedes in a wooden bathtub

Archimedes in his bathtub on a horse-drawn cart moving at constant speed

Archimedes in the bathtub on a rotating platform

the flow in the tub

move it at constant speed: nothing changes

rotate it: the flow rotates along

mirror it: the flow is mirrored

Build it in: share weights over the 48 (G-CNN), or use invariants and a tensor basis (TBNN).

A cube grid keeps 48 rotations and reflections. A closure should respect all of them.

The equations of fluid motion do not care how you look at them. Put Archimedes' tub on a cart moving at constant speed [click]: nothing inside changes. That is Galilean invariance, and it is why a closure may not depend on the velocity itself. Rotate the tub [click], and the flow inside rotates with it. Mirror it [click], and the flow is mirrored. [click] A cubic grid cannot keep every rotation, but it keeps 48: the rotations and reflections of the cube. A closure should be equivariant: rotate the input gradient and the predicted stress rotates the same way. Two standard ways to build this in: a group-equivariant network that shares its weights over the 48 transformations, or a tensor-basis network that works with invariants and a fixed tensor basis. Classical closures like Smagorinsky are equivariant by construction; a generic neural network is not.

II · Learned closures

Two ways to build the symmetry in

G-CNN · share the weights

TBNN · invariants in, a tensor basis out

zeta of g equals phi of the sum over h in O_h of the shared kernel k of h inverse g times xi of h, plus b

M of Abar equals Delta squared times norm of Abar squared times the sum over k from 1 to 7 of alpha k of lambda 1 to lambda 5, times T k of S and W

48 copies of every neuron share one kernel k. Turning the input only moves the copies around, so the output turns with it.

The network sees only five invariants, which no turn changes. Seven fixed tensors built from S and W carry the turn.

Random weights, exact symmetry · G-CNN: Cohen & Welling 2016 · TBNN: Ling et al. 2016, basis of Stallcup et al. 2022

Two ways to build the 48 symmetries in. Left: a group-convolutional network. Every neuron is copied 48 times, one copy per rotation or reflection of the cube, and all copies share one kernel: the response of copy g sums the copies h with the shared kernel evaluated at h-inverse g. Watch what happens when the input gradient is turned: its nine entries move and change sign, and the 48 copies only swap places, so the output turns exactly as the input did. It is the idea of an ordinary convolutional network, where sharing weights over shifts gives translation equivariance, with the shifts replaced by the 48 cube symmetries. [click] Right: a tensor-basis network. Instead of constraining the weights, it constrains what goes in and what comes out. In: five invariants of the normalized gradient, such as the trace of S squared, numbers no rotation changes. Out: seven coefficients that weight seven fixed tensors built from the strain S and the rotation W. Turn the input, and the invariants and coefficients do not move; the basis tensors carry the turn. The tensor-basis form goes back to Pope in 1975 and to Ling and colleagues for RANS; we use the minimal basis of Stallcup and colleagues. The animations use random weights: what you see is exact equivariance, not a trained model. Both are exactly equivariant to the 48; the tensor basis even to every rotation, as a pointwise map. Next: what that buys.

II · Learned closures

Every architecture hits the same floor

25×

fewer parameters when the symmetry is built in

…but no network gets below the floor.

The floor: the error of the best possible function of Ā, set by the input, not the architecture.

Median over 5 seeds; floor from a nonparametric estimate · Agdestein & Sanderse, arXiv:2603.05325

Here is what symmetry buys, on a flow the networks never saw. Error against network size. The plain network, without built-in symmetry, needs about three thousand parameters to get good. [click] The two equivariant networks are there with about 120. [click] That is 25 times fewer parameters. [click] But look at the dashed orange line. All three networks saturate at the same level, and making them bigger does not help. That floor is a nonparametric estimate of the error of the best possible function of the local gradient: the conditional mean. So the remaining error is not a failure of the architecture; it is a property of the input, something the local gradient cannot know. (Scope: this floor was scored on projected predictions, so it is not an exact local conditional variance.) What is missing from the input?

II · Learned closures

Blind to the Reynolds number

Archimedes in a large bathtub next to a small copy of himself in a small bathtub

x goes to a x, t goes to b t, implies nu goes to a squared over b times nu

A closure that sees only Ā and Δ must, by dimensional analysis, take the form

so it cannot tell apart flows that differ only in viscosity.

add the one dimensionless number the symmetry allows:

Re_Delta equals Delta squared times norm of Abar over nu

Symmetry hands the closure its missing input.

So what is missing? Another symmetry, one that is usually ignored: scaling. Shrink Archimedes and his tub: the flow looks the same only if the viscosity is rescaled too. Stretch space by a and time by b, and the viscosity must go to a squared over b times ν. [click] Now take a closure whose only inputs are the local gradient and the filter width. Dimensional analysis forces it into this form: a magnitude times a function of the gradient's direction. It cannot tell apart two flows that differ only in viscosity. [click] The scaling symmetry allows one more dimensionless input: the filter-scale Reynolds number, Re_Δ, the filter width squared times the gradient magnitude of the whole flow, divided by the viscosity. [click] So symmetry does not just constrain the closure; here it hands the closure its missing input.

II · Learned closures

Now the dissipation stays calibrated

Reynolds-blind: 1.36 → 0.75

drifts from too much to too little

With Re_Δ: 1.08–1.27

stays calibrated, also beyond the training range

Six unseen test flows, median over 5 seeds · Agdestein & Sanderse, arXiv:2603.05325

$Three turbulence simulations with identical large scales whose filtered fields drift apart$

III · Beyond the mean

One coarse state, many futures

Three simulations of 3D turbulence: identical large scales, tiny differences in the small eddies.

Same resolved initial state. Different futures.

Taylor–Green vortex, Re = 1600, 256³ grid. Top: vorticity. Bottom: filtered to what a coarse model sees.

Part three. So far the closure was a function: one input, one output. Should it be? This film says no. Three simulations of the Taylor–Green vortex at Re = 1600 on a 256-cubed grid. Same deterministic equations, exactly the same large-scale flow at the start. Only the small eddies differ: a perturbation of 3 percent of the velocity amplitude, less than 0.1 percent of the energy, at wavenumbers 16 to 24. The bottom row is the coarse view, filtered at wavenumber 12, which cannot see that perturbation at all. Watch the bottom row: the coarse fields start identical and drift apart. Same resolved initial state, different futures. The coarse state does not determine the closure term; it determines a distribution of them. (The film is 20 seconds, silent, and loops. These three runs illustrate possibilities, not a calibrated distribution.)

III · Beyond the mean

The best closure is an average

Orange: the law of τ at the dashed line, and draws from it.

One slice of filtered DNS, 900 of 16,384 points · Ideal LES as a conditional average: Adrian (1990)

Now the newer part: what the floor is made of. Every dot is one grid point from a real filtered simulation. Horizontally, the shape of the local gradient, from rotation-dominated to strain-dominated. Vertically, the energy the subfilter scales drain from the resolved flow; below zero is backscatter, energy flowing back up. [click] The blue curve is the conditional mean: the best any closure can do in mean square, given this input. In the formula, M-star is the best of all closures M. That is what all the networks converged to. But the truth is the cloud. [click] The dashed line is one value of the input; on the right is the law of the stress there, and the orange dots are draws from it. The idea of an ideal closure as a conditional average goes back to Adrian in 1990, conditioning on the whole resolved field. Conditioning on the local gradient only, as here, is a locality choice.

III · Beyond the mean

So measure the whole law

A mean closure keeps the blue term and drops the orange one: all the spread.

The object

the law of the stress, given the resolved gradient

Made small by symmetry

5 → 5

five invariants of Ā in, five stress components out

Measured on new data

13 runs

new simulations up to 1200³, scored on held-out data

Next: how much is left over, what shape it has, and whether it survives in a running LES.

So we stop at the mean no more, and measure the whole conditional law. Why bother? The expected squared stress splits into the squared best closure plus the conditional variance. A mean closure keeps the blue term and drops the orange one: all of the spread, including the backscatter events. The object is the law of the stress given the local gradient. Symmetry makes it small: the traceless gradient has eight components, rotation removes three, so five invariants in; the deviatoric stress has five components out. A five-to-five law is something we can actually estimate. And we measured it on new data: thirteen new simulations on the Dutch national supercomputer, up to 1200 cubed, Taylor-scale Reynolds numbers 190 to 455, scored by likelihood on data held out from fitting, against a nonparametric reference. Three questions for the next slides.

III · Beyond the mean

Five numbers in, five out

8 − 3 = 5

A traceless gradient has eight numbers. Three of them only say how it is turned: look along the strain's own axes and five remain.

Write τ in the same axes: its five components are the output.

So the law is a 5 → 5 conditional density.

Each axis is fixed only up to its sign: a vorticity convention picks one, and Ā and τ flip together.

Chart: Meneveau 2011 · used for subgrid-stress models by Prakash, Jansen & Evans 2022

How does symmetry make the law small? Take the gradient at one point. Its symmetric part, the strain, is drawn as an ellipsoid with three principal axes; its antisymmetric part is the green vorticity arrow. In the lab axes, the nine entries of Ā, eight independent since the trace is zero, all change as the same local flow is turned. Now turn the axes onto the principal directions of the strain. The strain becomes diagonal: three rates that sum to zero, so two numbers. The vorticity in those axes: three more. Eight minus three rotations: five numbers in. [click] Write the stress in the same axes, and its five deviatoric components are the output: a five-to-five conditional density. The chart is in Meneveau's 2011 review of the velocity gradient, and Prakash, Jansen and Evans used it for subgrid-stress models; what we add is the conditional density on it. Two caveats. Each axis is fixed only up to its sign; a vorticity sign convention fixes it, and a change of convention flips Ā and τ together. And our cubic grid keeps only the 48 cube symmetries, so dropping the orientation is an approximation, which we checked with a spherical-mask control. Now the three questions.

III · Beyond the mean

The mean explains most of the stress

That remainder is what a stochastic closure is for.

Best closure given the local gradient · one viscosity, filter widths of 2 to 4 grid spacings

III · Beyond the mean

What remains is heavy-tailed

63–130

kurtosis, against 35 for a Gaussian

3.2–6.2

tail index of the Student-t fit; a Gaussian's is infinite

The long tail sits on the forward-transfer side: extreme stresses drain energy to the small scales.

Plots: one operating point · numbers: all eleven operating points

Second question: what shape is the part the mean leaves out? Left: the residual in two of its five coordinates, whitened; the white lines are contours of the measured density at 10^-2, 10^-3 and 10^-4. Right: the distribution along one coordinate, on a log scale. [click] A Gaussian with the same spread, dashed blue, falls off a cliff: five standard deviations out the turbulence is several hundred times more likely, and ten out, where the Gaussian says never, it is still clearly there. The kurtosis is 63 to 130 against 35 for a Gaussian. [click] A Student-t, dotted orange, with tail index between 3 and 6 follows the data over six decades. The contours are not circles, though: they are triangular. That shape appears in random-phase fields of the same spectrum too, so it is geometry rather than turbulence. [click] The physics is in the skew: the long tail sits on the forward-transfer side, so the extreme stresses are the ones that drain energy to the small scales, as flux measurements in the same flows find, and a skew-t captures it. (Scope: the skew statement is for the standard representation of the stress.)

III · Beyond the mean

…but two numbers carry most of it

92–93%

of what the law gains over “mean plus constant Gaussian noise” comes from a scale σ(Ā) and one tail index.

What the mean leaves out is low-order.

Share of the likelihood gain, terms added in the order shown, at each of the eleven operating points

So the residual is heavy-tailed and far from Gaussian. Does that mean we need a complicated generative model? Start from the mean plus constant Gaussian noise, add one ingredient at a time, and measure how much likelihood each buys, as a share of the total. The bar cycles through our eleven operating points: different viscosities and filter widths. Blue is a noise scale that depends on the gradient; orange is one global tail index. How the gain splits between the two changes from point to point: at wide filters the tail matters more. But together they always carry 92 to 93 percent. Everything else, the full covariance shape, the skew, a tail index per bin, adds little. [click] So what the mean leaves out is heavy-tailed but low-order. That is good news for anyone who wants to build a stochastic closure. (The individual shares depend on the order in which terms are added; the combined share does not. Earlier projected residual.)

III · Beyond the mean

Learn the law where it is used

A law learned on DNS data loses its whole advantage over mean + Gaussian noise within three snapshots of a running LES.

Relearned on the LES's own fields, its advantage even grows: LES errors make the residual more heavy-tailed.

Above 100% means a bigger lead over the Gaussian reference, not a better prediction: absolute accuracy still falls.

Quantifies a known effect: Zandonade, Langford & Moser (2004); Parthipan et al. (2023) · LES with a trained CNN closure

[Optional slide; skip if short on time.] Third question: does the law survive inside a running simulation? We start a large-eddy simulation from filtered DNS and let it run. The vertical axis is the law's likelihood advantage over the simplest stochastic model, the mean plus Gaussian noise, as a percentage of the advantage it had on DNS data. The law we learned on DNS data starts at 100 percent, and after three snapshots its advantage is gone, while the resolved field is only about 7 percent off. The mean itself is much more robust: it still keeps 76 to 96 percent of its skill at that error. That a-priori quality need not carry over a posteriori is a known effect: Zandonade, Langford and Moser saw it in 2004, and Parthipan and colleagues in 2023 on a likelihood score; what is new here is the number. [click] The good news: relearn the law on the LES's own fields and its advantage grows to about 140 percent, 131 to 149 percent across the closures we tried. That is not a better prediction than on DNS: the absolute likelihood still falls as the LES drifts. It means the Gaussian reference falls faster, because LES errors make the residual more heteroscedastic and heavier-tailed (in one run the refit tail index drops from 4.6 to about 3.7 by lead 8), which the law captures and a Gaussian cannot. So the law exists on the solver's own field; we should learn it there.

III · Beyond the mean

A one-point law is not yet a closure

Variance of the sum over T over Delta t steps of Delta t times xi_n equals T Delta t sigma squared, which goes to zero

Fresh noise at every step averages away as Δt → 0.

A stochastic closure without memory is deterministic in disguise.

From the one-point pdf f, “a time or length scale of turbulence cannot be deduced from f”: Pope (1983), p. 3448

III · Beyond the mean

A one-point law is not yet a closure

Noise with memory keeps its spread, whatever Δt.

The measured residual has memory: it decorrelates in about one strain time. The one-point law does not fix it.

Next: stochastic LES on Snellius. The first stage has run; no verdict yet.

From the one-point pdf f, “a time or length scale of turbulence cannot be deduced from f”: Pope (1983), p. 3448

IV · How this was done

How was all this done?

Of the 442 commits in this project since July, 346 list an AI agent as co-author.

IV · How this was done