Integration-point variables and a semi-Lagrangian history at the integration points

[lmoresi]

Summary

A field stored at the mesh integration points, and a semi-Lagrangian history built on it.

uw.discretisation.IntegrationPointVariable is a peer of MeshVariable and swarm variables: one value per quadrature point per cell, on an element whose basis is the identity on the mesh rule (mesh.integration_rule). The assembler reads its values at the integration points with no interpolation, so it carries values that are injected into the weak form: a semi-Lagrangian history sampled at the departure points of the integration points, or a material property reconstructed from a swarm with sub-cell resolution.

• The element is built on a UW3-registered PetscSpace type (uwdelta, PetscSpaceRegister), so it needs no PETSc patch and works on stock conda PETSc. PETSc's own PETSCSPACEPOINT cannot be tabulated anywhere but its own points in its own order, which breaks PetscFESetUp, face tabulation and boundary integrals.
• evaluate() is defined as the nearest integration point of the owning cell, the one extension under which a query agrees with what the assembler used at that point. Exact at the variable's own coordinates.
• Guards: the JIT refuses a derivative of the symbol (the tabulated gradient is zero, so it would be a silent zero); solvers check their element's rule against the mesh rule when they attach the auxiliary vector (off-rule the field reads as zero).
• Scalar components only for now.

uw.systems.ddt.IntegrationPointSemiLagrangian drops into AdvDiffusion as DuDt. Its slots are integration-point variables; the value the weak form sees at each integration point is the solution from k+1 steps ago evaluated exactly at that point's departure point. A delta field cannot be sampled off its points, so the slot-to-slot chain is replaced by nodal snapshots of the solution and cached velocities at the last order times, with slot k filled by tracing k+1 RK2 segments back from every integration point.

Both semi-Lagrangian schemes now take the RK2 mid-point velocity at the mid time, 1.5 v^n - 0.5 v^{n-1}; with v^n alone the trace is first order in an unsteady flow. The previous velocity is V_fn evaluated at the true nodes and cached per history level, so any expression (-v, v/2, c(t) v, v - v_mesh) is carried as it was at that time. On by default (midtime_velocity=False reproduces the earlier trace, for controls).

Measurements

Rotating Gaussian, P2, cellSize=0.05, half a revolution, L2 error, nodal SLCN vs integration-point (12 points per cell):

Courant	nodal	integration-point	peak nodal / IP
0.25	1.30e-2	6.3e-4	0.935 / 0.994
0.5	3.55e-3	9.0e-4	0.974 / 0.992
1	3.62e-3	3.24e-3	0.987 / 0.997
2	1.33e-2	1.33e-2	0.996 / 1.000

Energy (∫T²) over 63 steps at Courant 0.5: nodal −1.4 %, integration-point −0.03 %.

Blankenbach 1a (box, Ra 1e4, P2, cellSize=1/32, Courant 1, t = 0.3): interior transport and Vrms identical across the schemes to four digits; wall Nusselt numbers within the first-cell band (a 5e-4 temperature difference in the wall cells). Not a measured gain at cell Péclet of order one, as expected.

The rule must oversample the history space

The Galerkin step with a sampled load is algebraically the weighted least-squares fit of the samples on the rule's points; the shifted field is piecewise polynomial on the shifted mesh and has sub-cell kinks, and the fit contracts only in the sampled norm. One-step growth factors by power iteration at Courant 0.25 (P2 on triangles): 6 points per cell (qdegree=2) 1.028, 9 points (conical rule) 1.005, 12 points (qdegree=3) 1.0003; with diffusion at cell Péclet 100: 1.0095, 0.996, 0.996. Nodal SLCN: 0.9989. So the scheme does not dissipate and relies on oversampling or physics; the constructor raises when the rule has no more points per cell than the history space has local dofs and warns below 2x. P1 histories are 2x oversampled at qdegree=2; P2 needs qdegree=3 (or the conical rule). Raising the rule costs every solver on the mesh its assembly time.

Not in this PR

Vector/tensor histories (the viscoelastic stress is the intended next case), ALE and old-frame trace-back for the integration-point scheme, checkpoint state for it, the conical rule as a mesh option, a swarm-to-integration-point reconstruction.

Tests

tests/test_0064_quadrature_point_fe.py (element: identity tabulation on triangle/tet/quad/hex, cell-only layout, zeros off-rule), tests/test_0065_integration_point_variable.py (assembler reads the stored values: integral equals the hand quadrature sum with a one-point perturbation control; P2 projection exact; evaluate exact at own points and nearest-in-cell elsewhere; both guards), tests/test_0066_integration_point_slcn.py (exact departure-point values for one and two segments; under-sampled rule refused; rotating Gaussian beats nodal; mid-time velocity exact for v, -v, v/2, c(t) v with the v^n-only foot as control). Full suite: 1695 passed. Docs: docs/developer/subsystems/integration-point-variables.md.

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G

[lmoresi]

Adversarial review

We went looking for the ways this can be wrong before it is merged. Findings, with what we did about each.

1. The integration-point scheme is not dissipative; its stability comes from the rule, not the transport step. Power iteration of the one-step operator at Courant 0.25 gives 1.028 at 6 points per cell, 1.005 at 9, 1.0003 at 12 under pure advection (nodal SLCN: 0.9989). A run at 6 points with a P2 history blew up after 75 flat steps. Mitigation: the constructor raises at ≤1x oversampling and warns below 2x; the docs state the mechanism (least-squares fit on the sampled norm, aliasing of the shifted field's sub-cell kinks). Residual risk: 9 points is slowly unstable without diffusion (invisible over one revolution, not over ten); with cell Péclet ≤ 100 it is stable (0.996). Users on tets (14 points, 1.4x for P2) get a warning, not a measurement.

2. The wall Nusselt number is not a benchmark of the transport scheme at this resolution. A 0.9 % wall shift that appeared during the work was a 5e-4 first-cell temperature difference, and it was a bug (the previous velocity evaluated at nudged nodes), found by restarting each variant from the other's state. Interior transport and Vrms never moved. The README of the study records the discriminators. The wall values in the PR are reported within their band, not as a gain.

3. evaluate on an integration-point variable overrides the FE interpolation column-wise after the fact. The delta basis tabulates to zero off the rule, so without the override a query would silently return zero. The override uses the cached owning cells when the evaluator has them and the robust locator otherwise; exterior points fall back to the nearest point on the rank. Risk: a code path that reads a delta field through DMInterpolation directly, bypassing uw.function.evaluate, gets zeros. We found none in the tree; the docs say so.

4. Boundary integrals see zeros from an integration-point field. Face tabulation of the delta space is zero by construction. Correct for a history term, wrong for anything a user puts in a Neumann or Nitsche form. Documented; not guarded, since the assembler evaluates all auxiliary fields on faces regardless of use.

5. The mid-time velocity changes every existing SLCN result at the fourth digit. Default on, by the ruling that SLCN defaults move toward accuracy; midtime_velocity=False reproduces the earlier trace and was verified bit-for-bit against the published Blankenbach 1a row.

6. Parallel. The variable is cell-local by construction and evaluate uses the local owning-cell locator; the new tests were run under two ranks. The semi-Lagrangian history's trace goes through global_evaluate as the nodal scheme's does. No partition-count study of the answers beyond that.

7. Scope limits stated in code, not only in docs: scalar components raise NotImplementedError; the integration-point history has no checkpoint state and no ALE support, and does not register as a state bearer.

Underworld development team with AI support from Claude Code

[Copilot]

🟡 Changes recommended

There are correctness/robustness issues to address (notably a mutable default argument in the new IntegrationPointSemiLagrangian API and a misleading configuration error message).

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR introduces integration-point (quadrature-point) stored fields as a first-class discretisation primitive, and builds a new semi-Lagrangian time-history scheme that samples and stores the advected history directly at those integration points (reducing the “interpolant of an interpolant” error path in the weak form). It also upgrades both semi-Lagrangian schemes’ RK2 trace-back to use a mid-time velocity extrapolation for second-order accuracy in unsteady flows.

Changes:

• Add uw.discretisation.IntegrationPointVariable backed by a new PETSc “delta space” quadrature-point element, including evaluation semantics and guards.
• Add uw.systems.ddt.IntegrationPointSemiLagrangian and introduce mid-time velocity extrapolation support in the semi-Lagrangian trace-back.
• Add tests and developer documentation covering the new element, variable semantics, guards, and semi-Lagrangian behavior.

File summaries

File	Description
tests/test_0064_quadrature_point_fe.py	Validates the quadrature-point (“delta”) element tabulation and dof layout.
tests/test_0065_integration_point_variable.py	Tests IntegrationPointVariable layout, assembler read semantics, evaluate semantics, and guard rails.
tests/test_0066_integration_point_slcn.py	Tests integration-point SLCN departure-point sampling accuracy, oversampling refusal, and mid-time velocity behavior.
src/underworld3/utilities/_jitextension.py	Refuses derivatives of integration-point variable symbols at JIT/codegen time.
src/underworld3/systems/ddt.py	Adds mid-time velocity snapshot logic and new IntegrationPointSemiLagrangian implementation.
src/underworld3/systems/init.py	Exposes the new DDt class through the systems API.
src/underworld3/function/_function.pyx	Implements IntegrationPointVariable evaluation as “nearest integration point in owning cell” during interpolation.
src/underworld3/discretisation/enhanced_variables.py	Adds public IntegrationPointVariable wrapper and base-class delegation hook.
src/underworld3/discretisation/discretisation_mesh.py	Introduces mesh.integration_rule and solver-side verification against the mesh’s integration rule.
src/underworld3/discretisation/discretisation_mesh_variables.py	Adds _BaseIntegrationPointVariable implementing quadrature-point layout and nearest-point evaluation.
src/underworld3/discretisation/init.py	Exports IntegrationPointVariable from the discretisation package.
src/underworld3/cython/uw_delta_space.h	Implements and registers the uwdelta PetscSpace plugin (header-only).
src/underworld3/cython/petsc_quadrature_fe.pyx	Creates the delta FE and provides utilities (tabulation + cell quadrature point coordinates).
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Verifies solver quadrature matches mesh rule when integration-point variables are present.
setup.py	Adds the new Cython extension build for petsc_quadrature_fe.
docs/developer/subsystems/integration-point-variables.md	Adds subsystem documentation describing semantics, guards, and the integration-point SLCN.
docs/developer/index.md	Links the new subsystem documentation into the developer docs toctree.

Review details

• Files reviewed: 17/17 changed files
• Comments generated: 2
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[lmoresi]

Fixed: bcs=None with a fresh list per instance (the pre-existing SemiLagrangian has the same bcs=[] default; left as is here, outside this PR's scope).

[lmoresi]

Fixed: the message now scans PETSc's default rules for this cell type and reports the smallest qdegree with more points than the history's local dofs and the one reaching 2x oversampling (P2 on triangles: 3 and 3; P3: 3 and 5).

[lmoresi]

CI caught a units-frame bug in the new velocity cache (test_1056: units-active SLCN 63 % off the non-dimensional run). The cache variable was created without units and filled with dimensional values from evaluate, while .data is the non-dimensional store. Fixed in the last commit: the cache carries V_fn's units and every evaluated value is reduced with _to_nondim_ndarray before storing, the same idiom the nodal history uses for psi_star; the integration-point history's snapshots and slots get the same treatment, and test_1056 now runs both schemes.

Underworld development team with AI support from Claude Code

[lmoresi]

Merged development (#688, #687) into the branch. One conflict in _jitextension.py (ccode_patch_fns: kept both the integration-point derivative guard and #688's component_offsets). uw.systems.AdvDiffusion is now the composed Eulerian solver, so the tests and docs name AdvDiffusionSLCN explicitly.

We tried the integration-point history as the composed solver's transport manager. The base-class manager interface accepts it, but the composed solver's theta rule applies the diffusive flux to the history level, which differentiates the slot, and the JIT guard refused (a delta field has no gradient; silently assembling zero would have been the failure mode). So the integration-point history remains an AdvDiffusionSLCN plugin, where the diffusive history is a separate nodal DFDt. Documented.

Underworld development team with AI support from Claude Code

[lmoresi]

Reachability in the merged default solver, measured. The composed AdvDiffusion takes the integration-point history as its transport manager at order=2 (BDF2) and at order=1, theta=1, where no spatial term sits on the old level: rotating Gaussian, 16 steps, L2 1.75e-3 / peak 0.994 against 1.70e-3 / 0.991 for AdvDiffusionSLCN with the same history, and 2.6e-2 / 0.974 for the default SUPG. theta = 1 needed one change, kept inside the integration-point manager: its old-level Adams-Moulton weight is identically zero but was a runtime constant, so 0 * grad(psi*) was still differentiated at code generation and the guard fired on a dead term; the manager now returns that weight as a literal. The Crank-Nicolson flux (theta = 0.5) differentiates the old level, which a delta field cannot supply, and stays with AdvDiffusionSLCN; the guard's refusal is a test case. The base-class weights are unchanged (a first attempt there broke test_1055's weighted-sum probe).

Underworld development team with AI support from Claude Code

[Copilot]

🔵 Needs a closer look

It introduces new PETSc/Cython discretisation infrastructure and modifies core evaluation/assembly paths, which warrants careful human review beyond automated checks.

Review details

• Files reviewed: 18/18 changed files
• Comments generated: 3
• Review effort level: Lite

[lmoresi]

Measurements on this branch, prompted by two questions: does the growth get
worse at higher Courant, and is refusal the right response when the Courant
number is not something a user chooses?

The answers are no and no. The risk is at LOW Courant, and the required
stabilisation has a closed form.

The growth is a low-Courant phenomenon

Rigid rotation, P2, cellSize=0.08, qdegree=3 (12 pts/cell), 120 steps,
diffusivity 1e-8:

kind	Courant	peak/peak0	energy/E0
ip	0.25	1.0113	1.0029
ip	0.5	1.0063	1.0011
ip	1.0	1.0046	0.9977
ip	2.0	1.0114	0.9810
nodal	0.25	0.9528	0.9722
nodal	2.0	0.9542	0.9379

Energy is gained at Courant 0.25 and lost at Courant 2. Over 600 steps at
Courant 0.25 it reaches 1.0102 — a slow creep, not a blow-up, with an
adequately sampled rule. The design note's Courant 0.25 measurement is
therefore the worst case, not a representative one, and the note does not say
so.

This matters because the Courant number is chosen in the mean: in a real run
many cells sit at C much less than 1 whatever timestep is picked, so the
regime cannot be avoided by configuration.

It is not a sub-grid null space

The cell-local sampling matrix B (basis values at the rule's points, P2 on a
triangle) is well conditioned:

rule	points	dofs	cond(B)	smin/smax
3-point	3	6	—	rank-deficient
6-point	6	6	3.27	0.306
12-point	12	6	3.24	0.309

So the extra points are not helping by improving cell-local conditioning, and
there is no near-null direction within a cell. Only a rule with fewer points
than dofs has a genuine null space, which is the case worth refusing. This is
consistent with the design note's own explanation — the aliasing is between
the sampled norm and L2 for a mode at the grid scale, which is an INTER-cell
effect.

The required damping has a closed form

A grid-scale mode decays per step like 1 - kappa pi^2 dt / h^2. To beat a
per-step growth eps:

kappa_stab  =  eps * h^2 / (pi^2 * dt)          equivalently   Pe_max = pi^2 C / eps

The velocity cancels. Tested by sweeping cell Péclet at two Courant numbers,
200 steps, eps ~ 1e-3 for the 12-point rule:

Courant	Pe*	cell Pe	Pe/Pe*	energy/E0
0.25	2467	500	0.20	0.9814
0.25		2000	0.81	0.9986
0.25		8000	3.24	1.0030
0.25		100000	40.5	1.0043
0.10	987	500	0.51	1.0013
0.10		2000	2.03	1.0082
0.10		8000	8.11	1.0099
0.10		100000	101	1.0105

Pe* falls 2467 -> 987 as C falls 0.25 -> 0.1 and the measured crossover moves
with it; expressed as Pe/Pe* the two Courant numbers collapse onto one curve
with the transition near 1.

Earlier, at Courant 0.25 over 300 steps, cell Pe 1000 was already enough to
flip the scheme from gaining to losing energy (0.9887), while it still stayed
closer to conservative than nodal at the same Péclet (0.9238).

Suggestion

Replace the hard refusal at 2x with a computed stabilisation:

• add kappa_stab = eps h^2 / (pi^2 dt) per cell, eps a measured property of
  the rule and element (not a user knob);
• report the maximum kappa_stab / kappa_physical after setup, so a value that
  would matter is visible rather than silent;
• keep the refusal only for a rank-deficient rule, where no diffusivity helps.

Three properties make this behave: it is local, so slow cells are stabilised
and C~1 cells are not; it vanishes as dt grows; and kappa_stab/kappa_phys is
the inverse Fourier number, so wherever a cell is diffusive at all it is
invisible.

One honest difference from SUPG, and the reason for reporting rather than
applying silently: SUPG's stabilisation is consistent — proportional to the
residual, vanishing for the exact solution. This one is not. It adds to the
physical diffusivity unconditionally.

Not established

eps ~ 1e-3 is fitted from these runs rather than derived, and at C = 0.1 the
crossover sits nearer Pe/Pe* ~ 0.5 than 1, so eps grows somewhat as C falls —
the law's shape holds, the constant needs pinning per rule and element. One
geometry only (rigid rotation), scalar P2, order=1, theta=1.

Scripts are in this session's scratch; happy to fold them into
tests/test_0066_integration_point_slcn.py as a growth-factor check if that
is useful.

Underworld development team with AI support from Claude Code

[lmoresi]

Reviewing the growth-factor analysis above rather than the PR as a whole. It is
convincing on the mechanism and on the shape of the law; three things I think
are still open, and the third bears on whether the fix is worth taking.

1. Conditioning rules out one explanation without supplying another

The table shows cond(B) is 3.27 at 6 points and 3.24 at 12 — indistinguishable.
So cell-local conditioning cannot be why 12 points behaves and 6 does not
(1.0003 against 1.03–1.14). That is a good negative result.

But the conclusion drawn from it, that this is therefore an inter-cell effect,
does not follow on its own: if the mechanism were purely inter-cell, adding
points within a cell should not help, and measurably it does — monotonically,
6 -> 9 -> 12.

The explanation that fits both facts is sampling density rather than
conditioning: more points make the sampled norm a better approximation to L2, so
there is less for the grid-scale mode to alias into. That is cell-local in
origin even though the mode it aliases is not. Worth stating explicitly, because
it changes what the guard should key on — density relative to the mode being
resolved, not points-versus-dofs, and not purely a step-size question either.

2. eps drifting with C partly cancels the law it appears in

Pe_max = pi^2 C / eps has C in the numerator, and eps is reported to grow
as C falls (crossover at Pe/Pe* ~ 1 at C = 0.25, ~ 0.5 at C = 0.1, i.e. eps
roughly doubled). So the two C-dependences work against each other and the net
scaling of Pe_max with C is weaker than the formula suggests.

That is not a reason to reject the law — its shape clearly holds — but it does
mean the constant cannot be fitted once and reused. Two Courant numbers is
enough to see the drift and not enough to characterise it; a third point (C =
0.05) would say whether eps is heading somewhere finite.

3. Stabilising spends the conservation that motivates the scheme

This is the one I would want settled before adopting the computed diffusivity.

The integration-point history's case rests on conservation: the PR reports
energy over 63 steps at Courant 0.5 as -0.03 % against -1.4 % nodal, and
the maintainer's position is that the extra cost buys exactly that.

The proposed stabilisation is, as the comment says itself, not consistent — it
does not vanish for the exact solution, it adds to the physical diffusivity
unconditionally. The measurement above shows what that costs: at cell Pe 1000,
C = 0.25, 300 steps, energy goes to 0.9887. Still better than nodal's 0.9238
at the same Peclet, but that is a 1.1 % loss against the 0.03 % that made the
scheme worth its price.

So the trade is not "fix a slow creep for free" — it is "give back most of the
conservation advantage to remove a 1 % gain over 600 steps". Those may still be
the right terms, but they should be chosen deliberately rather than absorbed
into a default. The reporting suggestion (kappa_stab / kappa_phys after setup)
is right and necessary; I would go further and not apply it by default at all
until there is a case where the creep actually breaks a result.

On the test offer

Yes, and it should assert the scaling, not a number. A single growth factor
at one Courant is precisely what made the design note's worst case read as
representative — the comment's own first finding. A check across at least C =
0.1, 0.25, 1 that the trend is monotone and bounded would have caught that, and
would catch a future change to the rule or the element that shifts eps.

[lmoresi]

Follow-up, after Louis exercised the scheme: the stabilising term should ship
switched off.

The recommendation is to keep it as a documented option, inactive for all but
expert users — not as a default, and not as the hard refusal it replaces.

Three reasons it should not be on by default:

• it is not consistent (unlike SUPG it does not vanish with the residual), so
  it modifies the physics unconditionally;
• the effect it corrects is a slow creep on an adequately sampled rule — 1.0102
  over 600 steps at Courant 0.25 — not a failure mode most runs will meet;
• any real diffusion swamps it. kappa_stab/kappa_phys is the inverse Fourier
  number, and cell Pe 1000 was already enough to flip Courant 0.25 from gaining
  to losing energy.

Suggested shape, for whoever picks this up:

ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(
    mesh, T, V, degree=2, order=1,
    stabilise=False,      # default; True adds kappa_stab per cell
)

• stabilise=False by default. Off, and silent.
• stabilise=True adds kappa_stab = eps h**2 / (pi**2 dt) per cell and
  reports the maximum kappa_stab / kappa_physical, so a value that would
  matter is visible rather than silent.
• keep the existing refusal ONLY for a rank-deficient rule (fewer points than
  local dofs), where no diffusivity helps. The current 2x warn/refuse can then
  become a warning that names the option.

The derivation, the measured constants and the caveats are in the comment
above; the short version is that the velocity cancels, so the term is
eps h^2 / (pi^2 dt), local per cell, and it vanishes as dt grows — the C~1
cells pay nothing for the C much-less-than-1 cells being stabilised.

Worth documenting even though it is off: the reason to write it down is that
the failure it addresses is real, and the formula is not obvious to rederive.

Underworld development team with AI support from Claude Code