Exact along one axis

The diffusion that makes the answer exact is the sharpest result on this site. Add ξ = coth(Pe) − 1/Pe worth of artificial diffusion to a central-difference discretisation of a convection–diffusion problem, and the computed solution equals the exact one at every node, to 2.4·10⁻¹⁷, at every Péclet number. Both of the convection field’s earlier schemes turn out to be limits of the same curve.

That essay ended by noting the exactness belongs to the problem rather than to the scheme: on a manufactured smooth solution of the identical operator the tuned scheme is forty-two times worse than adding nothing. This essay is the other boundary, and it is a narrower one.

Everything in that result was one-dimensional.

The cone the tuning is worth having in narrows as the problem gets harder, and the slider is what says so.

Across ε = 0.05, 0.02, 0.015, 0.01 and 0.005 — cell Péclet numbers of 0.63, 1.56, 2.08, 3.13 and 6.25 — the tuned scheme’s error at 0° reads 1.67·10⁻¹⁶, 3.47·10⁻¹⁷, 2.6·10⁻¹⁷, 2.39·10⁻¹⁷ and 1.88·10⁻¹⁷: exact at every one, and if anything more exact as the problem hardens.

Its error at 45° reads 0.00938, 0.037, 0.0488, 0.0661 and 0.0905 — rising by a factor of ten. Upwinding’s at 45° reads 0.0342, 0.0188, 0.00926, 0.00189 and 1.4·10⁻⁵ — falling by a factor of 2,400.

So the two curves cross, and where they cross moves. At Pe = 0.63 the tuned scheme is better off the axis by 3.6×; at Pe = 1.56 it is worse by 2.0×; at Pe = 6.25 it is worse by 6,500×. The exactness at zero degrees is a constant of the method and the penalty for leaving zero degrees is not: it grows with the very quantity the tuning exists to handle.

That is a sharper statement of this essay’s title than the one it opens with. Exact along one axis is true at every Péclet number; what changes is how much it costs to be a few degrees off, and the cost rises exactly where the scheme is supposed to be earning its keep.

The problem, and why it has an exact solution at every angle #

−ε∆u + b·∇u = 0 on the unit square, with b a unit vector at an angle θ to the grid.

Any function of the streamwise coordinate s = x·b alone solves it, because the Laplacian of such a function is its second derivative along the flow and the equation collapses to the one-dimensional one. So

$u(x, y) = (e^{s/\varepsilon} - 1) / (e^{L/\varepsilon} - 1)$, with s = x b₁ + y b₂ and L = b₁ + b₂

is an exact solution at every angle, with the boundary values read off it. The layer sits at the outflow corner and the interior is exponentially small, exactly as in one dimension.

The L in the denominator is the whole of what makes the test problem usable, and the first version of this measurement did not have it. Normalised by $e^{1/\varepsilon}$ instead — the one-dimensional constant — the boundary data at the far corner of a 45° flow is e⁴¹, every scheme is compared on numbers of order 10¹⁷, and the resulting table says nothing about schemes.

Exactness at one angle #

At θ = 0 the flow is along a grid line, the problem is one-dimensional, and the two-dimensional assembly reproduces the one-dimensional result exactly:

Worst nodal error on a 15×15 grid. That the tuned column is at the level of rounding at three Péclet numbers, from an assembly that knows nothing about one dimension, is the check that the construction is the one the previous essay measured.

Now rotate the flow by five degrees.

Twelve to fifteen orders of magnitude, for five degrees. And the fall is steeper the harder the problem is, which is the direction that matters: the sharper the layer, the more the exactness was worth and the less of it survives.

Nodal exactness is a property of the alignment, not of the scheme. In one dimension there is only one direction and the distinction cannot arise, which is why the result reads as being about ξ. The same sentence with smoother in place of scheme is a direction the smoother cannot see, and with coarsening it is coarsening in one direction only — three methods whose competence was put in by hand along an axis, and three that lose it at the same five degrees.

The cone #

That the exactness goes is not by itself a criticism — a scheme can be inexact and still be the best available. The question is where it stops being the best available.

Measured by sweeping the angle in five-degree steps and finding where another scheme takes the lead. The cone narrows as the layer sharpens, and at the sharpest layer drawn it is ten degrees wide.

The cone has a law in it #

Three numbers at five-degree resolution are three numbers. Bisecting the edge instead, on the same 15 × 15 grid:

The product settles at about 37, so the cone’s half-width is ≈ 37/Pe degrees — inversely proportional to the Péclet number over a twentyfold range, with only the mildest problem in the table off the line.

The tensor predicts exactly that. The cross term the rotation introduces is τb₁b₂ ≈ τθ for small θ, with τ = ξh/2 → h/2 at large Pe, and it competes against the physical diffusion ε = h/(2Pe). Their ratio is (h/2)θ ÷ (h/(2Pe)) = θ·Pe, with the mesh spacing cancelling — so the cone should be a function of the Péclet number alone and not of the grid. It is: at Pe = 3.125 the edge is 12.38° on a 7 × 7 grid and 12.38° on a 15 × 15 one, agreeing to four digits.

That turns the table into a rule, and the rule says something the table could not. At Pe = 12.5 the cone is 2.9° wide — narrower than the five-degree step that measured the original table, which would have reported the cone as empty. A scheme whose region of advantage shrinks as 1/Pe is a scheme whose advantage disappears exactly as the problem becomes the one it was built for, and the sweep that produced the first table is coarse enough to miss the last decade of that.

There is a repair for the cone and it is measured in the direction the diffusion does not go: adding a crosswind term buys back most of what the rotation costs, at every angle, and costs nothing where the flow is aligned. What it cannot do is widen the cone by an order — the law above is about which of three schemes wins, and a fourth scheme is a different comparison. The bound-versus-practice shape of all this — a guarantee that is exact on a constructed case and unreachable on the ones anybody has — is the bound that is never attained in the elimination field, where the constructed case is one matrix and here it is one angle.

Past it, the tuned scheme is not merely no longer best. At 45° and Pe = 3.13:

Thirty-five times upwinding’s, and twelve times the error of adding nothing at all — on a problem where at 0° the tuned scheme was exact and upwinding was wrong by 0.14.

That is a reversal, not a degradation, and the ratio between the two ends of it is 2·10¹⁷. The structural half of the same event — the sign pattern the scheme gives up off the axis — is the stencil that is not symmetric, and the algebraic half is aggregating what the matrix calls strong, where a hierarchy reading the rotated operator finds no direction in it at all.

Why the interior errors fall with the angle #

One column of the table above runs the wrong way at first reading: central differencing’s error falls from 5.2·10⁻¹ at 0° to 5.7·10⁻³ at 45°, as though rotating the flow made the problem easier.

It did not. The exact solution’s own size on the grid falls with the angle, because L grows from 1 to 1.41 and the layer moves further from every interior node: the largest exact value at an interior node is 1.9·10⁻³ at 0° and 1.5·10⁻⁴ at 45°. Relative to the solution, every scheme is worse at 45° — central differencing’s relative error is 268 at 0° and 39 at 45°, which is still an answer thirty-nine times the size of the thing it is approximating.

So the absolute errors are what the comparison between schemes rests on, and the relative ones are what says that none of the three is solving this problem at this resolution — which is a small residual is not a small error’s distinction, drawn between two normalisations of one measurement rather than between two quantities. Both are reported for that reason.

What is actually added #

The scheme adds a tensor, not a number, and the shape of it is where the rest of this essay comes from.

ε I → ε I + τ b bᵀ, with τ = ξ h / 2

In one dimension b bᵀ is the scalar 1 and the whole thing is a number. In two it is a rank-one matrix: τb₁² added to the xx diffusion, τb₂² to the yy, and τb₁b₂ to a cross term the five-point stencil does not have — so the discretisation grows a four-point diagonal stencil it did not previously need, and the matrix acquires entries at (i±1, j±1).

At θ = 0 that cross term is zero and the yy entry is zero: the added diffusion is entirely in x, which is the direction of the flow, and the scheme is the one-dimensional one applied along rows. At 45° the three coefficients are τ/2, τ/2 and τ/2, and the added diffusion is spread over five directions of which one is the flow.

The rank is the point. b bᵀ annihilates every direction perpendicular to b, so whatever τ is, the crosswind direction receives exactly nothing — which is the design, and which the next essay measures the price of.

Refining the grid fixes it, by removing the reason for the scheme #

The obvious question about any of this is whether it survives refinement. At ε = 0.02 and 45°:

Three things happen at once and only one of them is a convergence rate.

The Péclet number falls, because Pe = 1/(2ε(k+1)) and the grid is what is being refined. The exact solution’s size in the interior rises by three orders of magnitude, because the layer is now resolved and the domain is no longer almost entirely zero. And the three schemes converge on each other: at 31×31 the tuned scheme and central differencing are within 20% and upwinding is the worst of the three.

So the off-axis failure is a property of the under-resolved regime, and it disappears when the layer is resolved. That is not a reassurance. A resolved layer is the case where plain central differencing works and no artificial diffusion is wanted at all; the entire reason for a tuned scheme is the case where the layer is thinner than a cell. The failure lives exactly where the method is needed, which is the same sentence the depth phase wrote about the anisotropic V-cycle and is the reason both are worth measuring rather than bounding.

The third time #

An aligned repair failing off its axis is now a pattern on this site rather than an observation.

The depth phase. A line smoother takes the anisotropic V-cycle’s convergence factor from 0.9565 to 0.0370, and semi-coarsening to 0.1074. Turn the same problem’s axes round and both return 0.967 and 0.966 — what the repair encodes is not cope with anisotropy but cope with anisotropy along y, and the direction was supplied by whoever typed the string.

The synthesis phase. Smoothed aggregation returns 0.193 on the aligned problem and 0.789 on the rotated one, and the reason is readable off the stencil before any solver runs.

And this. A scheme exact at 0° and the worst of three at 45°.

The three have nothing in common mathematically. One is a smoother, one is a coarsening, one is a discretisation; two are solvers and one is not. What they share is that each was tuned using a quantity measured along a coordinate direction — a line of the grid, a column of the stencil, a cell Péclet number computed from h — and each therefore encodes an assumption about where the difficulty lies that nothing in the method states.

What this does not say #

Three limits on the claim, each of which would change the numbers.

The comparison is between three discretisations, not between three methods. All three use central differences for the convection term and differ only in what is added to the diffusion, which is the only way to attribute the difference to the addition. A practical upwind scheme differences the convection term itself, and is a different object from the one measured here even though it produces the same added diffusion in one dimension.

The flow is uniform. b is the same vector at every node, so “the angle to the grid” is one number rather than a field. In a real problem the flow turns, the alignment is good in one region and bad in another, and the tuned scheme is exact nowhere and best somewhere — which is a harder measurement and not obviously a worse situation, since the failure at 45° here is a failure over the whole domain at once.

And the exact solution is a function of one variable. That is what makes ground truth available at every angle, and it is also the most favourable possible case for a streamline method: the solution genuinely has no crosswind structure to represent. A scheme that adds nothing across the flow is being tested on a problem that has nothing across the flow, and it still fails at 45°.

That last one is worth sitting with. The failure is not that the crosswind physics was missed. There is no crosswind physics. It is that a discretisation which is correct along one direction and untouched across it produces an operator whose discrete solution has crosswind structure that the continuous one does not.

What a diagonal Péclet number would be #

There is an obvious repair to try, and it is worth saying why it is not the repair.

The coefficient τ is built from a cell Péclet number Pe = |b|h/2ε, and h there is the grid spacing. For a flow at 45° the distance between grid lines along the flow is h√2, so a natural fix is to use the streamwise mesh size instead — a larger h, a larger Pe, a larger τ.

That changes the size of what is added and not its direction. The tensor is still τ b bᵀ, still rank one, still with nothing across the flow, and the crosswind oscillation the next essay measures is untouched by any choice of τ. A coefficient repair cannot fix a rank deficiency, which is the shape of the argument and is why the standard repairs add a second term rather than adjusting the first.

Where the error goes instead #

The next essay is about the mechanism, and one figure from it belongs here because it answers the question this one leaves open: if the scheme is exact along the flow, what is it getting wrong?

The error is entirely across the flow. Along it the tuned coefficient is still doing what it does in one dimension; the direction it does not touch is the one where the discrete operator has been left with a central difference at a cell Péclet number above one, which is exactly the configuration the whole convection field exists to describe as unstable.

So the two-dimensional failure is the one-dimensional failure, in a direction the one-dimensional analysis has no name for. That is not a new mechanism; it is the same mechanism, unrepaired, because the repair was pointed somewhere else.

Why the exactness is a coincidence worth understanding #

It is easy to read the one-dimensional result as a piece of good fortune — a coefficient that happens to make a discretisation exact — and the two-dimensional measurement says what kind of fortune it is.

The one-dimensional problem has an exact solution that is a single exponential, and the tuned difference operator has the same exponential as an exact discrete solution. The coefficient ξ = coth(Pe) − 1/Pe is the number that makes the discrete characteristic root equal the continuous one. That is a matching of two one-dimensional objects, and it is exact because both objects are one-dimensional.

In two dimensions the exact solution along a grid line is still that exponential, so the aligned case inherits the match. Off the axis the discrete operator’s characteristic behaviour is a two-dimensional object — a stencil with five or nine entries and a cross term — and there is no coefficient that makes it match, because there is no single number for a coefficient to be.

That is a stronger statement than the exactness does not extend. It says nothing of that shape can extend: a scalar τ cannot match a two-dimensional operator to a two-dimensional solution except in the cases where one of the two is degenerate. Which is why the practical schemes stopped looking for a better τ and started adding terms.

What the drag does #

The slider is the diffusion, which is the Péclet number in disguise. Two things move with it and one does not.

The cone narrows: from about 35° at Pe = 0.63 to about 10° at 3.13. And the fall at five degrees gets steeper, because the exactness at zero is worth more when the layer is sharper.

What does not move is the exactness itself. At every position of the slider the tuned scheme reproduces the solution at the nodes to the level of rounding when the flow is along a grid line — which is what makes the cliff at the second point a statement about alignment rather than about the layer, and is the assertion the figure carries at every frame.

Exactness that is an artefact of the alignment #

A property that holds exactly in one orientation and fails off it is a property of the description. This site’s phase found the same shape in an orthogonalisation, where the description is the order of the columns.