Your grid is a filter

An interactive companion to our paper:

The subgrid stress

might be the most-modeled object in computational physics. Sixty years of closure models — Smagorinsky, dynamic, scale-similarity, and lately a wave of neural networks — all aim at it. The paper’s claim, stated bluntly: for any simulation you can actually afford to run, this is the wrong target. There is an exact expression for what a coarse simulation leaves unresolved, we can write it down, and it looks noticeably different — it knows your numerical scheme, it is not symmetric, and it is not even local. This post is the tour of how that happens, with the paper’s wall of equations replaced by things you can click.

The finite volume method never made an error

Everything starts from an identity that deserves to be more famous. Take the top-hat filter of width $h$, $\bar{u}^h(x) = \frac{1}{h} \int_{x-h/2}^{x+h/2} u(y) \, \mathrm{d}y$, and a finite difference over the same $h$:

This is not an approximation — both sides are the same number, exactly (it is the fundamental theorem of calculus wearing a hat). A coarse finite difference of the true field equals the true derivative of a filtered field. I call this the filter-swap property, and it means the grid is not where accuracy goes to die: the grid is a filter, and it tells you exactly which smoothed field your discrete operators are talking about.

Push a conservation law $\partial_t u + \partial_x r(u) = 0$ through this identity and you get

which is the finite volume method — flux difference across a cell, no approximation anywhere. The FVM is derived by filtering, not by truncating Taylor series. The only sin committed so far is that the flux $r(u)$ still involves the unfiltered field $u$, which the coarse grid does not carry. All the error in a coarse simulation enters in one single place: when you replace that unclosed flux by something computable. Which raises the question of when in the pipeline you commit that sin.

Two routes to a runnable simulation

Classical LES filters the continuous equations, models the subgrid stress while everything is still continuous, and only then discretizes — a second approximation, applied on top of a closure that never knew a grid existed. The paper’s route applies two filters instead: the LES filter $f^\Delta$ picks the physics you want to keep, and the grid filter $f^h$ makes the equation discrete — both steps exact. Click around:

The routes end at equations of identical shape, but they are not the same equation. On the discretize-first route, the single approximation step has a well-defined, exact target: the residual flux

Look closely at the second term: it contains $r^h$, the numerical flux — your actual scheme, interpolations, central differences and all. The discretization error is not floating around outside the closure problem; it is a term inside the stress, where a closure model can finally see it. If your scheme is dissipative (upwind), $\tau^{\Delta,h}$ automatically asks the closure for less dissipation; if it is a crisp central scheme, it asks for more. Implicit and explicit LES stop being different philosophies and become different values in the same expression.

Does it hold up? The paper runs the brutal test: a 13 500-cell Burgers DNS computes $\tau^{\Delta,h}(u)$ at every time step and injects it into a 300-cell coarse simulation, over 1000 random initial conditions. With the classical stress as the injected “closure”, the coarse run drifts to 7% error (at $\Delta = 4h$; 16% in the implicit-LES limit). With $\tau^{\Delta,h}$: error at machine precision, every grid, every filter width. Not small — zero. It is the closure problem’s rarest commodity: a ground truth.

How much of your “turbulence” is actually your grid?

Splitting $\tau^{\Delta,h}$ into the classical part plus the two discretization parts (numerical flux and discrete divergence) tells you how big the blind spot of a classical closure is. The answer depends entirely on the ratio of filter width to grid spacing:

At $\Delta = 32h$ the classical theory is vindicated — and nobody runs there, because it wastes a factor 32 of resolution. Where simulations actually live, between $\Delta = 0$ and $4h$, a quarter to nearly three quarters of the thing you are modeling is your own discretization. The classical framework handles this with folklore (“let the filter width absorb the grid effects”); the exact expression just accounts for it.

Then 3D happens: pressure is geometry

Burgers is a warm-up; the incompressible Navier–Stokes equations bring a constraint: $\nabla \cdot u = 0$. The pressure is not really physics here — it is a Lagrange multiplier, and the cleanest way to see it is geometric. Think of every velocity field as a point in a huge vector space. The divergence-free fields form a subspace — in this cartoon, a line. The pressure gradient is precisely the correction that moves your field onto that subspace, orthogonally. Drag the field around:

The paper uses this projection $\pi$ (and a cousin that acts on stress tensors) to eliminate the pressure entirely and write Navier–Stokes as a self-contained conservation law for the velocity — which is what lets the whole filter-swap machinery from 1D go through. But the projection charges a price: computing it means solving a Poisson equation, so it couples every point in the domain to every other point. And on the coarse grid, the filtered field is not discretely divergence-free, so it must be re-projected with the discrete projector, which is baked into the exact stress. The consequence is worth italics: the exact unresolved stress is non-local in the velocity field. A closure model that looks only at its own grid neighborhood cannot represent it, no matter how many parameters it has.

The exact stress is not even symmetric

The second 3D surprise comes from the filter-swap identity itself. In 3D the grid filter splits into face averages: swapping $\partial_j$ for the discrete difference $\partial^h_j$ forces the stress component $\xi_{ij}$ to be averaged over the cell face orthogonal to direction $j$:

Same underlying stress, averaged over two different surfaces: $\xi_{12} \neq \xi_{21}$. The exact unresolved stress in a finite-volume LES is non-symmetric — it has nine independent components, not six. Meanwhile, essentially every structural closure model in existence is built symmetric, because the continuous theory said so. Symmetric models are not approximating the exact target imperfectly; they are confined to a subspace that provably does not contain it.

The ladder to machine precision

Talk is cheap; the 3D experiment is not. A $500^3$ decaying-turbulence DNS runs alongside a $100^3$ coarse simulation, and the DNS computes the “closure” term injected into the coarse run — using five candidate stress expressions, from the textbook one to the paper’s exact one. Each rung of the ladder adds one ingredient from this post:

Two readings of this chart. Pessimistic: even with full DNS data in hand, the textbook stress expression cannot steer a coarse simulation better than 21% error — the target itself is off. Optimistic: the gap between 2.86% and zero is the non-symmetric part, the gap between 6.42% and 2.86% is the filter-swap face-averaging, and the jump from 21.2% to 6.42% is simply admitting your numerical scheme into the stress. Every ingredient is identifiable, and the ladder bottoms out at exactly zero.

Fine, but I can’t run a DNS next to my LES

True — DNS-aided LES only certifies the target; it is not a model. The practical payoff is what happens when you fit an actual closure to the right target. The paper redoes the humble Smagorinsky model, fitting its one coefficient by least squares against either the classical stress or the exact one. Fitted to the exact stress, the coefficient comes out consistently larger — it inherits the extra dissipation that the discretization demands — and the resulting simulations have lower errors and cleaner spectra, most visibly in the implicit-LES regime where discretization effects dominate. That is the cheapest possible model gaining accuracy purely from better target data. For data-driven closures with thousands of parameters, trained directly on residual data, the choice of target is the whole game — that was the lesson of our earlier JCP paper, and the present paper supplies the exact target it called for.

What I want you to remember

• The grid is a filter, and the FVM is exact. Discretization error is not a separate tax on top of the closure problem — written correctly, it is a term inside the unresolved stress.
• The exact target exists. Feed it back and a coarse turbulence simulation tracks the filtered DNS to machine precision, in 1D and 3D.
• It is non-symmetric and non-local. The symmetric, local closures we have been building for sixty years cannot reach it even in principle — at best they reach the 2.86% rung of the ladder.

One more confession is still owed. This whole story treats space as guilty and time as innocent: the equations above are discrete in $x$ and blissfully continuous in $t$. But nobody integrates ODEs exactly either — and the time step, it turns out, is also a filter. That, however, is a story for another paper.

All numbers in this post are from the paper: the Burgers ensemble (13 500-cell DNS, 1000 initial conditions, central scheme) and the 3D decaying-turbulence experiment (500³ DNS onto a 100³ grid, Gaussian filter of width Δ = 2h). Everything is reproducible with ExactClosure.jl.