Your time integrator
is a filter

An interactive companion to the preprint:

Part one and part two carried a quiet assumption. The grid was coarse and the scheme was imperfect, but time was continuous: the exact equation for the cell average began with $\partial_t \bar{u}$, a real time derivative. No code solves that equation. A code takes steps.

The usual view is that the time step adds a second error on top of the first, to be kept small enough not to matter. This post takes the other view, the one that has worked twice already. The time step is one more box.

A brick in space-time

Part one integrated a derivative over an interval and was left with its two end points. Part two integrated over a cube and was left with its six faces. Now draw a rectangle in the $(x, t)$ plane, one grid cell wide and one time step $\Delta t$ tall. Integrate the conservation law $\partial_t u + \partial_x F(u) = 0$ over it, and what remains sits on its four sides:

The top and the bottom of the brick carry $\bar{u}$, the cell average of part one, at the end and at the start of the step. The left and right sides carry $\langle F(u) \rangle$, the flux at a cell face, averaged over the step from $t$ to $t + \Delta t$. As before, $\partial_x^h$ is the difference between the two faces, divided by $h$. Nothing is approximated. What the cell gains between bottom and top is what crossed its two sides in between.

Here is a flow, drawn in space-time, and a brick. Make it as tall as you like.

The first two numbers are computed separately, one from cell averages and the other by adding up fluxes, and they agree for every brick: any height, any place, straight across a shock.

Now read the equation again. Its first term is the forward-Euler difference. So forward Euler is not a first-order approximation of the conservation law. It is an exact statement of it, for a time step of any size, provided the flux is the one averaged over the step. No $\mathcal{O}(\Delta t)$ is hiding in it.

The first-order error enters with the third number. At time $t$, a code cannot know what will cross the sides during the step. It knows the flux now, so it holds that value for the whole step. This is the left-rectangle rule for an integral in time, and the third number is what it gives. Make the brick short and the third number joins the other two. Make it tall and it goes its own way.

One more bracket

Part one replaced the flux $F(u)$ by a numerical flux $F_h(\bar{u})$ and called what was thrown away the residual flux $\tau$. Do the same here. The exact flux is now $\langle F(u) \rangle$, and the numerical flux is evaluated at the start of the step:

The first three brackets are those of part one, taken at the start of the step. The fourth is everything the time step adds.

Look at what is missing. The field itself is never averaged in time: $\bar{u}$ at the top and at the bottom of the brick are snapshots. The time integrator filters the flux, not the solution. So there is no new unresolved stress to model along the time axis, and no new commutator. There is one quadrature error, of a flux you already had.

Exact at any time step

This is the twin simulation of part one, with a clock. The DNS runs with its own small time step. The coarse run takes one step for every 100 of them, and it is fed whichever brackets you leave switched on.

With all four brackets on, the dots land on the bars at every coarse step, to round-off. Three things to try:

• Push the time step all the way up. The CFL number now starts above two, far past any stability limit. The dots still land on the bars. Nothing in the identity knows about a stability limit: such a limit belongs to a scheme that only knows the present, and $\langle F(u) \rangle$ has seen the whole step.
• Switch off your time step and leave the other three on. This is the best that part one can offer, a closure that is perfect in space. At the smallest step it is still exact, since the coarse run and the DNS then share a clock and the fourth bracket is zero. Raise the step one notch at a time. The perfect spatial closure drifts by a tenth, then by a quarter, then by more than half, and at the two largest steps the run blows up.
• Switch to upwind, turn all four brackets off, and raise the time step again. This time the error goes down with every notch but the last. The next section says why.

The new bracket is an old acquaintance

Expand the average over the step to first order in $\Delta t$, and use the conservation law to trade the time derivative for a space derivative:

(The last step is for the convective part of the Burgers flux.) This is a diffusive flux, with a viscosity $\tfrac{1}{2} u^2 \Delta t$ that cannot be negative. To leading order, forward Euler is short of this much diffusion, which is why it is unstable with a central scheme.

The term has a name. It is the diffusion of the Lax-Wendroff scheme, from 1960, and it is what makes that scheme stable. Here nobody added it. It is the leading part of a residual that was defined before any expansion was made, and it comes with its coefficient, one half.

Written with the CFL number, $|u| \Delta t / h$, the viscosity is

and $\tfrac{1}{2} |u| h$ is the numerical viscosity of a first-order upwind scheme. So the time step hands back a fraction of an upwind scheme’s diffusion, and that fraction is the CFL number. This is the third experiment above. An upwind scheme dissipates too much, forward Euler too little, and a larger time step lets the second error eat more of the first, until the CFL number passes one and there is nothing left to eat.

The new bracket grows with the time step. The other three do not care about it. In the paper’s experiment with a central scheme, at a CFL number of 0.89, it is the largest of the four, larger than the other three together. A closure that only knows about space never sees it.

A closure with a clock

The twin needs a DNS running beside it. A model has no DNS, but it does have $\Delta t$. So take the simplest closure there is, the Smagorinsky term fitted to the three spatial brackets as in part one, and put the expansion next to it:

Below, no DNS data enters the coarse runs. One uses the first term only. The other has the clock term as well.

Slide the time step down and the two closures agree: at small steps the clock term does nothing, and the error that is left belongs to Smagorinsky. Slide it up and they part. The space-only closure cannot supply the diffusion that the time integrator asks for, and its run fills with wiggles. The other one does not mind the step. Then drag $c$ at the largest step, and look for the value with the smallest error. In these runs it sits close to one half.

What this asks of a closure

A closure is tied to its time step. Fit or train a model at one $\Delta t$ and run it at another, and part of what it learned was the time integrator. For closures trained on data, this is one concrete way in which the training setup and the simulation can disagree.

The target exists at any step. The four brackets give a residual flux that reproduces the filtered DNS exactly at coarse time steps. A model can be trained against it directly.

Diffusion is not all of it. In the paper, the time bracket is most of the residual at large steps but accounts for only a small part of the energy the residual removes. Much of it does not dissipate at all, and no eddy viscosity can represent that part.

This is forward Euler. A higher-order integrator is a better quadrature rule, and its bracket shrinks accordingly. Multi-stage and implicit integrators, and the Navier-Stokes equations in three dimensions, have not been worked out yet.

Three boxes

An interval gave a residual flux that knows your scheme. A cube gave a stress with nine components and a long reach. A brick in space-time gave one more bracket, and it turned out to be Lax and Wendroff’s. All three came from the same move: integrate a derivative over a box, and what remains lives on its boundary. In each case the discretization was not an error added afterwards. It was the box.