Docs: the integration-point history's low Courant number limit, and the choice between transport managers

docs/advanced/eulerian-advection-diffusion.md

@@ -80,6 +80,50 @@ is resolved in time (a fraction of a feature width per step), the Eulerian solve
 is cheaper and more accurate; if the step is deliberately long relative to the
 transported features, the semi-Lagrangian solver is the one that survives it.
 
+## The three transport managers
+
+`DuDt` selects the transport, and three managers are worth considering for a
+scalar field. The choice turns on the Courant number the model runs at and on
+whether the model is carrying particles for another reason.
+
+**Eulerian SUPG** (`uw.systems.ddt.EulerianSUPG`, the default) is the general
+choice. Its error falls as $\Delta t^2$, it puts no lower limit on the Courant
+number, and it conserves the integral of the transported field to solver
+tolerance. On the LeVeque deformation test at its standard period it matches the
+integration-point history's accuracy at a third of the cost, and holds the
+enclosed volume to 4e-5 against that scheme's 5e-3. Use it unless something
+below applies.
+
+**Semi-Lagrangian on the integration points**
+(`uw.systems.ddt.IntegrationPointSemiLagrangian`) is the accurate choice at
+larger Courant numbers. Its error is flat between Courant 0.5 and 2, so a model
+that takes long steps keeps its accuracy where the Eulerian scheme loses it, and
+it loses 45 times less of the second moment than the nodal scheme does. It
+carries a history at the integration points rather than at nodes, which is the
+reason to prefer it for tensor transport, where the extra sub-cell resolution
+has more to represent. It has a low Courant number limit: the fit it performs is
+not contractive under pure advection, and below about Courant 0.5 a mode grows.
+The growth is suppressed by physical diffusion and is unreachable when the
+timestep comes from the Courant condition, but it is a real limit for
+advection-dominated flow with sharp interfaces. Adding diffusivity to damp it
+makes the answer worse at every strength, so the limit is a reason to choose a
+different manager rather than something to correct. The measurements are in
+`docs/developer/subsystems/integration-point-variables.md`.
+
+**Lagrangian on a swarm** (`uw.systems.ddt.Lagrangian_Swarm`) transports the
+field on particles, and is **in development**: it is being measured as part of
+the particle-in-cell work and its accuracy and cost are not yet settled, so
+treat what follows as the intended role rather than a recommendation. It is
+aimed at models that already carry a swarm for material tracking, where the
+transport rides on particles the model is advecting anyway. We would not expect
+to introduce particles in order to use it.
+
+**Semi-Lagrangian at the nodes** (`uw.systems.ddt.SemiLagrangian`) remains the
+historical default of `AdvDiffusionSLCN`. It re-interpolates once per step, which
+costs it accuracy at small Courant numbers, and on a deforming flow with a sharp
+interface it diverges below a Courant number that depends on the problem. Prefer
+one of the three above.
+
 ## Choosing the time scheme
 
 Measured on a rotating Gaussian, one revolution, relative $L_2$ error; the full

docs/developer/subsystems/integration-point-variables.md

@@ -261,6 +261,64 @@ The trace-back samples twelve points per cell rather than the P2 nodes, and
 the snapshot evaluation at the moving feet misses the locator cache every
 step, so the update costs about three times the nodal one.
 
+### The low Courant number limit in practice
+
+The growth measured above sets a lower limit on the Courant number at which
+this history is useful, and the limit matters most for advection-dominated
+flow carrying a sharp interface. Two measurements bound it.
+
+On the rotating Gaussian at Courant 0.1 with the 12-point rule, pure
+advection, one revolution (503 steps), the energy finishes 0.92 % above the
+analytic value with a -3.7 % undershoot. Adding physical diffusion removes
+it: at cell diffusion number `D = κΔt/h² = 1e-3` the undershoot is -0.04 %
+and at `D = 1e-2` the energy ratio is 0.9934, below one. The two dimensionless
+numbers are not independent within a run, since `Pe = C/D` and `D` carries no
+velocity, so a cell at low Courant number sits at proportionally low Péclet
+number. A timestep taken from the Courant condition fixes
+`D = C_max/Pe_max` and cannot reach the growth region; a timestep imposed by
+something else, such as a coupled free surface or an elastic timescale, can.
+
+On the LeVeque single-vortex deformation at its standard period `T = 8`
+(unit box, `h = 1/48`, P2, tanh level set of one cell, no diffusion, no
+reinitialisation) the filament thins below the cell and the growth becomes
+visible. At Courant 2 the returned field has a relative L2 error of 1.54e-1,
+a peak of 1.080 and an undershoot of -0.068, with ringing behind the returned
+blob. Eulerian SUPG at Crank-Nicolson reaches 1.50e-1 on the same problem and
+holds the enclosed volume to 4e-5 against this scheme's 5e-3. Both are
+approaching the error the mesh can resolve at that spacing, about 1.45e-1, so
+on a sharp interface at this resolution the two schemes are equivalent in
+accuracy and SUPG is the cheaper of them.
+
+### Added diffusion does not remove the growth
+
+A stabilising diffusivity `κ = ε h²/(π² Δt)`, a cell diffusion number
+`D = ε/π²`, has been suggested for the growth. We measured it on the case
+that shows the overshoot, LeVeque `T = 8` with this history at Courant 2:
+
+| ε | D | cell Péclet | L2 rel | peak | undershoot | band mass |
+|---|---|---|---|---|---|---|
+| 0 (control) | — | — | 1.543e-1 | 1.0797 | -6.77e-2 | 1.168 |
+| 1e-3 | 1.0e-4 | 19700 | 1.684e-1 | 1.0522 | -5.35e-2 | 1.275 |
+| 1e-2 | 1.0e-3 | 1970 | 3.190e-1 | 0.8804 | -2.85e-3 | 1.896 |
+| 1e-1 | 1.0e-2 | 197 | 7.053e-1 | 0.4275 | -1.59e-7 | 3.002 |
+| 1 | 1.0e-1 | 19.7 | 9.000e-1 | 0.1619 | 0 | 3.120 |
+
+The L2 error rises with `ε` throughout, so the unstabilised run is the most
+accurate of them and the damping does not pay at any strength. At `ε = 1e-3`
+the overshoot moves a third of the way while the error and the interface
+width each grow by 9 %. There is no intermediate setting: between `1e-3` and
+`1e-2` the peak goes from 1.052 to 0.880, from overshooting to over-damped,
+and by `ε = 0.1` the field is bounded to 1.6e-7 with a peak of 0.43 and an
+interface three times its original width. The damping is applied to the whole
+field and so reaches the physically sharp interface, which is most of what
+this problem contains, more strongly than it reaches the spurious mode.
+
+The limit belongs in the choice of scheme rather than in a correction term.
+Use this history at Courant numbers of about 0.5 and above, where its error is
+flat and its conservation is good, and use the Eulerian SUPG solver for
+advection-dominated flow with sharp interfaces at small Courant number. See
+`docs/advanced/eulerian-advection-diffusion.md` for the comparison.
+
 ## Swarm proxy at the integration points
 
 A swarm variable normally reaches the weak form through a nodal proxy: