The smaller cell downstream

A Galerkin reduced model of a linear system $\dot x = Ax$ keeps its poles inside the numerical range of A, so when that range sits in the left half plane every such model is stable, whatever shifts built its basis. Half the conditions and a certificate found that certificate for a convection–diffusion system after a two-sided reduction of the same system had come back unstable in a model that cannot be run. A certificate written in coordinates then showed the certificate belongs to the inner product the state is written in, and the inner product the mesh already computed found the coordinates every code actually uses breaking it: nodal values on a mesh graded towards the outflow boundary, where the layer forms and where anyone resolving it puts the small cells. At Péclet 100 and a grading ratio of a thousand the range’s edge reached +8.94 and seven of nineteen shift placements returned unstable models. Weighting the projection by cell size gave every model back.

Every operator in that essay was centred. A centred difference at a cell Péclet number above two oscillates, and a code that grades its mesh towards a layer very often also upwinds its convection, replacing the centred difference with a one-sided difference taken from the side the flow comes from. A different equation on every grid showed that on a uniform mesh this is exactly the centred stencil for the same equation with diffusion $\varepsilon + h/2$. Extra diffusion is extra dissipation, and dissipation is what the certificate asks for.

So the obvious guess is that upwinding restores the nodal certificate, or at least pushes the crossing to a grading nobody uses. It does neither. It lowers the edge by about a quarter and leaves it positive, and it moves the failures to a different order rather than removing them. The reason is one line of arithmetic about each row of the stencil, and it says exactly where the artificial diffusion is measured and why it is measured against the wrong cell.

The same meshes, upwinded #

The system is the one the previous measurement used: thirty interior nodes on the unit interval, $n + 1$ cells growing geometrically so that the largest is a chosen multiple of the smallest, the small cells at the outflow boundary $x = 1$, velocity one, diffusion $\varepsilon = 1/\text{Pe}$, and the same random input and output vectors. The only change is the convection stencil. Centred, row i carries $(u_{i+1} - u_{i-1})/(h_l + h_r)$; upwind, with the flow towards $x = 1$, it carries $(u_i - u_{i-1})/h_l$, where $h_l$ is the width of the cell on the node’s upstream side.

For each grading ratio from 1 to 1,000 and Péclet numbers 10, 40 and 100, the right edge of the numerical range — the largest eigenvalue of the symmetric part of A — is computed twice: in nodal coordinates, and in the inner product weighted by cell size, $\langle x, y \rangle = x^{\mathsf T} M y$ with $M$ the diagonal of half-sums of neighbouring widths. Then one-sided rational Krylov models of orders 1, 2, 3 and 6 are built at nineteen shift placements, and each is checked for a pole in the right half plane.

At Péclet 100 the centred nodal edge runs $-0.10, 0.48, 1.38, 2.71, 4.69, 6.75, 8.94$ across ratios 1, 3, 10, 30, 100, 300 and 1,000. Upwinded it runs $-0.26, 0.30, 1.05, 2.08, 3.63, 5.19, 6.78$. The upwind line lies under the centred one at every grading, by roughly a quarter of its height at the fine end, and crosses zero at the same place, between a ratio of one and three.

The other two Péclet numbers say the same thing. At Péclet 40 upwinding takes the edge at 1,000 from 4.18 to 3.22 and still crosses zero before a ratio of three. At Péclet 10, where diffusion is strong enough that the centred edge only crossed at a ratio of ten, upwinding moves the crossing to between ten and thirty and ends at 0.80 instead of 1.19. That is the most it moves anything: one step along the grading axis, at the Péclet number where the certificate mattered least.

In the weighted inner product the story runs the other way. The centred weighted edge sits at $-0.10$ at every grading at Péclet 100, because in that inner product the centred convection is exactly skew and the edge is the diffusion’s own. Upwinded, the weighted edge starts at $-0.26$ on the uniform mesh and moves further left as the mesh is graded, to $-0.45$ at a thousand. Upwinding makes the weighted certificate stronger and the nodal one only slightly less wrong.

The convection alone #

To see where the positive part comes from, take the diffusion out. The next figure is the right edge of the convection stencil by itself, in nodal coordinates, on the same meshes.

On the uniform mesh the centred convection is skew, with edge exactly zero, and the upwind convection is dissipative, with edge $-0.159$. That is the uniform-mesh statement that upwinding adds diffusion. On every graded mesh both edges are positive and both grow in proportion to the grading. Centred: 0.86 at a ratio of three, 3.47 at ten, 33.2 at a hundred, 300 at a thousand. Upwind: 0.43, 1.59, 13.1 and 109. Upwinding divides the edge by a factor that creeps from two at a ratio of three to two and three-quarters at a thousand, while the edge itself grows three-hundredfold over the same range: it takes a constant-sized bite out of a quantity whose growth it leaves alone.

In the cell-size-weighted inner product the same two stencils behave as the continuous operator does. The centred one is skew to rounding at every ratio — its symmetric part is of order $10^{-13}$ — and the upwind one is dissipative, with edge from $-0.159$ on the uniform mesh to $-0.220$ at a thousand. The weight turns the convection back into what it discretises, and upwinding then adds exactly what it was supposed to add.

The whole operator’s nodal edge is far smaller than the convection’s alone — 6.78 against 109 at a ratio of a thousand — because the diffusion stencil, which is dissipative in nodal coordinates at every grading measured, absorbs most of it. Its contribution grows like $\varepsilon/h^2$ in a cell of width h while the convection’s grows like $1/h$, so in the smallest cells the diffusion wins, and what is left over is decided in the cells where the two are comparable. That is an account of why the sum is so much smaller than the part, not a measurement of where the edge’s eigenvector lives.

One row of the stencil #

Why is a one-sided difference, which dissipates on a uniform mesh, not dissipative on a graded one? Write out its symmetric part in nodal coordinates.

Row i of the upwind stencil has $-1/h_i$ on the diagonal and $+1/h_i$ just below it, where $h_i$ is the width of the cell upstream of node i. The symmetric part averages the matrix with its transpose. Its diagonal is still $-1/h_i$. Its entry to the left of the diagonal in row i is half of row i’s own coupling, $1/(2h_i)$. Its entry to the right is half of the next row’s coupling back to node i, and the next row belongs to node $i + 1$, whose upstream cell is $h_{i+1}$ — so that entry is $1/(2h_{i+1})$. The row sums to

$$-\frac{1}{h_i} + \frac{1}{2h_i} + \frac{1}{2h_{i+1}} = \frac{1}{2}\left(\frac{1}{h_{i+1}} - \frac{1}{h_i}\right).$$

On a uniform mesh that is zero, and the interior rows contribute nothing positive; the boundary rows, which lack a neighbour, make the whole thing dissipative. On a mesh whose cells shrink in the direction of the flow, $h_{i+1} < h_i$, and every interior row sums to something positive. The dissipation upwinding adds sits on the diagonal and is measured against the node’s own upstream cell. The coupling that the symmetric part hands back is measured partly against the next cell downstream, which is smaller. The smaller cell downstream wins.

At a ratio of a hundred the row sums run from 0.67 near the inflow to 42.5 in the last interior row, and they agree with the formula to rounding. The formula also says what the measurement confirmed about direction in the centred case: with the small cells at the inflow the cells grow downstream, every interior row sum is negative, and the nodal range stays in the left half plane. The grading the physics asks for — small cells where the layer is, at the outflow — is the one that makes each row’s upstream cell larger than its downstream one.

The row sums are an upper bound on the edge rather than the edge: a positive row sum makes the edge positive only together with the rest of the matrix. At a ratio of a hundred the convection’s actual edge is 13.1 while the largest row sum is 42.5. What the row sums give that the edge does not is the mechanism, located node by node, and the fact that no choice of upwind weight measured against the upstream cell alone can cancel it.

Weighted by cell size, row i of the upwind stencil is multiplied by $m_i = (h_i + h_{i+1})/2$, and becomes $(m_i/h_i)(u_i - u_{i-1})$. The factor $m_i/h_i$ is close to one when neighbouring cells are close in size, which is what a geometric grading with thirty nodes makes them, and the row is then close to the uniform upwind row. That is why the weighted edge is negative, and why it grows more negative rather than less as the grading grows.

The failures move rather than go #

A positive edge removes the guarantee; it does not make any model fail. The count of models that actually fail is the measure the previous essay used, and it is the one that matters to someone running the reduced model.

At Péclet 100 the centred operator returned $0, 0, 0, 4, 5, 7, 7$ unstable models across the seven gradings, 23 in all; the upwind operator returns $0, 0, 0, 0, 6, 7, 5$, 18 in all. At Péclet 40 the centred failures were two at a ratio of 300 and four at a thousand; the upwind ones are four at a ratio of a hundred and none at the two finer gradings. At Péclet 10 neither scheme returns any unstable model.

The totals are close, and the pattern is not a shift in one direction. At Péclet 40 upwinding makes the finest gradings safe and a coarser one unsafe. At Péclet 100 it clears a ratio of thirty and makes a ratio of a hundred slightly worse. What changes more than the count is which models fail. Of the 23 centred failures at Péclet 100, 21 are at order 2 and 2 at order 6. Of the 18 upwind failures, 17 are at order 1 and one at order 2. Upwinding moves the failures from one order to another.

And it makes them shallower. The worst centred pole at Péclet 100 is 1.37, at a ratio of 300. The worst upwind pole is 0.46 at the same ratio, 0.28 at a hundred, and 0.06 at a thousand. A reduced model with a pole at +0.06 is unstable, and a simulation long enough to reach $t = 50$ grows by twenty times, but its instability is a small fraction of the system’s own decay rate.

In the weighted inner product nothing fails, in either scheme, at any grading, Péclet number or order.

One order, three projections #

The previous essay drew the rightmost pole at every shift placement for order 2, where the centred failures were. For the upwind operator the failures are at order 1, and the same picture drawn there shows where along the placements they sit.

The centred nodal models at order 2 fail at the low placements, from $s = 0.1$ to about 0.7, with poles up to 1.23. The weighted models on the same basis stay near $-1$ or below everywhere.

At order 1 the centred nodal models never fail: their rightmost pole stays near $-3$ or below. The upwind nodal models fail at the five lowest placements, from $s = 0.1$ to 0.46, with the rightmost pole at 0.28 and falling steadily to zero as the shift grows. The same bases projected in the weighted inner product put the rightmost pole between $-3.4$ and $-4.7$ over that stretch.

The low placements are where the one-sided basis is built from the solves with the smallest shifts, which see the slow part of the dynamics: the transport of mass across the domain and into the layer. That is the part of the operator in which nodal coordinates overstate the state most, because a vector concentrated in the small outflow cells has many nodes and little mass. An order-1 model is a single direction; if that direction leans into the small cells, the nodal Rayleigh quotient is dominated by the positive row sums there.

Upwinding is a different system as well #

There is a second difference between the two operators, and it is not about the projection. A different equation on every grid found that on a uniform mesh the upwind scheme is the exact discretisation of the same problem with diffusion $\varepsilon + h/2$: the scheme solves a different equation, and which one is chosen by the mesh. On a graded mesh the added diffusion is $h_i/2$ in each cell. At a ratio of a thousand the widths run from about $2 \times 10^{-4}$ at the outflow to 0.21 at the inflow, so the added diffusion is ten times the physical $\varepsilon = 0.01$ upstream and a hundredth of it inside the layer.

That is precisely where a code wants it: large where the cells are too coarse to resolve anything, and negligible in the layer the mesh was graded to resolve. A parameter that is also a price measured what the fitted fraction of that extra diffusion costs on a problem with no layer at all. The full system’s slowest decay rate moves with it — at Péclet 100 from $-19.2$ centred to $-10.9$ upwinded on the uniform mesh, and from $-3.75$ to $-4.15$ at a grading of a thousand — so the reduced models of the two schemes are approximating different transfer functions. The count of failures above compares each scheme’s reduced models against its own full system, which is the comparison a practitioner makes.

It is also why upwinding’s effect on the nodal edge is smallest exactly where it would be needed. The added diffusion is proportional to the cell width, and the positive row sums are largest where the widths are smallest. The dissipation arrives where the cells are large; the non-normality lives where they are small.

What the nodal edge is measuring, again #

The previous essay read the nodal edge as the rate at which the Euclidean norm of the nodal vector can grow — a norm that counts nodes rather than mass, and so weighs the many small cells at the outflow as heavily as the few large ones upstream. A spectral radius that grows first and the eigenvalues that are not there described the same transient growth from the side of the powers and of the resolvent.

The row sums above give that reading its arithmetic. Upwinding transports each node’s value into the next node downstream, and on a mesh where the next cell is smaller, the same mass arriving there is a larger nodal value. The Euclidean norm of the nodal vector counts the larger value and not the smaller mass, and sees growth. Any stencil that moves mass downstream into shrinking cells will show it, however much dissipation it carries on its own diagonal, unless that dissipation is measured against the smaller cell — which is what weighting by $m_i$ does.

What a code should take from it #

Upwinding is not a stability fix for the reduced model; it is a fix for oscillations in the full one. A code that upwinds a graded convection–diffusion operator and reduces it by Galerkin projection in nodal coordinates has the same problem the centred code had, at almost the same gradings, with failures of a similar number moved to a lower order. The lower order is the worse place for them in one respect: an order-1 model is the cheapest one to try, and the first one a code tests a pipeline with.

The weight is still the fix, and upwinding makes it stronger. The mass matrix of a finite-element discretisation, or the diagonal of cell sizes of a finite-difference one, is already computed, and projecting in it gives every one-sided model of either scheme a certificate at every grading measured, with a margin that grows as the grading grows.

What this rests on #

Thirty interior nodes, geometric grading from 1 to 1,000 towards the outflow, Péclet numbers 10, 40 and 100, orders 1, 2, 3 and 6, nineteen shift placements from 0.1 to 100 at each order, with shifts placed geometrically from s to 4s; the first-order upwind stencil and no other. A second-order upwind stencil, or a streamline-diffusion finite-element method, adds dissipation of a different form and is not measured. The row-sum formula is exact for the first-order stencil on any mesh; the edges and failure counts are measurements on these meshes.

The claim that has to fail #

The claim is that upwinding, by adding diffusion, makes the nodal operator dissipative on any mesh. At Péclet 100 on the mesh graded a thousand to one the upwind operator’s nodal edge is 6.78 and five of its 76 reduced models are unstable. The refusal is fed the claim that the nodal edge is negative, and fails.

Still open: a stencil weighted against the downstream cell, and the two-sided model #

An upwind stencil that measures against the smaller cell. The positive row sum is the difference between the downstream cell’s inverse width and the upstream one’s. An upwind stencil that divided by the downstream width, or by $m_i$, would make each row’s symmetric part non-positive in nodal coordinates; it would also be a different discretisation, with a different truncation error at the outflow. Whether any consistent first-order stencil can be dissipative in nodal coordinates on a mesh graded towards the outflow, and what it costs in the layer, is the direct next measurement.

The two-sided model. The upwind two-sided models at order 1 fail at the highest placement, $s = 100$, at a ratio of a hundred, as the centred ones did at order 2. The one-sided construction has a certificate in the weighted inner product for both schemes; whether a two-sided projection whose left space is chosen in that inner product inherits any of it is still the open question the previous essay left.

Where the edge’s eigenvector lives. The account above places the nodal edge in the cells where diffusion and convection are comparable. Drawing the eigenvector of the symmetric part that attains the edge against the mesh, for both schemes, would say whether that is right.