Eulerian SUPG advection-diffusion with symbolic BDF and Adams-Moulton orders (#657 follow-on)

[lmoresi]

Credit

The SUPG weak form used here, writing the Petrov-Galerkin test-function perturbation as a flux so PETSc needs no modified test space, and its first working implementation on PetscDS with P2 elements are @NengLu's, from the levelset branch of issue #657, together with the LeVeque swirling-flow comparison against SLCN and the conservative level-set pipeline (reinitialisation and mass correction) that motivated it. This PR keeps that idea and its tau parameter, builds the time integration on the symbolic history machinery, and adds the measurements; it would not exist without that prototype. A follow-on PR brings the level-set solver from the same branch onto this class.

Follows issue #657 (@NengLu's SUPG and level-set branch). This PR builds the Eulerian transport solver on the symbolic time-history machinery so that the scheme's order is a construction argument, adds the analytic oracle and the tests that measure it, and fixes two defects found on the way, one of them general.

What is here

uw.systems.AdvDiffusionSUPG(mesh, T, V_fn, order=2, integrator="bdf") (systems/advection_diffusion_eulerian.py). Fully implicit Eulerian advection-diffusion with SUPG stabilisation. Every past time level is a mesh variable held by the Eulerian history manager, so BDF 1-3 and Adams-Moulton 1-3 (theta at order 1, so am, 1, 0.5 is Crank-Nicolson) are assembled from the same stored history, including the diffusive flux of past states. The SUPG term is the test-function perturbation written as a flux, as in the prototype. The timestep, multistep coefficients and the tau weights are runtime constants of the compiled kernels: a change of timestep costs nothing (the prototype recompiled on every change, 1.2 s against 0.03 s for a step). Diffusivity is set on the constitutive model like every scalar solver. The class reproduces @NengLu's Crank-Nicolson prototype to four digits.

uw.analytic.RotatingGaussian: a Gaussian carried round the origin by rigid rotation while diffusing, exact at every time, so a transport scheme's error is measured rather than eyeballed. Also fixes AnalyticSolution.error(norm="integral") for scalar variables, which had never been exercised.

Tests (named to match the CI batch globs): test_1055 (API, no-recompile contract, adapt-child solve), test_1100_advdiff_supg_rotating_gaussian (temporal order 1 and 2, band invariance, round trip), tests/parallel/test_1077 (np=2 equals serial to 1e-12), test_1058 (the aux-layout fix below).

Design note docs/developer/design/eulerian-supg-transport.md with the measurements, and an example docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py.

What was measured

Rotating Gaussian, P2, unstructured simplex box, one revolution, relative L2 error. The implicit scheme is stable at any cell Courant number, and cells the scalar does not need are free: a band refined to h/9 at local Courant 13 changes the error in the third digit only. Its accuracy is set by u dt against the feature width (dt^2 for the second-order schemes). SLCN's error is flat in dt but accumulates at small Courant (21% against 0.6% at Courant 0.5 on the same mesh) and costs 4 to 6 times more per step in serial; its limit is the trace-back arc, about 10 degrees per step.

The integrator study (res 32, Courant 0.25 to 8, pure advection and kappa 1e-3) decides the default: Adams-Moulton above order 1 blows up from Courant 1 (bounded stability region), BDF3 fails from Courant 4, Crank-Nicolson is three to four times more accurate than BDF2 at the same timestep but rings once the feature is under-resolved in time, and backward Euler carries 20 to 40% error at any practical timestep. Cost per step is identical across schemes. BDF2 is the robust default; the note says when to pick the others.

Defects fixed

1. Auxiliary fields packed and indexed by position (Mesh.update_lvec, _jitextension). mesh.vars holds variables weakly; a variable released by uw.reset_default_model() (the test suite does it between tests) or by the statistics helpers' deliberate del mesh.vars[...] leaves its PETSc field in the DM. Every later variable was then packed into, and read from, the wrong slots: measured as a P0 cell-size field landing in a P2 slot as garbage, NaN residuals in one run and a subtly wrong answer in the next depending on when the collector ran. Now packed by field name (orphans zeroed) and offsets read from the DM's own field list. test_1058 fails 2 of 3 checks without the fix. This predates the PR and affects any solver with auxiliary fields; it could be cherry-picked on its own.
2. Mesh-owned multigrid pickup on a solver that owns its preconditioner (custom_mg.build_transfers): installing the adapt child's PCMG hierarchy on an additive-Schwarz PC segfaulted inside PETSc. _pc_option_prefix is None now counts as the explicit choice it is.

Not in this PR (recorded in the design note)

ALE mesh velocity (phase 1 keeps the default REMAP transfer with the material velocity, which is already correct to interpolation accuracy), a streamline element length from a mesh-owned metric tensor, a corrected discontinuity-capturing term, and the trajectory-based estimate_dt for SLCN (separate small PR). The level-set wrappers and scripts stay on the levelset branch for their author to rebase onto this solver.

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

[lmoresi]

Adversarial review of this PR, by the same session that wrote it. Findings first, evidence with each.

1. The strong residual has no diffusion term, so SUPG is inconsistent for P2 fields with kappa > 0. PETSc kernels see first derivatives only; the class documents this and the study ran at cell Peclet ~40 without visible harm, but nobody has measured the P2 error with kappa large enough for the missing term to matter. A lagged projected Laplacian is the known repair. Not fixed here.

2. tau uses mesh.cell_size(), which is volume^(1/d), not an edge length, and is partition-dependent by ~1e-3 (documented TODO in _assemble_cell_size). The weights (2, 2, 4) were kept from the prototype and shown to reproduce its numbers; they were not calibrated against the streamline length. The np=2 test passes to 1e-12 on a uniform mesh, which does not exercise the seam dependence.

3. BDF3 silently falls back to variable-step BDF2 when consecutive timesteps differ by more than 5% (ddt._bdf_coefficients). A user who varies dt gets second order without being told. Documented in the docstring, not warned at runtime.

4. Adams-Moulton orders 2 and 3 are accepted by the constructor and blow up from Courant 1 on advection. We kept them with a docstring warning rather than refusing them. A runtime check is impossible without knowing the flow; a warning at construction would be cheap. Judgement call, flagged.

5. Orphaned DM fields are zeroed at packing but never reclaimed. The fix makes the layout correct; the DM still grows by one field per dropped variable for the life of the mesh. A script that creates variables in a loop and drops them will see its aux vector and its kernels grow. That is the pre-existing memory behaviour, now merely correct.

6. The default order=2 changes nothing for existing users (new class), but theta is silently ignored at AM order 2 and 3. We raise if theta != 1 there, so it is loud, not silent. Verified by test_scheme_arguments_are_validated.

7. The rotating-Gaussian convergence slopes for BDF1 are 0.80 and 0.88, not 1. The test tolerance is 0.35 below the expected order. The pre-asymptotic behaviour is real (backward Euler saturates fast on this problem); a tighter test would need smaller timesteps and a finer mesh than a level_2 test can afford.

8. The study's res-64 rows were still running at PR time. The res-32 table is complete and decides the default; res 64 is confirmation and will be appended to the design note when it lands.

What we checked and found sound: the two families agree to every digit where they coincide (BDF1 = backward Euler); Crank-Nicolson reproduces the #657 prototype to four digits at two Courant numbers; a timestep change keeps the JIT cache key and changes the answer by 2e-11 against a 2e-2 negative control; the aux-layout regression test fails 2 of 3 checks with the fix reverted; core batches test_00xx to test_02xx pass (398) with the layout fix in place.

[Copilot]

🟡 Changes recommended

The new Eulerian solver currently uses mesh.dim where the scalar solver infrastructure expects embedded-coordinate (mesh.cdim) vector sizes, which will break manifold meshes, and it drops the verbose flag when delegating to the base SNES solve.

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

Pull request overview

Adds a fully implicit Eulerian advection–diffusion SUPG transport solver built on the existing symbolic time-history (ddt.Eulerian) machinery, plus an analytic “rotating Gaussian” oracle and a suite of tests/notes validating temporal order, timestep-constant (no-recompile) behavior, mesh-adaptation invariance, and parallel partition-independence. It also fixes two infrastructure defects uncovered during development: auxiliary-field packing/offseting after dropped MeshVariables, and a multigrid pickup edge case when a solver explicitly owns its PC.

Changes:

• Introduce uw.systems.AdvDiffusionSUPG (SNES_AdvectionDiffusion_SUPG) with BDF/Adams–Moulton order selection and runtime-constant timestep/coefficients.
• Add uw.analytic.RotatingGaussian and tests that quantitatively measure temporal convergence, refinement invariance, and serial/parallel consistency.
• Fix auxiliary-field layout robustness (Mesh.update_lvec, JIT aux offsets) and guard mesh-owned MG pickup for solvers with _pc_option_prefix is None.

File summaries

File	Description
src/underworld3/systems/advection_diffusion_eulerian.py	New Eulerian implicit SUPG advection–diffusion solver implementation and API.
src/underworld3/systems/__init__.py	Exports AdvDiffusionSUPG.
src/underworld3/systems/ddt.py	Exposes BDF/AM coefficient symbols for solvers assembling their own weighted history sums.
src/underworld3/systems/solvers.py	Refactors shared adv/diff timestep estimation into reusable helpers.
src/underworld3/utilities/custom_mg.py	Prevents MG pickup when solver explicitly owns its PC (_pc_option_prefix is None).
src/underworld3/utilities/_jitextension.py	Fixes auxiliary-field component offsets to account for orphan DM fields.
src/underworld3/discretisation/discretisation_mesh.py	Fixes Mesh.update_lvec() to pack aux data by DM field name and zero orphans.
src/underworld3/analytic/transport.py	Adds RotatingGaussian analytic transport oracle.
src/underworld3/analytic/_base.py	Fixes AnalyticSolution.error(norm="integral") for scalar mesh variables (1x1 Matrix symbol).
src/underworld3/analytic/__init__.py	Exports RotatingGaussian.
tests/test_1055_advdiff_supg_api.py	API/no-recompile/argument validation and adapt-child PC ownership test.
tests/test_1058_dropped_meshvariable_aux_layout.py	Regression tests for dropped-variable orphan DM field layout bug.
tests/test_1100_advdiff_supg_rotating_gaussian.py	Accuracy/temporal order, refinement invariance, and round-trip tests against analytic oracle.
tests/parallel/test_1077_advdiff_supg_parallel.py	Parallel partition-independence test against a serial reference error.
docs/developer/design/eulerian-supg-transport.md	Design note capturing decisions and measurements.
docs/developer/index.md	Adds new design note to the developer docs toctree.
docs/advanced/semi-lagrangian-time-integration.md	Cross-references the new Eulerian alternative.
docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py	Example script demonstrating solver usage and measured error against the oracle.

Review details

Suppressed comments (3)

src/underworld3/systems/advection_diffusion_eulerian.py:219

• The velocity coercion should be done in embedded coordinate dimension (mesh.cdim) to match SNES_Scalar's expectations for flux/gradients; using mesh.dim will break manifold meshes and can mismatch F1's required size.

        self.f = sympy.Matrix.zeros(1, 1)
        self._integrator = integrator
        self._time_order = order
        self._theta = theta
        self._V_fn = _as_row_vector(V_fn, mesh.dim)

src/underworld3/systems/advection_diffusion_eulerian.py:369

• _advection() uses mesh.dim to form u·∇phi, but SNES_Scalar kernels (and UW's variable gradients) are defined in embedded coordinates (mesh.cdim). On manifold meshes this will either be wrong or fail when V_fn has cdim components.

    def _advection(self):
        dim = self.mesh.dim
        u = self._V_fn
        total = sympy.Integer(0)
        for w, phi in zip(self._spatial_weights(), self._states()):

src/underworld3/systems/advection_diffusion_eulerian.py:380

• _diffusive_flux() builds a (1, mesh.dim) vector, but SNES_Scalar reshapes fluxes to mesh.cdim. This will raise at runtime on manifold meshes (dim!=cdim). Build the flux in embedded dimension (cdim) instead.

    def _diffusive_flux(self):
        r"""``(1, dim)`` flux :math:`\sum_k w_k\,\nabla\phi^{(k)}\cdot\kappa` from the constitutive tensor."""
        dim = self.mesh.dim
        c = self.constitutive_model.c
        total = sympy.zeros(1, dim)

• Files reviewed: 18/18 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]

Res-64 rows of the integrator study are in (design note updated in the last commit). BDF2 and Crank-Nicolson track their res-32 errors at the same u dt. BDF3 is not safe for pure advection at any Courant number: its stability region misses the imaginary axis near the origin, so the finer mesh's low-frequency modes grow slowly (31x the exact field after 590 steps at Courant 1); it behaves with kappa = 1e-3. Review item 8 is closed; the BDF2 default stands.

[lmoresi]

The auxiliary-field layout fix (commit e84dea9 and its test) is extracted as #676 so it can land on development ahead of this PR; this branch will rebase onto it once merged.

[lmoresi]

Interface change after discussion: AdvDiffusionSUPG now takes the semi-Lagrangian solver's signature and semantics so it is a drop-in replacement (change the class name, nothing else). order=1 with theta=0.5 is Crank-Nicolson and the default, as for SLCN; order=2 takes theta=1 (BDF2) unless 0.5 is asked for explicitly, which is refused for the reason the SLCN docs give; the trace-back-only arguments are accepted and ignored with a warning; delta_t is settable and solve() reuses it. New user page docs/advanced/eulerian-advection-diffusion.md with the swap table and when to use which solver. Review items 4 and 6 stand as written; the default moved from BDF2 to Crank-Nicolson by the drop-in contract, and the design note records why that is also the better default on the measurements.

[lmoresi]

Removed the integrator argument: order and theta already reach every scheme that is safe (Crank-Nicolson, backward Euler, BDF2, BDF3 with diffusion), and the only settings the argument added were the Adams-Moulton orders 2 and 3 that the study shows blowing up on advection from Courant 1. The family now follows the order; the higher Adams-Moulton assembly remains reachable only by switching it on the instance, for the record. Review item 4 is closed by removal.

[lmoresi]

Timestep: estimate_dt() on the Eulerian solver is now accuracy-based. It returns the step at which the field changes by fraction (0.02) of its range, from the advective rate before the first solve and from the realised rate of the last step after it, so it does not depend on the mesh (the band test checks that the resolution estimate collapses 3x on the refined child while this one moves under 25%). The cell-crossing time SLCN reports is neither a stability nor an accuracy limit here; it stays available as basis='resolution'. At the default fraction the Crank-Nicolson round trip on the rotating Gaussian is under one per cent (test_1100). Credit to @NengLu is now in the PR description, the module docstring and the design note.

[lmoresi]

Linear solver settled (bd0b1a5): matched tolerances, Schwarz by default, preconditioner="fmg" as a real switch.

We set out to measure geometric multigrid against the solver's additive-Schwarz ILU on the level-set advection step (256² and 512² quad box with a refinement hierarchy, P2, Crank-Nicolson; script and logs in ~/+Simulations/supg_vs_slcn_657/parallel/). The first result read as a 1.65x win for multigrid (0.94 s against 1.56 s per step in serial, one F-cycle against three GMRES iterations). -log_view said otherwise: the Schwarz run was doing two Newton iterations per step on a linear operator, because the Krylov default tolerance (1e-5) does not reach the SNES tolerance (1e-8), and that second Jacobian assembly cost more than every linear solve of the step. The linear solve is under a tenth of a step either way; assembly is the rest.

With the Krylov tolerance at 1e-9 (now the default) both do one Newton step, and the comparison at 256² is:

Courant	ASM-ILU its (np 1 / 8)	s/step	FMG cycles (np 1 / 8)	s/step
1/2	5 / 5	0.913 / 0.121	1 / 1	0.943 / 0.145
2	8.9 / 8.6	0.925 / 0.141	3.4 / 3.6	1.079 / 0.178
8	16.6 / 16.5	0.971 / 0.146	12.8 / 12.8	1.657 / 0.299
32	37 / 37.8	1.128 / 0.172	23.8 / 24.1	2.347 / 0.457

Schwarz is cheaper at every Courant number and its count is the same on one and eight ranks. The multigrid's cycle count grows with the Courant number nearly as fast (the Galerkin coarse operators inherit the fine-grid tau; four levels equal three at Courant 1/2, so they are not under-stabilised there), and a cycle costs about three Schwarz iterations. At 512² with four levels the same holds: 3.51 / 0.48 s (Schwarz) against 3.62 / 0.54 s (multigrid). GAMG as a control was slowest (5 iterations, 2.07 s serial). Smoother variants made no difference. So the default is GMRES + ASM-ILU on every mesh, and solver.preconditioner = "fmg" hands the block to the managed multigrid route for the rank count where a one-level method has no coarse space, which we have not measured beyond eight ranks.

Two defects found on the way:

• The solver's solve() pre-ran the three setup stages, which marks the solver set up, so the base _build returned before resolving the preconditioner choice and preconditioner = "fmg" was silently inert. It now builds through _build. The semi-Lagrangian solvers share the pattern and the defect: preconditioner="fmg" is inert on the semi-Lagrangian solvers: their solve() pre-runs setup and skips the base _build resolution #683.
• The parallel timings in the design note were taken with the unmatched tolerance and carry the second Newton step; the note says so now.

Tests: test_1055 gained the switch test (Schwarz and multigrid agree to the solve tolerance, and switching back rebuilds); the adapt-child test is unchanged in what it asserts. 18 serial and the np=2 parallel test pass.