A smoother that stops being one

Two constants hold up this site’s whole multigrid field. Weighted Jacobi at ω = 2/3 multiplies the oscillatory half of the error by exactly one third, and no other weight does better. Everything built on top of them inherits them: the smoothing factor whose derivation is the field’s first rung, the V-cycle rate that does not notice the grid size, and the two-grid correction whose coarse operator is the same discretisation one level down.

Both constants are properties of symmetry rather than of Jacobi. Put a first derivative into the operator and the Fourier symbol goes complex, the closed form stops applying, and what damps is a modulus. The essay that found that measured it at one diffusion coefficient and reported one number: the best any weight achieved was 0.84, against exactly ⅓.

One number is not a reading. Swept across the cell Péclet number, the failure turns out to be continuous, to be a function of exactly one variable, and to have three landmarks in it — one of which sits on the threshold that field already found by an entirely different route.

The factor depends on one number, and it is not the diffusion coefficient #

The symbol of one weighted Jacobi sweep is the diagonal-scaled operator subtracted from one, and on this operator it can be written out. The row is (−ε/h² − 1/2h, 2ε/h², −ε/h² + 1/2h), the diagonal is d = 2ε/h², and at frequency θ the sweep multiplies a mode by

g(θ)  =  1 − ω(1 − cos θ)  −  i ω Pe sin θ,        Pe = h/2ε

The real part has no ε and no h left in it. Both off-diagonals carry −ε/h², their sum against the diagonal is exactly −1, and everything about the mesh and the diffusion cancels. What survives is in the imaginary part, and it survives as a single dimensionless group: the two off-diagonals differ by 1/h, that difference divided by d is h/2ε, and that is the cell Péclet number.

So the damping is a function of ω and Pe alone. Checked against the library rather than left as algebra, the imaginary part of the symbol agrees with −ω·Pe·sin θ to twelve digits at every frequency tried, and the real part with 1 − ω(1 − cos θ) to the same.

The consequence is measurable and it is the reason the sweep below is a sweep in one variable. Holding Pe fixed and changing the grid from 15 points to 255 — a factor of sixteen in h, with ε moved to match — leaves the smoothing factor at ω = 2/3 identical to eight decimal places at all five sizes: 0.34923030 at Pe = 0.15625, 0.74535599 at Pe = 1, 5.21898910 at Pe = 7.8125. A quantity that does not move when the grid moves by a factor of sixteen is not a property of the grid.

The closed form the field lost comes back, with Pe in it #

The symmetric case has a closed form — max(|1 − ω|, |1 − 2ω|) — and losing it is what forced the convective case to be scanned. It can be recovered. Writing s = 1 − cos θ, the squared modulus is

|g|²  =  1 − 2ωs + ω²s²(1 − Pe²) + 2ω²Pe² s

a quadratic in s, and the oscillatory half of the spectrum is exactly s ∈ [1, 2]. For Pe < 1 the quadratic is convex, so its maximum is at an endpoint; for Pe > 1 it is concave, and its vertex lies below s = 1 for every ω ≤ 1, so the maximum is at an endpoint there too. The two endpoints evaluate to (1 − ω)² + ω²Pe² at θ = π/2 and (1 − 2ω)² at θ = π, which gives

μ(ω, Pe)  =  max( √( (1 − ω)² + ω²Pe² ),  |1 − 2ω| )

At Pe = 0 the first branch is |1 − ω| and this is the symmetric closed form, unchanged. The convection enters in exactly one place: it adds ω²Pe² under the square root of the branch that governs the middle of the oscillatory band, and it leaves the branch at θ = π alone.

That second branch is worth pausing on, because it says where the damage is not. At θ = π a mode alternates from one grid point to the next, sin θ is zero, and the imaginary part of the symbol vanishes at every Péclet number — so the most oscillatory mode the grid carries is damped by exactly |1 − 2ω| whatever the convection does. At ω = 2/3 the library returns 0.333333333333 for that mode at Pe = 0.016 and 0.333333333333 at Pe = 7.81, where the smoothing factor over the band is 5.22. The mode a smoother is most obviously built to remove is the one convection leaves untouched.

So the failure sits at the edge of the band rather than at the top of it. The worst mode is at θ = π/2 at every Péclet number above zero, which the scan confirms by reporting its argument as well as its value, and it is there because sin θ is at its largest and 1 − cos θ at its smallest in the same place. That mode is the boundary between the half of the spectrum the smoother is responsible for and the half the coarser grid is supposed to take, and it is precisely the mode whose damping the division of labour depends on.

Held against the library’s scan over 0 < ω ≤ 1 at fifteen Péclet numbers from 0.01 to 20, the largest disagreement is 3.6·10⁻¹⁵, which is rounding. That is two routes to a number applied to a quantity this field had been computing only by scanning, and the second route is worth having for a reason beyond confirmation: a scan reports the smallest value on its own grid, which is an upper bound on the true minimum and not the minimum. The figures’ 41-point weight scan puts the best at ω = 0.654 and reads 0.3464 at Pe = 0.0156; the closed form puts it at 0.66672 and reads 0.33344. The gap is the grid, not the function.

At the field’s own weight the formula collapses to something short. Setting ω = 2/3 makes the second branch exactly ⅓ and the first branch always at least that, so

μ(2/3, Pe)  =  ⅓ √(1 + 4Pe²)

which reproduces 0.333496, 0.349230, 0.618368, 0.745356, 1.093700 and 5.218989 at the six Péclet numbers drawn here, to every digit the scan reports.

Three landmarks, and each one is exact #

Minimising the closed form over ω gives the best a weighted sweep can do, and it does so in two regimes with a clean boundary between them.

Below a threshold the two branches cross before either reaches its own minimum, so the optimum sits at the crossing: ω* = 2/(3 − Pe²) and μ* = (1 + Pe²)/(3 − Pe²). At Pe = 0 that is ω = 2/3 and μ = ⅓, which is the constant the field was built on, arriving as the zero-convection end of a formula rather than as a special case.

Above it the first branch’s own minimum is already the larger of the two, and the optimum is ω* = 1/(1 + Pe²) with

μ*(Pe)  =  Pe / √(1 + Pe²)

The boundary between the regimes is at Pe = 1/√3, and it is the first landmark: there the optimal weight is exactly 3/4 and the best smoothing factor is exactly 1/2, measured at 0.500000000 against the scan. The second is Pe = 1, where ω* = 1/2 and μ* = 1/√2 = 0.70710678 — attained at 0.500000 and 0.707106781 at n = 15, 31, 63, 127 and 255, identically. The third is where the field’s own weight stops contracting at all, and it comes later.

Written as an angle the two regimes are one statement. Pe is the ratio of the convective half of the off-diagonal, 1/2h, to its diffusive half, ε/h², so it is the tangent of an angle the stencil makes, and past 30° the best achievable smoothing factor is the sine of that angle. At 30° the sine is 1/2; at 45° it is 1/√2. The classical third is what the sine formula has not taken over yet.

Read that way the convection-dominated end is not a cliff. The sine of an angle approaching a right angle approaches one from below, so the best achievable smoothing factor rises towards one and never passes it — 0.9488 at Pe = 3, 0.9919 at Pe = 7.8, 0.9965 at Pe = 12. What collapses past the threshold is not that bound but the weight that attains it, which is a distinction worth keeping separate and is taken up further down.

The two thresholds are one threshold #

Pe = 1 is not a new landmark on this operator. It is the one the convection essay measured by bisection, where the upper off-diagonal −ε/h² + 1/2h changes sign, the matrix stops being an M-matrix, and the computed solution starts leaving the interval the equation guarantees. That was a statement about the discretisation, arrived at without any solver in it.

The smoothing landmark is a statement about relaxation, arrived at by minimising a modulus over a weight. The two coincide, and the angle says why: the sign of the off-diagonal and the size of the imaginary part are the same ratio read twice. While the diffusive half dominates, the off-diagonal is negative and the symbol’s imaginary part is smaller than its real deficit; past 45° both statements reverse together.

At exactly Pe = 1 the arithmetic is unusually clean, and it is worth taking literally. With n = 63 and ε = 1/128 the upper off-diagonal is exactly zero in binary — −ε/h² and 1/2h are both 32 — so the matrix is lower bidiagonal. A single forward sweep solves any system it appears in to a relative residual of 0, exactly, since forward substitution is what a Gauss–Seidel sweep on a lower triangular matrix is. And the boundary-layer problem’s right-hand side, which is where the u(1) = 1 condition enters, is identically zero at that ε: the stencil has stopped looking downstream, the only boundary that carries information cannot reach the interior, and the discrete solution is zero at all 63 points against a true value of 0.135 at the last one.

So the threshold has three independent signatures on the same operator at the same setting — a sign change, a smoothing factor of 1/√2, and a stencil that has gone one-sided. None of them was derived from the others.

That figure is where the 1/√2 comes from. At Pe = 1 and ω = 1/2 the worst oscillatory mode sits at θ = π/2, where the real part of the symbol is 1 − ω = 0.500 and the imaginary part is ωPe = 0.500. The convective smoothing factor is the hypotenuse of the symmetric one and the imaginary part, and at this Péclet number the two legs are equal. 1/√2 is not an approximation of anything; it is what a right angle with equal legs gives, and the equal legs are the definition of Pe = 1.

The cost of that is easy to state and it is the reason the landmark matters. A factor of ⅓ a sweep clears an order of magnitude of oscillatory error in 2.1 sweeps. A factor of 1/√2 needs 6.6 — and the smoother is the expensive half of a V-cycle, so a three-fold increase in sweeps is close to a three-fold increase in the cycle.

Past √2, the field’s own weight amplifies #

The third landmark is where the constant the field ships stops being a contraction. Setting ⅓√(1 + 4Pe²) = 1 gives 4Pe² = 8, so

Pe  =  √2  =  1.414213562

and bisecting on Pe with the library’s scan finds the crossing at 1.414213562373, which is √2 to twelve digits. Above it, a sweep at ω = 2/3 multiplies the worst oscillatory mode by more than one.

Those twelve digits are a check on the whole sweep rather than a decoration on one number. The root is algebra — ⅓√(1 + 4Pe²) = 1 has the single positive solution √2, obtained without running anything — and the crossing is a bisection on a scanned maximum over 20,001 frequencies at 63 grid points, which knows nothing about the algebra. Where they agree, the curve every figure below is a point on can be read as a measurement rather than as a formula plotted through its own predictions.

That is a stronger failure than a slow smoother and it is worth naming precisely. A V-cycle assumes the smoother removes the modes the coarse grid cannot hold and hands the rest down, which is the division of labour the coarse-grid correction depends on. At Pe = 1.56 the smoother returns those modes larger than it found them, so the correction is applied to an error that the sweep has grown. Nothing in the cycle notices: the residual is still computed correctly, the coarse solve is still exact, and the rate is simply wrong.

This is the setting the site had already looked at once, as a closing remark in the essay on the tuned discretisation, and the one number it reported — 0.8425 — is the best-over-ω entry above. Read alone it says a smoother is poor. Read as one point on a curve it says something else: the sweep at the field’s own weight has already crossed one, and the 0.84 is what a weight chosen for this problem rescues.

Where the best over every weight stops meaning anything #

The remaining question is whether choosing ω per problem is a repair, and the closed form answers it in a way a scan cannot.

Pe/√(1 + Pe²) is below one at every Péclet number, so on the face of it some weight is always a contraction. The weight that attains it is 1/(1 + Pe²), which falls like Pe⁻²: 0.291 at Pe = 1.56, 0.0615 at Pe = 3.91, 0.0161 at Pe = 7.81, 0.0069 at Pe = 12. A sweep at ω = 0.0069 changes the iterate by seven parts in a thousand of what a full sweep would, and its smoothing factor is 0.9965. The bound is met by doing almost nothing, and almost nothing is not a smoother.

Restricted to weights of 0.05 and above — the range the figures scan, and a generous floor for anything anyone would call relaxation — the best crosses one at Pe = √39 = 6.245, measured by bisection at 6.244998. Above that there is no usable weight at all, which is why the badge on the last figure reports 1.0272 where the unrestricted formula reports 0.9919. The two numbers are the same finding stated against different floors.

What follows for the field, case by case #

A smoothing factor is not a property of a smoother. It is a property of a smoother and an operator, and the ⅓ that appears throughout this collection carries an unstated hypothesis about the second. The same lesson arrived once before from a different direction, when a strong direction in the Laplacian took the factor to 0.9993 with every component of the method unchanged. This is that finding on a non-symmetric operator, with the difference that the factor here goes above one rather than close to it.

The repair is the same shape as the anisotropy repairs, and so is its defect. Solving whole grid lines exactly rather than sweeping over them restores the one-dimensional third on an anisotropic problem; a smoother built to follow the flow does the corresponding thing here. Both encode a direction somebody had to know, and both return a different number when the direction changes — which is exactly what the tuned scheme in two dimensions measures when the flow is five degrees off the mesh.

Checking Pe is checking two things at once, and it costs nothing. It is h/2ε and needs no arithmetic. Below one, the discretisation obeys a maximum principle and a weighted sweep is a contraction at every weight worth using. Above √2 the field’s own weight amplifies. A code that computes Pe before believing its answer is also, without doing anything else, deciding whether its V-cycle has a smoother in it.

And refining the mesh repairs both. Pe falls with h, so the same refinement that removes the oscillation moves the smoothing factor back down the sine curve. That is the repair which is always correct and often unaffordable, and its cost is now readable as one number rather than two.

The coarse grids are not exempt. A V-cycle halves h at every level, so Pe doubles going down, and a hierarchy whose finest grid sits comfortably at Pe = 0.4 has a level at 1.6 and a level below it at 3.2. The smoothing factor at those levels is 0.43, 1.12 and 2.16 at ω = 2/3. The coarse problem is a different problem from the one the recursion assumed in two dimensions for a structural reason; here it is a different problem in one dimension, for an arithmetic one, and the direction of the change is always the wrong way.

A weight tuned per level is what the closed form actually recommends, and it is cheap, since 1/(1 + Pe²) needs no scan. It is also the pattern an optimum that turned out not to be a choice describes: a formula somebody wrote down for one situation, found to be the minimiser of something once the situation is stated properly. The 2/3 is the Pe = 0 member of that family and was never a constant.

And none of it is visible in the spectrum. The operator here is strongly non-symmetric — 1.34 and 1.95 in relative norm at the last two settings — and a spectrum on such a matrix predicts nothing about a Krylov method. The Fourier symbol used throughout this essay is not an eigenvalue and never was; it is what one sweep does to one grid mode, which is a statement about the sweep rather than about the operator’s eigenvectors. That the analysis survives the loss of symmetry at all is because it never needed the modes to be eigenvectors — only the sweep to act on them one at a time, which a constant-coefficient stencil does.

What is not measured here #

The smoothing factor is half of a two-grid argument and this essay measures only that half. What the coarse-grid correction does to a non-symmetric operator — whether the Galerkin product is still the coarse discretisation, whether interpolation built for a symmetric problem still spans what the smoother leaves behind — is a separate measurement, and the one-dimensional identity that held entry for entry on the Laplacian has no reason to survive an operator with a first derivative in it.

The other omission is deliberate. Every number here is for the central-difference operator, which past Pe = 1 is solving a problem whose answers leave the interval the equation permits. The upwind operator is an M-matrix at every Pe, so its smoother behaves — and it is an exact discretisation of a different equation, with a boundary layer (1 + Pe) times too thick. Measuring a smoothing factor on it would report a healthy number about a problem nobody posed, which is the trade this whole subject keeps making and the one worth stating plainly rather than burying in a rate.