The mesh that makes its own mass

The transport theorem holds for a region moving at any velocity, not only for one made of fluid. This collection has an essay about the identity itself, and about how three boxes drawn in one flow give three different answers that all reconcile.

That freedom is not decorative. It is what makes a moving-mesh calculation possible at all — a wing that flaps, a valve that closes, a piston, a free surface, a rotor passing a stator — and every one of those is computed by writing the conservation laws on regions that move at a velocity somebody chose.

This essay is about the price of discretising it.

The identity, and the geometric statement inside it #

For a region $V_w(t)$ whose boundary moves at $\mathbf{w}$,

$$\frac{d}{dt}\int_{V_w} f\,dV = \int_{V_w} \frac{\partial f}{\partial t}\,dV + \oint f\,(\mathbf{w}\cdot\mathbf{n})\,dS,$$

which is the theorem with the material velocity replaced by whatever the region is doing.

Three independent computations here: the left side by differencing the integral over the moving cell, the first term on the right by quadrature of the partial derivative over the same cell, and the flux by five-point Gauss along each edge. Two of them add to the third to seven figures, with an unsteady scalar and a mesh velocity that has nothing to do with it.

Now put $f = 1$. Every term about the fluid disappears and what is left is a statement about nothing but geometry:

$$\frac{dV}{dt} = \oint (\mathbf{w}\cdot\mathbf{n})\,dS.$$

A region’s volume changes at the rate its boundary sweeps volume. It is obvious, it is exact, and it is where the trouble is: a scheme has to satisfy it in its own arithmetic, not merely in the limit of a small step.

The three forms of the theorem that matter are the same identity with three integrands, and this collection has already used all of them. Put $f = \rho$ and it is continuity — mass has nowhere to go. Put $f = \rho\mathbf{u}$ and it is the momentum theorem, which is how every force on this site is computed. Put $f = \rho(e + u^2/2)$ and it is the energy equation. Three founding accountings, one identity, and one geometric statement underneath all of them.

What happens if it does not #

Consider a fluid at rest with a uniform density, and a mesh that moves.

Nothing is flowing. The scalar is one everywhere and should stay one for ever, at any step size, under any mesh motion whatever, because nothing is happening.

A finite-volume scheme does not know that. It holds a mass in each cell, updates it by the fluxes through the cell’s moving faces, and divides by the new volume to recover the scalar. If the fluxes it computes do not add to the volume change it computes, then the mass it ends with does not match the volume it ends with, and the scalar moves.

The rule the identity actually asks for is the swept volume: between one step and the next each edge sweeps a quadrilateral, and the signed area of that quadrilateral is exactly the integral of $\mathbf{w}\cdot\mathbf{n}$ over the step. Nothing is approximated. The four of them add to the change in the cell’s area to the last bit of double precision, because it is a shoelace identity rather than a numerical result.

The rule anybody would write down instead is the face’s mean velocity over the step, dotted into its mid-time normal. It is second-order accurate in the step, it needs no geometry beyond what the scheme already has, and it is entirely reasonable.

On a rigid translation the second rule is exact. That is not luck: a mesh that only slides gives every face the same velocity, and the outward normals of a closed polygon sum to zero, so the fluxes cancel identically. Which is why nobody notices. Almost every first test of a moving-mesh code is a translation.

On a mesh that turns it leaves a quarter of the volume change unaccounted at a step of 0.4.

The reason is worth having in words rather than as a number. The swept quadrilateral is bounded by two versions of the edge — where it was and where it went — and by the two paths its endpoints travelled. The midpoint rule replaces that quadrilateral by a parallelogram: one edge, one mean velocity, multiplied. When the edge does not change length or direction the two are the same figure. When it rotates and stretches they differ at second order in the step, and the difference is exactly the amount of volume the scheme fails to account for.

So the error is a geometric one and has nothing to do with the fluid. It is the same for water, for air, for a scalar and for a full Navier-Stokes solution, and it is present at exactly the same size whether the flow is fast or the domain is empty.

The measurement #

The swept rule leaves the scalar at one to within a part in 10¹⁵ at every step of the two time units. The midpoint rule moves it by 2.2 parts in ten thousand — on a mesh motion that begins and ends in exactly the same place, in a fluid that is not moving.

Two parts in ten thousand sounds small. It is not, in context: it is a change of a quantity that is not changing, and it is larger than most of the physical effects a simulation of that kind is run to resolve. Worse, the excursion is smooth. It rises, peaks, and comes back, which is what a physical transient looks like, and there is nothing in the output that says it came from the grid.

The mesh motion is periodic and the cell’s area returns to its starting value, which is what makes the test unambiguous. A scheme that has changed a uniform scalar after one period has not made a small error in a changing quantity; it has changed a quantity that did not change.

There is a detail in the measurement that is easy to skip and is the reason the number can be trusted. The two rules are given the same vertex displacements: the mesh velocity is the mesh’s own motion divided by the step, not a differenced approximation to a continuous velocity field. So the comparison is between two ways of turning a known geometry into fluxes, with no differencing error anywhere in either, and the whole of the difference is the accounting.

That matters because an earlier version of the midpoint rule here used a finite difference to get the face velocity, and it left a residue of 6·10⁻¹⁰ even on a rigid translation — where the correct answer is exactly zero. Six parts in ten billion is small enough to be invisible and large enough to obscure the very contrast being drawn, and it came from the instrument rather than from the rule.

An error, not an identity #

The midpoint rule’s free-stream error falls as the square of the step — measured exponent 1.977 over four halvings — which is 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.

The swept rule has nothing to refine. It sits at round-off at every step size tried, which is the same signature this collection found in the area a two-dimensional incompressible flow must not change: an integrator that conserves by construction returns the invariant at step sizes far too coarse to be accurate, and one that conserves by accuracy returns it as a convergent error.

The parallel is exact and worth stating, because the two problems look unrelated. There the conserved quantity was the area of a patch of dye and the scheme was an integrator; here it is the volume of a cell and the scheme is a face rule. In both, the failure is invisible in the picture, the error obeys the scheme’s own order, and the only instrument that separates the two schemes is a quantity carried alongside deliberately.

Where else the same accounting appears #

Free-stream preservation is stated for moving meshes because that is where it was first noticed, but the requirement is more general and it is worth naming the other places it shows up.

A sliding interface between a rotor and a stator has faces on one side that do not line up with faces on the other. Whatever rule reconciles them has to conserve volume across the interface exactly, or the two sides disagree about how much fluid passed and the disagreement appears as a source.

A cut cell, where a body’s surface passes through a Cartesian grid, has a volume and a set of face areas computed by separate pieces of geometry. If they do not satisfy the divergence theorem discretely — the face area vectors of a closed cell summing to zero — then a uniform flow through the cell is not uniform when it leaves.

A free surface whose height is advanced by one equation and whose cell volumes are recomputed by another has exactly this problem, and it is the reason a volume-of-fluid method computes the flux volumes geometrically rather than from a velocity times an area.

In each case the pattern is the same: two pieces of a code compute a geometric quantity by different routes, both routes are accurate, and the answer depends on their agreeing exactly rather than accurately.

Why this is not an argument for the swept rule everywhere #

The swept rule is exact for the geometry and costs a little more arithmetic. That is not the interesting part, and treating this as a recommendation between two formulas would miss what the measurement is about.

What it is about is the class of statement being made. “The scheme is second-order accurate” is a statement about a limit. “A fluid at rest stays at rest” is a statement about an identity, and the two are different in kind: the first can be satisfied by any consistent discretisation and the second cannot. A scheme’s conservation properties are decided by what it satisfies exactly, at the step size it is actually being run at, and refinement studies do not reach that question because every consistent scheme passes them.

There is a second reason the distinction is not a formality. The free-stream error scales with the mesh velocity and not with the flow, so it does not shrink when the interesting physics is slow. A rotor spinning fast in a nearly quiescent fluid is the worst case, and it is a common one.

The freedom the identity leaves #

The other half of the theorem is the reason all of this is allowed. Because the identity holds for any $\mathbf{w}$, the physics cannot depend on which one is chosen — and under the swept rule three completely different mesh motions give the same answer to four parts in 10¹⁵.

That is a freedom of exactly the kind an exact constraint always leaves: the mesh velocity is a whole vector field that the answer must not depend on, and a scheme’s correctness is partly the statement that it does not. A scheme failing the geometric conservation law has made the answer depend on the mesh velocity, quietly, and the amount it depends on it is the number measured above.

It is also why the failure has the reputation it has. A code that violates the law will give different answers on two grids that describe the same motion differently, and the difference will look like grid dependence — which is a familiar thing, is usually blamed on resolution, and is usually chased by refining. Refining does make it go away, at second order, which confirms the diagnosis and hides the cause.

What a reader of a moving-mesh result should ask #

The measurement suggests three questions, and none of them is about the physics being modelled.

Was the case with nothing in it run? A uniform field, a moving mesh and no flow. It costs one run and it is the only test that isolates this.

How fast is the mesh compared with the flow? The error scales with the mesh velocity, so the worst case is a fast-moving boundary in a slow fluid — a valve seating, a rotor windmilling, a body being started. A code validated on a fast flow past a slowly deforming boundary has not been tested where it is weakest.

And does the reported transient have a timescale the mesh motion also has? A free-stream error tracks the mesh’s own motion, so it appears at the frequency the boundary is moving at. In a periodically flapping or rotating case that is exactly the frequency the result is being measured at, which is the least convenient possible coincidence.

What to carry away from it #

The general form is short and is not about meshes.

An accounting identity that is exact in the continuum has to be exact in the discrete accounting as well, or it is not an accounting at all. The mass in a cell and the volume of a cell are computed by two different pieces of a scheme, and if those two pieces disagree about geometry the disagreement appears as physics. The same statement covers a conservative scheme’s flux at a shared face, a moving-boundary condition, and the geometry of a cut cell.

The test for it is to run the case where nothing should happen. That is the whole of free-stream preservation: put in a uniform field, do not flow it, and see whether it comes out uniform.

And the failure is not detected by accuracy. It is second order, it converges, and every convergence study a careful person would run reports success. What finds it is a quantity that is known exactly in advance, which is the same discipline every figure in this collection is drawn under and the reason a smooth picture of a flow proves nothing on its own.

The identity is older than the difficulty #

One more thing is worth saying, because it changes how the requirement should be read.

The geometric conservation law is not a numerical trick invented for finite volumes. It is the $f = 1$ case of a theorem that predates computing entirely, and every continuous derivation on this site has been using it silently. When a control volume is drawn round a body and the momentum theorem applied, the volume is fixed and the term vanishes; when the volume is material, the term is the divergence of the velocity and the incompressibility condition kills it. The only situation in which it has to be carried explicitly is the one where the region is moving at a velocity nobody chose for physical reasons — which is exactly the moving-mesh case, and is why the difficulty arrived with the method.

That also settles what kind of requirement it is. It is not an extra condition imposed on a scheme from outside; it is the scheme being asked to satisfy the identity it is a discretisation of. A method that fails it is not solving a slightly different problem accurately. It is solving a problem in which the geometry does not close, and the flow it reports has a source term in it whose size is a property of the grid’s motion.

The distinction shows up in what fixes it. Better accuracy does not, because the failure is not an accuracy failure; a smaller step reduces it because the source is second order in the step, which is the same thing every truncation error does and tells nobody anything. What fixes it is computing one geometric quantity once and using it in both places, which costs nothing and is not a numerical improvement at all.

A note on what “conservative” means here #

The word conservative is used loosely enough in this subject to be worth pinning down, since this page turns on a distinction it hides.

A scheme is usually called conservative when it is written in flux form: each cell’s content changes by the fluxes through its faces, and a flux leaving one cell is exactly the flux entering its neighbour. That property guarantees that the total of a conserved quantity over the whole domain changes only by what crosses the outer boundary, whatever the scheme does inside — which is a strong and useful guarantee, and it is what makes shock capturing work.

It does not guarantee anything about a moving mesh. Flux form says the fluxes agree between neighbours; it says nothing about whether the fluxes agree with the volumes. A scheme can be perfectly conservative in the flux sense and still create density out of grid motion, because the mass it moves is right and the volume it divides by is not.

So there are two separate properties, and both are needed. Flux consistency makes the total right; geometric consistency makes the density right. The first is checked by summing over the domain and the second by running the case where nothing happens, and a scheme can have either without the other.

That is worth carrying because the two are often conflated in a code’s documentation, and a claim of conservation without a free-stream test is a claim about only one of them.

What is not claimed #

One cell is not a scheme. Everything here is measured on a single quadrilateral with prescribed corner motion. A real code has cell-to-cell fluxes, a time integrator, and a geometry that is assembled rather than given, and the geometric conservation law in that setting is a condition on how the swept volumes are computed consistently between neighbours as well as within a cell.

Second order is not the only order the failure can have. The 1.977 measured here belongs to this face rule and this motion; a different rule fails at a different order, and some fail at first, which is far worse and far more obvious.

The identity does not make a scheme conservative. Satisfying the geometric conservation law is necessary and not sufficient: a scheme can preserve a uniform field exactly and still lose mass through an inconsistent flux at a shared face. The test used here checks one thing and is silent about the rest.

And nothing here says the midpoint rule is wrong to use. It is second order, it converges, and at a fine enough step its free-stream error is below every other error in the calculation. What the measurement establishes is where that step is, and that it has to be established rather than assumed.