A parameter that is also a price

The essay that picked the diffusion on purpose ends with a scheme and a warning. The scheme adds ε·ξ·Pe of artificial diffusion, with ξ = coth(Pe) − 1/Pe, and on the boundary-layer problem it was derived from the computed solution is the exact one at every grid point, to 2.4·10⁻¹⁷. The warning is that the exactness belongs to the problem rather than to the scheme: put a smooth source on the same operator, at the same ε and on the same grid, and the tuned scheme’s error is 6.69·10⁻² against 1.60·10⁻³ for adding nothing at all.

Forty-two times worse is a fact about one situation. It was read off one diffusion, on one manufactured problem, and offered as a warning rather than as a law — which leaves the question a code actually has open. A code deciding whether to stabilise an element does not want to be told that the tuning can be harmful. It wants to know how harmful, and as a function of what.

There is an exact answer and it needs no new measurement, because the quantity that answers it is the parameter that was already computed. The tuned scheme adds ξ·h/2 of diffusion. Upwinding adds h/2, which is the same quantity with ξ set to one. On a solution smooth enough that neither addition repairs anything, the error each scheme commits is dominated by what it added and is linear in it — so the two errors should stand in the ratio ξ to 1, and the number that buys the exactness on one problem is the invoice on the other.

The two schemes differ by one factor and their errors differ by the same one #

The measurement is the ratio of the tuned scheme’s worst nodal error to upwinding’s, on the manufactured smooth solution, at a fixed grid of 127 points and five diffusions. Against it stands ξ evaluated at the cell Péclet number each of those diffusions produces. Nothing is fitted; the right-hand column is a closed-form function of Pe and the left-hand one is two solves.

The two columns agree to within 0.71% at every one of the five and within 0.3% at four of them — −0.71%, +0.13%, +0.27%, +0.26% and +0.18%, in that order. The deviation has no fixed sign, which is what rules out the obvious way of being wrong about this: a systematic offset would be a missing constant, and there is none to find. The ratio is ξ, and the error a stabilisation commits on a problem with nothing to stabilise is the stabilisation parameter times the error of the crudest scheme in the field.

That is worth separating from the two statements it sits between, because both were already known and neither is this one. That upwinding is inaccurate on a smooth problem is ordinary. That the tuned scheme is inaccurate there was the finding the essay that derived it closed on. What is new is that the second is the first multiplied by a number the code has already computed for a different purpose — so the damage does not have to be measured, estimated or bounded. It is read off the same formula that produced the parameter, at the same argument.

One Péclet number, three grids, and the agreement sharpens #

A ratio that holds along one row of a table is a coincidence until it is checked in the direction the table does not run. ξ is a function of Pe alone, so the claim carries a prediction that has nothing to do with the grid: two different meshes at two different diffusions that produce the same Péclet number must give the same error ratio.

Pe = 0.7813 arises at ε = 0.02 on 31 points, at ε = 0.01 on 63, and at ε = 0.005 on 127. The ratios there are 0.2548, 0.2523 and 0.2511, against a single ξ of 0.2504. They agree with it to 1.77%, 0.76% and 0.27% — the same number three times, from three problems that share no mesh and no diffusion.

The agreement is not merely holding. It is sharpening, monotonically, and the reason is the term that has been neglected. Central differencing’s own truncation error does not vanish; it is simply small beside what the two schemes added, and its share of the total falls as the grid is refined. At Pe = 3.9063 the same check reads 0.7487, 0.7461 and 0.7452 on 63, 127 and 255 points against ξ = 0.7448 — 0.53%, 0.18% and 0.06%. The identity is exact in the limit and the approach to it is visible at three grids, which is the strongest form the claim can take on a page.

That also settles what the remaining discrepancy is made of. It is not a defect in the identity and not a rounding: it is the second-order error the unstabilised scheme would have committed anyway, riding along underneath a first-order one and diluting the ratio by a fraction of a per cent. Where the two are closest in size the deviation is largest. Across twenty-four cases in the convection-dominated range — four grids and six diffusions, every one of them a case where the tuned error exceeds central differencing’s tenfold — the largest departure from ξ is 2.02%, at 31 points and ε = 0.01, which is the coarsest grid at the middle of the range.

The curve is a price list, and the unit is one element #

Once the ratio is the parameter, the curve of the parameter stops being a description of a limit and becomes a table of charges. It was drawn in the essay that derived ξ to make a point about the two schemes it contains at its ends. Read the other way round — as the fraction of upwinding’s damage that a tuned element commits on a problem it was not derived for — the same curve answers a question no scheme comparison has asked.

The reading is direct and needs no interpolation. At a cell Péclet number of 0.1 the tuning costs 3.3% of upwinding’s error on a problem with no layer in it. At Pe = 0.3 it costs 9.9%. At Pe = 1 — the threshold where central differencing stops being an M-matrix and starts leaving the interval the equation guarantees — it costs 31.3%. At Pe = 3 it costs 67.2%. At Pe = 10 it costs 0.900000, which is 1 − 1/Pe to six places because coth has already reached one at that argument; at Pe = 30 it costs 96.7%; and by Pe = 100 it costs 99.0% and the scheme is upwinding with a more expensive derivation behind it.

Two things follow from having the price as a curve rather than as an anecdote, and both are about where the money goes. The first is that the charge is negligible exactly where the parameter is negligible: the small-argument behaviour is ξ ≈ Pe/3, so at Pe = 0.01 the price is a third of a per cent, and a code stabilising a well-resolved region is paying almost nothing for the insurance. The second is that the charge rises fastest through the decade around Pe = 1 — from 3.3% to 67.2% as Pe goes from 0.1 to 3 — which is precisely the range in which a code cannot tell by inspection whether the element needs stabilising. The parameter is cheap where the decision is easy and expensive where it is hard.

This is the same shape as a knob whose setting is itself the accuracy, where the tolerance handed to a compression turns up unchanged as the residual of the solve that uses it. There the number a caller types is the error it will get. Here the number a formula computes is the fraction of the worst available error that will be charged, and in both cases the useful consequence is that no separate error estimate has to be run: the setting and its consequence are one quantity read twice.

The price has a ceiling, and the ceiling is upwinding’s own error #

ξ runs from zero to one and never reaches one. That is an elementary property of coth(Pe) − 1/Pe and it is asserted in the field’s own machinery at Pe = 10, 30, 100 and 300. Multiplied through the identity above it becomes a statement nothing in this field has made: the harm the tuned scheme can do on a problem it was not derived for is bounded, and the bound is the harm upwinding does on the same problem.

Swept over four grids and twelve diffusions — 48 cases, with cell Péclet numbers from 0.0098 to 312.5 — the largest tuned-to-upwind ratio measured is 0.998415, at fifteen points and ε = 10⁻⁴. The ratio approaches one and stays under it, at every case tried.

The bound is worth something only if the quantity it bounds is itself finite, and it is, for a reason that belongs to the other scheme. Upwinding’s added diffusion is h/2, which does not depend on ε at all, so on a fixed grid upwinding’s error on the smooth problem stops moving once the flow dominates. At 31 points it reads 8.605·10⁻², 9.280·10⁻², 9.500·10⁻², 9.611·10⁻², 9.704·10⁻², 9.730·10⁻² and 9.739·10⁻² as ε falls through 0.05 to 2·10⁻⁴ — a factor of 250 in the diffusion and 13% in the error, converging on a number fixed by the mesh.

So the tuned scheme’s damage on the wrong problem saturates too, and against the same ceiling. Measured as a multiple of what adding nothing would have cost, at 31 points, it runs 6.0, 15.2, 27.7, 41.9, 53.3, 57.1, 59.0 and 60.2 across those same eight diffusions, while upwinding’s own multiple runs 58.5, 59.7, 60.0, 60.2, 60.5, 60.6, 60.6 and 60.6. The tuned column is climbing towards the upwind column and flattening against it. It is not growing without bound, and a code cannot be hurt by this tuning by more than it would have been hurt by the crudest scheme available.

The ceiling is a consolation with a sharp edge on it, and the edge is worth stating in the same breath. Bounded by upwinding is not the same as small. Sixty times the error of the scheme that adds nothing is what the bound permits on this grid, and the whole of the argument for stabilising was that upwinding was too blunt to use. A method whose worst case is the method it was built to replace has a guarantee, and the guarantee is not the one it was sold on.

Below one diffusion the price is negative #

The other end of the curve does something the warning this essay started from does not allow, and the carrier generator found it by refusing to draw the comparison at all.

Every one of the five carrier figures asserts, before it draws anything, that the tuned scheme’s error exceeds ten times central differencing’s at each grid size — the claim the picture is a picture of. At ε = 0.05 that assertion fails, reporting 0.00875 against 0.00147. At ε = 0.1 it fails again, at 0.00345 against 0.00132. And at ε = 0.2 it fails with the two quantities in the opposite order: 9.13·10⁻⁴ against 0.0011, the tuned scheme ahead of adding nothing. At ε = 0.4 the same comparison reads 5.26·10⁻⁴ against 9.17·10⁻⁴.

So the sentence the tuning harms every problem it was not derived for is false, and it is false at the diffusive end rather than in some pathological corner. Where the diffusion is large enough, ξ·h/2 is a small enough addition that the error it introduces is smaller than the second-order truncation error it happens to damp, and the stabilised answer is the better one. The refusals are the measurement rather than an obstacle to it: three diffusions at which no drawing is produced, each because the inequality the picture exists to show has stopped being true.

Bisecting the comparison gives the crossover to six figures. At 31 points the tuned scheme and central differencing commit the same error at ε = 0.183128; at 63 points at 0.183063; at 127 points at 0.183046. The figure’s own ten-times threshold, which is what actually decides whether a drawing is allowed, sits at ε = 0.030623, 0.030948 and 0.031056 on the same three grids.

The crossover is fixed in ε and not in the Péclet number #

Those three crossover diffusions are the same number to four significant figures, and the cell Péclet numbers at which they occur are not: 0.08532, 0.04268 and 0.02134, halving with each refinement exactly as h does. A fifteen-point grid gives 0.182278 and Pe = 0.17144, completing the pattern in the other direction.

That is the opposite of what the identity above would suggest on a careless reading. If the tuned scheme’s error is ξ times upwinding’s, and ξ is a function of Pe, then everything about this comparison ought to be a function of Pe — and the sign of the comparison against central differencing is not. It is decided by ε alone, and the mesh has nothing to do with it.

The reconciliation is that the identity is a ratio against upwinding and the crossover is a ratio against central differencing, and those two denominators behave differently under refinement. At small Péclet number ξ ≈ Pe/3 = h/6ε, so the tuned scheme’s added diffusion is ξ·h/2 = h²/12ε — second order in h, with ε in the denominator. Central differencing’s own error on a smooth solution is second order in h too. Two second-order quantities have a ratio that does not depend on h, and what is left in it is the ε.

That prediction is checkable and it holds. Holding ε fixed and refining, the tuned-to-central ratio converges rather than growing: at ε = 0.2 it reads 0.830, 0.829, 0.830, 0.830 and 0.830 at 31 through 511 points; at ε = 0.1, 2.610 then 2.615 four times over; at ε = 0.02, 15.213, 15.647, 15.760, 15.791 and 15.798; at ε = 0.005, 41.919, 55.408, 61.395, 63.226 and 63.703.

The ratio grows with refinement and it grows to a limit, and the limit is a function of ε. That is a correction to the natural reading of the original measurement, which reported 41.9, 55.4 and 61.4 at three grids and attributed the growth to the tuned scheme being first order where central differencing is second. It is first order at those grids, because they have not reached the small-Péclet regime; two more refinements take it to second order and the growth stops. The scheme does not become arbitrarily bad relative to doing nothing. It becomes a fixed multiple worse, and the multiple is set by the diffusion in the problem rather than by anything the solver chooses.

What a code can do with a price list #

The measurements above turn a warning into an arithmetic, and the arithmetic is cheap enough to run per element. Case by case:

The parameter is already computed, so the estimate is free. Every practical streamline-diffusion code evaluates ξ, or a cheaper approximation to it, from a local Péclet number. Multiplying it by an estimate of what upwinding would cost gives what the stabilisation costs where it is not needed, and both factors are in hand before any solve is done.

The comparison worth printing is the one against upwinding, not against nothing. Running the same problem unstabilised is the honest check and it costs a second solve. The identity says what that second solve would have found, to within a per cent, at no cost — and it says it in the ratio a reader can act on, because a scheme that is 25% of upwinding is one a reader has a reference point for and a scheme that is fifteen times central differencing is not.

The number is per element and it does not average. Cell Péclet numbers in a real problem span decades across a mesh, and ξ spans nearly its whole range across those decades, so a single global figure for what the stabilisation cost is the wrong statistic: a quantity varying by two orders of magnitude across the population it is summarised over has a mean that describes nothing in it. What must not stand in for it is the residual, which measures how well the tuned system was solved rather than whether it was the right system — a small residual is not a small error, and on the manufactured problem at ε = 0.005 the tuned scheme solves its own system to a relative residual of 1.6·10⁻¹⁵ while its answer is 42 times further from the truth than doing nothing.

And the price is charged whether or not the goods are delivered. ξ looks at the cell Péclet number and nothing else. It does not look at the right-hand side, and the right-hand side is what decides whether there is a layer to stabilise, which is the whole content of the reading that names which problem a scheme answers. An element with a large Péclet number and a smooth solution running through it is charged 90% of upwinding’s error for insurance against an oscillation that was never going to happen.

What the price does not cover #

Three limits on the arithmetic above, each of which is a measurement made elsewhere rather than a caveat.

It is a price in one dimension and along the flow. The whole identity is stated on a problem whose velocity is constant and aligned with the grid. Five degrees off the grid line the tuned scheme’s relative error goes from 1.2·10⁻¹⁴ to 6.9, and by twenty degrees it is worse than the upwinding it was built to improve on — which means the ceiling established here is a ceiling for the aligned case and is not one in general.

It prices the wrong-problem error and not the crosswind one. Streamline diffusion adds a rank-one tensor that annihilates every direction across the flow, which is the design and also the defect: the price list says what the added diffusion costs where it acts, and says nothing about the interval excursions in the direction it does not act at all.

And it is a price for a parameter that was optimal somewhere else. The whole of ξ’s authority comes from a one-dimensional constant-coefficient problem with no source term, which is the same inheritance a regularisation parameter chosen on a smaller problem carries and the same one four knobs on one problem measured across four unrelated methods. The difference here is that the inheritance has a receipt. ξ is not a parameter that might be wrong by an unknown amount on a problem it was not derived for; it is a parameter that is wrong by exactly itself, times upwinding’s error, and that is the rarest thing a transferred parameter can have.

The last observation is about which competence is being bought. The tuned scheme knows the answer is a geometric sequence because somebody put that knowledge into the formula, which is the shape a smoother built for one direction has in the multigrid field — 0.037 on the problem it was aimed at and 0.967 on the same problem turned through a right angle. The difference between the two cases is the one this page is about. That smoother’s failure is unpriced: 0.967 is a number somebody had to measure. Here the failure has a closed form, it is the same closed form as the success, and reading it costs a call to coth.