The area that must not move

This collection’s standing claim is that a wrong flow field is beautiful: streamlines are smooth whatever nonsense produced them, and neither a reader nor an author can tell a solution of the equations from a plausible function that satisfies none of them. Every generator here therefore computes an invariant and refuses to draw if it is not satisfied.

This essay is that claim one level down, about the integration rather than the field, and about the one quantity that can settle it.

The law #

Advection by an incompressible flow is exactly area-preserving. Let $\mathbf{X}(\mathbf{x}, t)$ be where the particle starting at $\mathbf{x}$ has got to, and $J$ the determinant of $\partial\mathbf{X}/\partial\mathbf{x}$. Then

$$\frac{dJ}{dt} = J\,(\nabla\cdot\mathbf{u}),$$

and $\nabla\cdot\mathbf{u}$ is exactly zero. So $J = 1$ for all time, in every incompressible flow, however violently it is stirred. A patch of dye has exactly the area it started with — not approximately, not on average, exactly — which is mass having nowhere to go written for a material region.

The flow used here is the time-periodic double gyre, written as the curl of a stream function stated in closed form, so its divergence is zero analytically and at round-off numerically. That matters: a field with a small divergence would change area at a small rate that could be mistaken for the integrator’s error.

And here is the difficulty. The patch is stretched into a filament and folded through itself, and its area does not change at any point of that — which is not visible in the picture and cannot be made visible. A scheme losing eight per cent an orbit would produce a picture indistinguishable from this one.

Why the area and not something else #

There is a reason this particular invariant is the one worth carrying, and it is worth separating from the fact that it happens to be conserved.

A two-dimensional incompressible flow is a Hamiltonian system with the stream function as its Hamiltonian: $\dot{x} = \partial\psi/\partial y$, $\dot{y} = -\partial\psi/\partial x$ is Hamilton’s equations with $x$ and $y$ as the conjugate pair. So advection in the plane is not merely like a mechanical system, it is one, and area in the $(x, y)$ plane is phase-space volume.

That brings a hundred and fifty years of machinery with it. Liouville’s theorem is the area conservation; the Poincaré recurrence theorem says a bounded steady flow brings almost every parcel arbitrarily close to where it started; and the whole apparatus of Hamiltonian chaos — islands, tori and the stochastic layers between them — applies to a stirred fluid without modification.

It also explains why the numerical question has an answer. Integrators that preserve the symplectic structure of a Hamiltonian system were developed for celestial mechanics, where the same problem arises: a planet’s orbit computed for a million years by a scheme that does not conserve phase-space volume slowly spirals, and the spiral is the arithmetic rather than the physics. The implicit midpoint rule is the simplest member of that family, and using it on a fluid is borrowing the fix.

Two schemes #

The implicit midpoint rule returns one to nine figures at every step size, including one far too coarse to be accurate. Fourth-order Runge–Kutta — the workhorse, more accurate about where the particle is by a wide margin — does not: at $h = 0.2$ its determinant is 0.999937, six parts in a hundred thousand adrift.

Neither of these is a straw man. Both are standard, both are used for exactly this job, and the one that fails is the more accurate of the two.

The reason is algebra rather than accuracy. The midpoint rule’s tangent map is the Cayley transform $(I + M/2)(I - M/2)^{-1}$ with $M$ the step times the velocity gradient. For a two-by-two traceless $M$ — and an incompressible velocity gradient is exactly traceless — both determinants come out at $1 + \det(M)/4$. They are the same number, and their ratio is one identically.

So the conservation is not a consequence of the scheme being good. It is a consequence of the scheme’s update being a particular algebraic form, and it holds at step sizes where the scheme’s trajectory is badly wrong.

Runge–Kutta’s area error falls as the fourth power of the step — its own order of accuracy showing up in a quantity it was never designed to conserve. That is the signature of an error rather than of an identity: refine, and it goes away; refine the midpoint rule and nothing happens, because there was nothing to improve.

And the pictures do not settle it #

Side by side at twelve time units, the two patches are the same patch. Every fold is in the same place, every filament has the same shape. One of them has conserved its area exactly and the other has not.

Even measuring the area does not settle it. The polygon through 720 markers on a filamented curve under-measures the area it encloses, and the shortfall grows as the folding proceeds — four parts in ten thousand by twenty time units, under the scheme that is conserving exactly. The invariant is exact and the instrument for reading it is not, which is why the determinant is carried alongside the trajectory as a variational equation rather than estimated from the markers afterwards.

What eight per cent an orbit would look like #

It is worth putting a number on how much area a scheme can lose before anybody notices, because the answer is “more than it ever will”.

The polygon measurement above drifts by four parts in ten thousand over twenty time units under a scheme that is conserving exactly, so a measured area is good to about that. A scheme losing area would have to lose more than that before the measurement could see it — and Runge–Kutta at a step of 0.2 loses six parts in a hundred thousand over two hundred time units, which is an order of magnitude below the measurement’s own noise.

So the honest position is worse than “the pictures do not settle it”. The obvious measurement does not settle it either. The only instrument that separates the two schemes is the variational equation, and carrying it is a deliberate act rather than something that falls out of the calculation.

That is the general shape of every check in this collection. An assertion has to be designed: it has to be a quantity that is exactly known, computable to much better precision than the thing being checked, and sensitive to the failure mode of interest. Area conservation qualifies on all three, which is why it is worth the extra four multiplications a step.

What is left once the positions stop meaning anything #

The gap between the two schemes’ particles grows exponentially, at the flow’s own Lyapunov rate, and reaches the width of the domain by fifty time units.

That divergence is not a bug in either scheme. It is the flow amplifying the difference between two approximations, at the rate an hour for every tenfold is about, and no amount of accuracy removes it — halving the step buys about seven more time units of agreement, and then the same thing happens.

So at long times “where is the marker” has no answer that either calculation can supply, and asking which scheme’s picture is right is asking the wrong question. The determinant still has an answer. It is one, for one of the schemes, at every horizon, and that is the only statement about the two calculations that survives the chaos.

That is the practical content of the whole essay. A chaotic flow destroys pointwise accuracy on a timescale set by its own physics; what it cannot destroy is a conservation law, because the law is about the map rather than about any particular trajectory.

How the determinant is carried #

The determinant above is not estimated from the markers; it is integrated alongside the trajectory, and the distinction is the reason the numbers are nine figures rather than four.

Alongside $\dot{\mathbf{x}} = \mathbf{u}(\mathbf{x}, t)$ runs the variational equation for the tangent map $\mathsf{J}$,

$$\dot{\mathsf{J}} = \nabla\mathbf{u}(\mathbf{x}, t)\,\mathsf{J}, \qquad \mathsf{J}(0) = \mathsf{I},$$

which says how an infinitesimal displacement is carried by the flow. Its determinant is the local area factor at the particle, exactly, with no polygon and no finite difference anywhere.

For the comparison to be fair, each scheme’s variational equation has to be integrated by that scheme — the quantity of interest is the derivative of the discrete map, not of the exact one. For the midpoint rule that is the Cayley transform above; for Runge–Kutta it is the same Butcher tableau applied to the augmented system, which is exactly the derivative of the RK4 step and not an approximation to it.

The first version of this calculation differenced the discrete step with a perturbation of $10^{-6}$ instead, and at $h = 0.01$ over twenty thousand steps the accumulated round-off of that differencing reached eight parts in ten million — larger than the scheme’s own area error, and reported as it. The determinant then got worse with refinement, which is the signature of measuring the instrument.

The identity, and what it does not say #

The conservation is about the determinant, and it is worth being clear how little that constrains.

A two-by-two matrix of determinant one still has two singular values, $\sigma$ and $1/\sigma$, and $\sigma$ can be anything. So an area-preserving map may stretch a patch by a factor of ten thousand in one direction as long as it squeezes by ten thousand in the other — which is exactly what a chaotic flow does, and is why the dye patch above becomes a filament while keeping its area.

That means the conservation law says nothing at all about mixing, stretching, folding, or how far apart two particles get. It is one number out of three that describe the local map, and it is the only one that is fixed.

Which is why the check is worth so much and costs so little. A quantity that is fixed regardless of what else the flow is doing can be monitored through violent, chaotic, badly resolved motion and still mean something — and a quantity that varied with the flow could not distinguish a wrong integration from an interesting one.

Which scheme to use, then #

The answer is not “always the symplectic one”, and it is worth being careful.

For a short integration where the trajectory matters, Runge–Kutta is better: it is more accurate about position by roughly two orders at these step sizes, and its area error over ten time units is below anything that would be noticed.

For a long integration of a passive tracer, the area-preserving scheme is better, and not because it is more accurate. It is because its errors do not accumulate in the one quantity that has a physical meaning — a scheme whose determinant drifts to 0.98 has, after long enough, produced a dye patch that is not a dye patch of anything.

For a compressible flow the question changes rather than goes away. There the determinant is not one but $\rho_0/\rho$, and a scheme that respects that is respecting mass conservation along a particle path — which is the same statement incompressible is not a property of the fluid makes about what the condition actually asserts.

And for anything whose conclusion is a statistic of the flow map — a mixing rate, a residence-time distribution, a Lyapunov exponent, a finite-time coherent structure — the conserved quantity is the whole business, because those statistics are computed from the map’s stretching and the determinant constrains it. A scheme that quietly gains area reports mixing that is not there.

The site’s own version of this problem #

Every figure on this site is a flow field drawn from a computation, and the reason each generator ends by asserting something is precisely the argument above. Four real errors were caught that way before a single essay existed — an unrotated doublet that made an aerofoil surface leak at incidence, a tangency check that tested the wrong thing, an unconverged Poisson solve, and a circulation whose sign made the wing fly downwards — and not one of them looked wrong on screen.

The generalisation is worth stating as a rule rather than as a war story. A picture of a flow is a picture of a smooth vector field, and smoothness is cheap: almost any function is smooth, and almost no function is a solution. So the information content of a flow picture is nearly zero, and everything that makes it worth drawing is in the assertions behind it.

That is the same argument as a smooth picture proves nothing makes about streamlines and pathlines being confused, and the same argument the flow with the least energy makes about a whole family of fields that look like solutions and are not. This essay adds the integration to the list: the flow can be right and the carrying of a particle through it wrong, in a way no drawing distinguishes.

Keeping the check readable, and what it gives away free #

The determinant is a difference of products of the tangent map’s entries, and those entries grow exponentially, so past a few tens of time units the arithmetic is subtracting numbers of order $10^{20}$ to obtain one. The invariant stays exact and the way of reading it does not — which is the instrument failing rather than the law, and it has a standard repair.

Factor the tangent map at every step into an orthogonal part and an upper triangular one with positive diagonal. The orthogonal part carries the orientation and has determinant one by construction; the triangular part carries all the stretching in its two diagonal entries. Reset the map to the orthogonal factor, and accumulate the logarithms of those two entries.

The determinant is then the exponential of a running sum of small numbers instead of a difference of enormous ones, and the conservation law becomes the statement that the two accumulated logs are equal and opposite. That is checkable at any horizon whatever, because nothing in it ever gets large.

And the same factorisation hands over the quantity this essay said the determinant could not constrain. The first accumulated log divided by the elapsed time is the finite-time Lyapunov exponent — the stretching rate. One decomposition supplies both the invariant that must not move and the number that says how violently everything else is moving.

What is not claimed #

Nothing here is a statement about the flow’s own vorticity. In two dimensions vorticity is carried unchanged on a particle, so a scheme that preserves area also preserves the vorticity a parcel is holding — which is a second conserved quantity worth monitoring, and is the spin a parcel keeps rather than the area it occupies.

Two dimensions. The Cayley argument works because $M$ is two by two and the traceless condition makes both determinants equal. In three dimensions the midpoint rule is still symplectic for a Hamiltonian system and still volume-preserving for a divergence-free one, and the proof is different.

The flow is a model. The double gyre is a standard test flow for exactly this kind of question and is not a solution of the Navier–Stokes equations. Nothing here depends on it being one — the identity holds for any solenoidal field — and nothing here says how a real turbulent field behaves.

The Newton iteration is solved to round-off. The midpoint rule’s conservation is a property of the exact implicit map, and a half-solved iteration is a different map with no such property. At a step of 0.4 the iteration does not converge and the determinant comes out at 4.0; the results above are at 0.2 and below, where twelve Newton steps reach round-off.

Nothing here is about the flow field’s own accuracy. The double gyre is exact by construction, so every error discussed is the integrator’s. A computed velocity field carries its own errors, and a scheme that conserves area exactly while advecting a field that is not solenoidal conserves an area that means nothing — which is the reason every figure here is checked against a conservation law before it is drawn at all.

And the determinant is read directly here, so it becomes unreadable past about forty time units in this flow. The section above says what to do instead; none of the numbers in this essay were obtained that way.