From 761a708b6dbb601c2ef2399d232f9d2ee96b9aef Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 9 Sep 2026 08:47:21 -0700 Subject: [PATCH 1/2] Docs: the integration-point history's low Courant number limit, and the choice between transport managers The growth measured in the oversampling section sets a lower limit on the Courant number at which the integration-point history is useful, and nothing said so. Two measurements bound it: on the rotating Gaussian at Courant 0.1 under pure advection the energy finishes 0.92 % high with a -3.7 % undershoot, and the growth is gone by cell diffusion number 1e-2; on the LeVeque deformation test at T = 8, where the filament thins below the cell, the returned field overshoots by 8 % and Eulerian SUPG reaches the same L2 error more cheaply and with 100 times less volume drift. Both schemes are near the error the mesh can resolve at that spacing, so on a sharp interface they are equivalent in accuracy. Also records that the stabilising diffusivity kappa = eps h^2/(pi^2 dt) suggested for the growth does not work: the L2 error rises with eps throughout, the proposed eps = 1e-3 moves the overshoot a third of the way while costing 9 % of accuracy and 9 % of interface width, and there is no intermediate setting, since the peak goes from 1.052 to 0.880 between eps = 1e-3 and 1e-2. The damping reaches the physically sharp interface more strongly than the spurious mode. docs/advanced/eulerian-advection-diffusion.md gains a section naming the three managers worth considering and what each is for: Eulerian SUPG as the general choice, the integration-point history for larger Courant numbers and for tensor transport, and the swarm history when the model already carries particles for material tracking. The nodal semi-Lagrangian history is described as superseded; it diverges on a deforming flow below a Courant number that depends on the problem. Measurements and scripts: ~/+Simulations/slcn-ip-real-cases/ Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv --- docs/advanced/eulerian-advection-diffusion.md | 41 +++++++++++++ .../subsystems/integration-point-variables.md | 58 +++++++++++++++++++ 2 files changed, 99 insertions(+) diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md index 2b7a1baf..0dbad7da 100644 --- a/docs/advanced/eulerian-advection-diffusion.md +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -80,6 +80,47 @@ 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. It is worth using when the model already carries a swarm for +material tracking, so the transport rides on particles it is advecting anyway. +We would not 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 diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index d8ee4cfc..36d9ec39 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/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: From b605fa7f3b73bac389a43f50b96554d5f6990cfa Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 9 Sep 2026 09:10:34 -0700 Subject: [PATCH 2/2] Docs: flag the swarm transport manager as in development Its accuracy and cost are still being measured as part of the particle-in-cell work, so the entry states the intended role rather than making a recommendation. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv --- docs/advanced/eulerian-advection-diffusion.md | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md index 0dbad7da..a4556d63 100644 --- a/docs/advanced/eulerian-advection-diffusion.md +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -111,9 +111,12 @@ 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. It is worth using when the model already carries a swarm for -material tracking, so the transport rides on particles it is advecting anyway. -We would not introduce particles in order to use it. +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