What symmetry buys a learned turbulence model

An interactive companion to my paper:

Comparison of Data-Driven Symmetry-Preserving Closure Models for Large-Eddy Simulation

Syver Døving Agdestein and Benjamin Sanderse · arXiv preprint · 2026

preprintcodeblog post

Everyone building neural closures for turbulence has heard the advice: bake the symmetries of the Navier–Stokes equations into your network, and it will be more accurate and generalize better. It sounds so reasonable that it is rarely tested. So we tested it — three networks, same inputs, same outputs, same training data, spanning "no symmetry constraints at all" to "exactly equivariant by construction". This post is the paper with the knobs exposed. If you only remember three things:

1. All three architectures — constrained or not — converge to the same accuracy. And that shared floor is not a mystery: it is the optimal one-point closure, which we can measure directly from data, no training required.
2. Symmetry constraints don't lower that floor. They buy you the trip to it with 25× fewer parameters, and make training reproducible where the unconstrained network is a lottery.
3. Generalizing across Reynolds number turned out to be less about architecture and more about one missing input: a filter-scale Reynolds number. Feed it to the network and the calibration problem largely disappears.

The problem, in one paragraph

Simulating every eddy in a turbulent flow (direct numerical simulation, DNS) is gloriously accurate and hopelessly expensive. Large-eddy simulation (LES) simulates a filtered — locally averaged — flow on an affordable grid, and pays for it with an unclosed term: the sub-filter stress $\tau (u)$, the effect of all the eddies the filter erased. A closure model $m(\overline{u})$ predicts that stress from the resolved field alone. Classical closures (Smagorinsky, Clark, …) are physics-derived; the modern alternative is to train a neural network on filtered DNS data. This paper is about what happens when you do that with — and without — symmetry constraints.

What "respecting a symmetry" means

The Navier–Stokes equations don't care how you orient your coordinate axes. Rotate the flow and the physics rotates with it — so the exact sub-filter stress rotates with it too. A closure model respects this if the two orders of operation commute: predicting then rotating gives the same answer as rotating then predicting. That property is called equivariance. Try it yourself on a toy closure (the mismatch is computed honestly from the two paths, not scripted):

Rotate the world by

Closure

Rotated velocity field

Predict, then rotate

Rotate, then predict (outline: other path)

Difference between the two paths: 0 (the two orders commute)

An unconstrained network is free to learn functions like the one behind the second toggle: perfectly fine on the training data, but reading raw tensor components off a fixed coordinate frame. Rotate the world, and its prediction doesn't rotate along — the model has invented a preferred direction that the physics doesn't have. (The toy exaggerates for visibility; our actual trained MLP violates equivariance by a subtler but very measurable ≈ 7%.)

The grid breaks the circle

Here is the wrinkle that shapes the whole paper: the continuous equations are symmetric under every rotation, but your simulation doesn't live in the continuum — it lives on a grid. Drag the slider and watch which rotations map the grid onto itself:

Rotation angle θ = 35°

✗ 99% of the rotated points miss the grid — a grid-sampled field cannot realize this rotation exactly.

Only multiples of 90°, plus mirror flips, survive. In 3D these grid-compatible roto-reflections form the octahedral group: exactly 48 elements out of the infinitely many rotations of the continuous world. Any closure that acts on grid-sampled data can be exactly equivariant to those 48 — and to no more. That number, 48, will come back below as the size of a network layer.

The same bookkeeping applies to the other Navier–Stokes symmetries, and it tells us what a closure can and cannot promise:

Survive the grid ✓

Translations, time shifts, constant-velocity (Galilean) frame changes, and the 48 octahedral roto-reflections. These are exact symmetries of the discrete equations — enforcing them in the closure costs nothing and can only help.

Survive partially ~

Scaling: no fixed grid maps to itself under rescaling, but the closure can still transform correctly across rescaled systems. Our input/output construction handles the part that doesn't involve viscosity; the rest returns in the Reynolds section.

Broken ✗

Continuous rotations (reduced to the 48), and accelerating frames (broken by the fixed forcing that keeps the turbulence alive). A closure can still enforce full rotational symmetry as a continuum ideal — the question is whether it should.

Three networks, one contract

All three data-driven closures in the paper sign the same contract: input is the local, normalized velocity-gradient tensor $\overline{A}/|\overline{A}|$, output is the predicted stress times the prefactor ${\mathrm{\Delta }}^2|\overline{A}{|}^2$ (with filter width $\mathrm{\Delta }$). Using gradients instead of velocities makes every model Galilean-invariant for free, and the prefactor gives the output the units and scaling behavior of an actual stress. The three differ only in how they treat rotations and reflections:

MLP — enforce nothing

A plain multi-layer perceptron mapping the 9 gradient components to the 6 stress components. No symmetry constraints; whatever equivariance it has, it learned from data. The control group.

G-CNN — constrain the weights

A group-convolutional network: each hidden layer carries 48 copies of its channels, one per octahedral element, with weights tied so that rotating the input provably permutes the copies. Exactly equivariant in floating point, by construction.

TBNN — constrain the features

A tensor-basis network: the stress is expanded in 7 tensors built from the strain and rotation rate that rotate correctly no matter what, and the network only predicts their 7 scalar coefficients from 5 rotation-invariant inputs. The network itself can be anything.

The G-CNN and the TBNN are the two classic routes to equivariance — constrain the architecture or constrain the representation — and part of the fun of the paper is that they arrive at the same place from opposite directions.

Same floor, very different price

Now the main event. Train all three at a range of network sizes and plot the error against the parameter count:

MLP (unconstrained)G-CNN (equivariant)TBNN (tensor basis)optimal closure (measured)

Relative error of the predicted sub-filter stress against the exact one, on held-out snapshots. The dashed black line is the training-free conditional-mean estimate of the best possible one-point closure. Whiskers span ± one standard deviation over five training seeds. Hover or tab across the columns for exact values.

Data table

Parameters (≈)	MLP (unconstrained)	G-CNN (equivariant)	TBNN (tensor basis)
120	0.896 ± 0.014	0.452 ± 0.006	0.445 ± 0.001
400	0.675 ± 0.019	0.449 ± 0.001	0.445 ± 0.001
1200	0.460 ± 0.001	0.449 ± 0.001	0.446 ± 0.001
3000	0.449 ± 0.000	0.449 ± 0.001	0.446 ± 0.000
no model	1.000
dyn. Smagorinsky	0.964
Clark	0.466

Two things to see here. First, the flat green and blue lines: the TBNN and G-CNN are already saturated at ~120 parameters. Their built-in structure leaves so little to learn that tiny networks suffice. The unconstrained MLP starts out worse than the 50-year-old Clark model and needs ~3000 parameters — 25 times more — to catch up. (It's also a gamble at small sizes: its seed-to-seed spread is an order of magnitude wider than the constrained models'.) Second, and to me the more interesting finding: nobody gets below the dashed line. All three architectures, however different their insides, converge to the same error and stop.

That dashed line is not a fit to the networks — it is measured independently. The best possible predictor of the stress from a single point's velocity gradient is a conditional average: over all turbulent flows that share this local gradient, what is the mean stress? That is the one-point optimal closure in the sense of Langford & Moser, and because it is an average, you can estimate it by sorting the training data into bins and averaging — a histogram, no training, no architecture. It lands at an error of 0.439, and the saturated networks land within 0.6–2.3% of it. The shared floor is the optimal closure. The remaining 44% error is not a failure of the networks: filtering destroys information, and no function of a single-point gradient — neural or otherwise — can recover it.

Two entertaining footnotes from this measurement. The binned estimate barely improves as the bins get finer (0.442 with one global bin — i.e. seven constant coefficients — down to 0.439 with 65 536 bins): almost everything a one-point closure can know is already in seven numbers. And averaging the trained MLP's predictions over the 48 octahedral elements makes it exactly equivariant at no accuracy cost — its 7% symmetry violation was wasted degrees of freedom, not useful ones.

Does it work when you actually run it?

A-priori stress errors are a laboratory metric. The real test is a posteriori: plug each closure into an actual LES and compare against the filtered DNS. The a posteriori toggle in the chart above shows the learned closures essentially tied at a solution error of ≈ 0.35, ahead of Clark (0.36), dynamic Smagorinsky (0.39), and no model at all (0.43). But the aggregate number hides the most physical difference — what the closures do with energy:

exact stress (reference)dyn. SmagorinskyClarkTBNN (learned)

Distribution of the pointwise sub-filter dissipation rate over all test snapshots (kernel density estimate, log scale), every closure evaluated on the same filtered DNS fields. Negative values mean energy flowing back from small to large scales. Hover for exact densities.

The exact sub-filter stress mostly drains energy from the resolved scales (the big positive lobe), but at 22% of all grid points it pushes energy backwards, from the unresolved eddies into the resolved flow — backscatter, the shaded region. Eddy-viscosity models like Smagorinsky cannot do this even in principle: their density stops dead at zero. Clark captures the backscatter tail (slightly too enthusiastically) but misses the heavy forward-transfer tail. The learned closure reproduces both tails — underestimating backscatter (it covers about half the reference fraction), which is the fingerprint of it having learned a conditional average: means are milder than samples.

One number to rule the Reynolds number

Here is the quiet limitation of the contract all three networks signed: the input $\overline{A}/|\overline{A}|$ and prefactor ${\mathrm{\Delta }}^2|\overline{A}{|}^2$ contain no viscosity. The closures are Reynolds-blind — whatever dissipation calibration they learned is pinned to the training regime. The symmetry analysis even predicts the fix: the full scaling symmetry of Navier–Stokes involves the viscosity, and the only scaling-invariant scalar you can build from the available quantities is a filter-scale Reynolds number,

$${\mathrm{Re}}_{\mathrm{\Delta }}=\frac{{\mathrm{\Delta }}^2\sqrt{⟨|\overline{A}{|}^2⟩}}{\nu }.$$

So we trained each closure across several viscosities and filter widths, with and without ${\mathrm{Re}}_{\mathrm{\Delta }}$ as one extra input, and tested on held-out Reynolds numbers reaching well outside the training range. Flip the toggle:

MLPG-CNNTBNNtraining range

feed the closures Re_Δ

Median sub-filter dissipation divided by the true value (1 = perfectly calibrated, log scale). Reynolds-blind closures slide from over- to under-dissipation as the test Reynolds number grows; with the Reynolds input they hold their calibration. Faint lines show the +Re counterpart. The shaded band is the range of Re_Δ seen during training; everything to its right is extrapolation. Hover for exact values.

With the toggle off, the story is a slide: over-dissipative at low ${\mathrm{Re}}_{\mathrm{\Delta }}$, under-dissipative at high — drifting from ~1.35× down to ~0.75× the true dissipation as the flow leaves the training band. Flip it on and the calibration holds out into extrapolation. Meanwhile the stress error tab barely moves: the tensor structure of the sub-filter stress is nearly Reynolds-independent; what changes with Reynolds number is only how much energy lives below the filter. One scalar input carries exactly that information — and notice that the three architectures never separate: the gain comes from the input feature, not from equivariance.

The honest caveat: on a completely different flow (a decaying Taylor–Green vortex — laminar roll-up, transition, decay, nothing like the forced training soup), the closures over-dissipate by ~2× and the Reynolds input recovers only part of it (2.07× → 1.80× at the tightest filter). A single global scalar can tell the closure which Reynolds regime it is in, but not which kind of flow it is in. The models stay stable and keep the lowest errors through the transition, though — while Clark diverges outright at the wider filters.

What I take away

• Test the premise before building the cathedral. Equivariance is routinely motivated by accuracy and generalization. In this controlled comparison it delivered neither — because the unconstrained optimum was already (statistically) equivariant, even a plain MLP could reach it with enough parameters and data.
• What symmetry actually buys is efficiency and reliability. The floor at 25× fewer parameters, reproducibility across seeds, exact physical consistency in floating point. Those are real currencies — cheaper training, cheaper inference, and no spurious frame-dependent forces — just different ones than advertised.
• Know your floor. The conditional-mean estimate costs nothing, requires no training, and tells you the best any pointwise closure can do before you design one. If your fancy architecture is at the floor, more capacity is wasted; the only way down is richer inputs — neighborhoods, history, or scalars like ${\mathrm{Re}}_{\mathrm{\Delta }}$ that carry genuinely new information (the +Re TBNN measurably dips below the one-point floor out of distribution).
• Generalization can be a feature-engineering problem. A decade of transfer-learning machinery for crossing Reynolds numbers, and here most of the job is done by handing the network one well-chosen number.

---

The interactive charts use the exact evaluation data of the paper — forced isotropic turbulence at ${810}^3$ DNS resolution, closures trained across three viscosities and three filter widths, five training seeds each — extracted from the archived result files of the companion code. The paper additionally reports the equivariance errors, timing, energy spectra, and the Taylor–Green study.