Conservative level-set transport on the shared solver interface (NengLu's #657 work)

[lmoresi]

Credit

This is @NengLu's conservative level-set work from issue #657 (branch levelset), kept in its own commits at the base of this branch: the tanh-profile initialisation from a signed distance or a polygon, interface thickness from the local cell size, the Parameswaran and Mandal (2023) reinitialisation integrated with SSP-RK3, the Zhang, Zou and Greaves (2010) global mass correction, material properties blended across the interface, and the LeVeque swirling-flow comparison of SUPG against SLCN. The commits on top bring that work onto the shared solver interface of #673 and fix what stopped it running.

What changed on top

• One module, both transport solvers. The two copies of the pipeline (level_set_SUPG.py and level_set_SLCN.py) become uw.systems.LevelSetSolver(psi, velocity=..., epsilon=..., advection="supg"|"slcn") in systems/level_set.py, with initialise_psi, interface_thickness and material_property_field alongside. The Eulerian solver is the default; order and theta pass straight through to either transport solver (Crank-Nicolson by default, which preserves the profile's amplitude between reinitialisations).
• The prototype SUPG solver is retired in favour of AdvDiffusionSUPG from Eulerian SUPG advection-diffusion with symbolic BDF and Adams-Moulton orders (#657 follow-on) #673, which reproduces it to four digits, and the structured-grid ENO reinitialisation that was already switched off in the SUPG variant is dropped; the projection of the gradient magnitude is the general path.
• Fixed: the wrapper called the transport solver's solve() positionally after it became keyword-only, so the LeVeque script could not run as pushed; shapely is now an optional import with a clear message (signed_distance= needs none of it); deprecated data access (mesh.data, .data[:, 0], mesh.access) replaced.
• Examples and tests. The LeVeque comparison moves to docs/examples/convection/advanced/Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py in the repository's script conventions (uw.Params, no argparse). The two 1-D block variants duplicated the existing example and one depended on the retired discontinuity-capturing term, so they are not carried over. tests/test_1100_levelset_rotation.py (level_2) checks, for both transport choices, that a circle carried once round the box holds its volume to 1e-6 at every step, keeps its profile sharp and comes back; plus the polygon initialisation and the property blend.
• Docs: docs/advanced/level-set-transport.md.

Measured

LeVeque flow, 64 by 64, period 2, Courant 0.5: round-trip shape error 0.028 (SUPG) against 0.051 (SLCN) at half the wall time; the mass correction pins both to the same volume.

Base branch

Opened against feature/eulerian-supg-transport (#673) so the diff shows only the level-set work; to be retargeted at development once #673 merges.

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 session that wrote it. Findings with evidence.

1. The reinitialisation projects |grad psi| three times per pseudo-time step (one L2 projection per SSP-RK3 stage), and the projection solver is a full SNES solve. At 64 by 64 that is most of the wall time of a level-set step (about 2 s against 0.3 s for the advection). Inherited from the prototype; a nodal Clement gradient would be cheaper and is what the transport solver's estimate_dt uses. Not changed here.

2. The mass correction is a uniform shift of the whole field, so it repairs volume by moving the contour everywhere by the same amount, not where the error was made. Standard for the method (Zhang, Zou and Greaves), but it means the "volume drift" a script reports after correction is always zero and says nothing about the transport. The example prints the drift for that reason only as a check that the corrector ran.

3. The wall correction (adv_solver_bc) is box-only and copies the neighbouring interior row or column onto the wall nodes. It assumes a structured node layout along the wall and is silent (a warning) when it cannot match a node. It is off by default and documented as a box-mesh convenience; it should not be used on unstructured meshes, and nothing stops a user from trying.

4. interface_thickness loops over cells in Python calling computeCellGeometryFVM per cell. Fine at these sizes, slow at 10^6 cells. The mesh's own per-cell size field (mesh.cell_size()) is the obvious replacement but is partition-dependent by ~1e-3, so the loop stays for now.

5. The default reinitialisation frequency is an empirical formula (4.9e-3 * domain size / min epsilon - 0.25, rounded) inherited from the prototype, with no reference. It gives every step on coarse meshes and less often as they refine; the rotation test and the LeVeque example override it.

6. Parallel. Every collective in the wrapper (volume integral, projections, the allreduces in the wall correction and the epsilon minimum) is unconditional, and the rotation test runs serially only. No np=2 run was made for this PR; the transport solvers underneath have theirs.

7. Dropped, deliberately: the structured-grid ENO reinitialisation (the SUPG variant already had it switched off), the two 1-D block scripts (one used the retired discontinuity-capturing term, both duplicated the existing example), and the SLCN variant's old_frame_traceback=True (meaningless on a static mesh).

Checked and sound: both transport choices pass the rotation test (volume to 1e-6 per step, profile sharp, circle back); the module imports without shapely; the style gate is clean; the branch keeps the original commits and authorship.

[lmoresi]

Mass correction: the bisection (15 to 30 integrals per step, most of a level-set step's cost) is replaced by a bracketed secant iteration on the same integral, started from the transition-band area as the slope. Same volume to ten digits on a perturbed circle at 64x64 with 5 integrals against 32 (0.20 s against 0.59 s). A lumped-mass version (one numpy sum per trial) would remove the integrals altogether but needs the PETSc global-to-local ordering plumbing; left for later.

[lmoresi]

Added far_field=: the value of psi imposed on every mesh boundary. Without it a CG transport lets mass in through inflow boundaries (4% volume drift in 20 steps of a rotating circle that crosses the box walls; 8e-5 with it). With the walls handled, SUPG conserves the enclosed volume to solver tolerance and the reinitialisation changes it only at second order, so the global shift corrector is really there for SLCN's interpolation loss and the clipping. Rotation test now uses it; 4 tests pass with the secant corrector in place.

[lmoresi]

Three more things on this branch: (1) conserve_mass='auto' (on for SLCN, off for SUPG) with volume_drift reported; the field is clipped to [0, 1] after each stage regardless. (2) Band thickness matters more than anything else for the CG transport: interface_thickness(scale=0.35), the g-adopt default, is a band under one cell (eps = h/8) and SUPG rings at it (0.84% clipped volume per revolution with no corrector); scale 2.0 (eps = 0.7 h) removes the ringing and leaves 0.18% per revolution, which is the reinitialisation's curvature error. Documented; the test uses scale 2. (3) Parallel timings (128x128 and 256x256, 1 to 8 ranks) are in the design note: SUPG advection 7x cheaper in serial, about 5x at 8 ranks, partition-independent to ten digits; the SLCN answer moves with the partition, filed as #682.

[lmoresi]

Heads-up: #673 is closed as superseded by #688 (which targets development now). This branch needs the new names before it can merge: the composing scalar solver is uw.systems.AdvDiffusion (AdvDiffusionSUPG no longer exists) and its transport is the uw.systems.ddt.EulerianSUPG manager. It will be rebased onto development after #688 lands.

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

[lmoresi]

#688 is merged to development (1d5e8c1), so this PR now targets development; it has 10 conflicts against it and still uses names that no longer exist. What to rebase onto: the scalar solver is uw.systems.AdvDiffusion (the composing solver; AdvDiffusionSUPG is gone), its transport is the uw.systems.ddt.EulerianSUPG manager (supg_weight, tau_weights, tau_shape, peclet_weight pass through from the solver), and the residual is composed from DuDt.time_derivative() + DuDt.advection() and DuDt.stabilisation_flux(R) rather than assembled in the solver. Design note: docs/developer/design/eulerian-supg-transport.md, section "The DDt as the transport plugin".

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

[Copilot]

Copilot review overview

🟡 Changes recommended

Import failures, incorrect solver dispatch, unsupported material APIs, MPI defects, and free-surface method regressions currently block approval.

Review effort: Balanced
Findings: 14 · 4 · 2

Open (20)

• Collective integral evaluation is called only on rank zero · New
• MPI ranks concurrently overwrite the same surface CSV · New
• Fixed epsilon variable name causes cross-solver interference · New
• SUPG branch incorrectly instantiates the SLCN solver · New
• Local band widths cause inconsistent MPI interface shifts · New
• Sharp adjoint produces NaN on the level-set interface · New
• Deprecated access context fails the repository style gate · New
• Import references nonexistent systems.level_set module · New
• Missing MaterialSwarm API breaks particle composition handling · New
• Deprecated mesh-variable access pattern fails the style gate · New
• MaterialCLS solvers double-correct remapped mesh motion · New
• SUPG path silently uses the SLCN solver · New
• Mesh velocity is subtracted twice during remapping · New
• Relative velocity is written into remap-managed solvers · New
• Left-endpoint velocity sampling reduces transport accuracy · New
• Parallel surface output omits nonlocal boundary nodes · New
• Undocumented high-risk FreeSurface changes expand PR scope · New
• Global relaxation rate replaces required pointwise rates · New
• Examples violate the required uw.Params configuration convention · New
• Examples parse sys.argv instead of using uw.Params · New

What changed in this PR

Adds conservative level-set transport and Eulerian SUPG advection-diffusion, with analytic validation, documentation, and free-surface integration.

Changes:

• Introduces SUPG transport, level-set initialization, reinitialization, and mass correction.
• Adds solver, parallel, rotation, and auxiliary-field regression tests.
• Extends free-surface material transport and supporting examples.

File	Description
src/​underworld3/​level_set.py	Implements level-set and material APIs.
src/​underworld3/​systems/​advection_diffusion_eulerian.py	Implements Eulerian SUPG transport.
src/​underworld3/​systems/​free_surface.py	Extends free-surface material transport.
src/​underworld3/​systems/​solvers.py	Shares timestep estimation logic.
src/​underworld3/​systems/​ddt.py	Exposes multistep coefficients.
src/​underworld3/​systems/​__init__.py	Exports new solver APIs.
src/​underworld3/​discretisation/​discretisation_mesh.py	Corrects auxiliary-field packing.
src/​underworld3/​utilities/​_jitextension.py	Corrects auxiliary-field offsets.
src/​underworld3/​utilities/​custom_mg.py	Guards multigrid pickup.
src/​underworld3/​analytic/​transport.py	Adds rotating-Gaussian solution.
src/​underworld3/​analytic/​_base.py	Fixes scalar integral errors.
src/​underworld3/​analytic/​__init__.py	Exports rotating Gaussian.
src/​underworld3/​__init__.py	Exposes level-set module.
tests/​test_1055_advdiff_supg_api.py	Covers SUPG API behavior.
tests/​test_1058_dropped_meshvariable_aux_layout.py	Tests orphaned auxiliary fields.
tests/​test_1100_advdiff_supg_rotating_gaussian.py	Measures SUPG accuracy.
tests/​test_1100_levelset_rotation.py	Tests level-set transport.
tests/​parallel/​test_1077_advdiff_supg_parallel.py	Checks partition independence.
docs/​examples/​convection/​advanced/​Ex_AdvectionDiffusionSUPG_RotationTest.py	Demonstrates SUPG rotation.
docs/​examples/​convection/​advanced/​Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py	Compares level-set solvers.
docs/​examples/​free_surface/​advanced/​Ex_Crameri2012Case1_uw3_FreeSurf.py	Adds Case 1 example.
docs/​examples/​free_surface/​advanced/​Ex_Crameri2012Case2_uw3_FreeSurf.py	Adds Case 2 example.
docs/​advanced/​eulerian-advection-diffusion.md	Documents SUPG usage.
docs/​advanced/​level-set-transport.md	Documents level-set transport.
docs/​advanced/​semi-lagrangian-time-integration.md	Links the Eulerian alternative.
docs/​advanced/​index.md	Registers advanced guides.
docs/​developer/​design/​eulerian-supg-transport.md	Records design and measurements.
docs/​developer/​index.md	Registers the design note.

---

💡 Configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[NengLu]

Details about the comment 8a85ee2 have been provided here to enhance issue #657.

Multi-material free surfaces: MaterialCLS, swarm support in FreeSurface, and a correct relaxation rate

Summary

Three related changes to underworld3/level_set.py and
underworld3/systems/free_surface.py:

1. MaterialCLS — a third MaterialDistribution, putting materials on conservative level sets with the same interface MaterialSwarm and MaterialRegions already present. stokes.materials = materials now works identically for all three, with no special-casing.
2. FreeSurface accepts all three material representations — a single MeshVariable, a MaterialCLS, or a MaterialSwarm — dispatched once at construction, each transported the way its own type requires.
3. The surface relaxation rate is now a single global least-squares value instead of a per-node ratio. The per-node form manufactured node-scale artefacts out of two perfectly smooth fields; see Bug fixes below.

Benchmarked against Crameri et al. (2012) Cases 1 and 2 — results at the end.

---

1. MaterialCLS — materials on conservative level sets

MaterialCLS subclasses the existing MaterialDistribution, so .add(), materials["name"] = region, materials.<property> and stokes.materials = materials are inherited unmodified:

materials = uw.level_set.MaterialCLS(mesh, velocity=v_mat)
materials.add("mantle", shear_viscosity_0=1.0,   density=3300)
materials.add("slab",   shear_viscosity_0=1.0e3, density=3400)
materials["slab"] = mesh.X[1] > 0.53
stokes.materials = materials
stokes.bodyforce = -materials.density * mesh.CoordinateSystem.unit_e_1
materials.solve(dt)                       # analogue of MaterialSwarm.advection()

What is genuinely different from MaterialSwarm

	behaviour	rationale
Regions	need a signed distance, not a boolean mask	the tanh profile is built from how far a point is from the interface, not which side it is on
Affine conditions	mesh.X[1] > 0.53 auto-converts to its exact distance (ℓ−r)/‖∇(ℓ−r)‖	a hyperplane's distance field is itself linear, so this is exact — same syntax as MaterialSwarm for the common case
Non-affine conditions	refused, not approximated	no closed-form distance to an arbitrary condition's boundary; an inside=+1/outside=−1 substitute would leave the tanh profile no transition zone at all
Unpainted materials	exactly one allowed — the first declared, which becomes the implicit background max(1−Σothers, 0)	mirrors the index-0-is-default convention of MaterialSwarm/MaterialRegions, translated to a continuous field. A second is ambiguous for a continuous partition and raises
Mixing rules	adds "geometric" alongside arithmetic/harmonic	the log-mean is the standard choice for a property spanning orders of magnitude (viscosity), where arithmetic under-weights a thin weak inclusion

New public API

MaterialCLS.solve, .estimate_dt, .sync_velocity_from, .level_set,.interface_volumes, .total_volume_error, .mixing, .weight, .blend,.built.

estimate_dt(basis=...) forwards basis only to materials actually using advection="supg". Per-material advection schemes are allowed (dict-keyedkwargs), and "slcn"'s estimate_dt has no basis concept and would raiseTypeError on a uniform forward.

---

2. FreeSurface takes any of the three material representations

Previously composition was a single MeshVariable, always transported by one SLCN solver that FreeSurface built and owned. Now the type is resolved once into _comp_kind and every per-step path reads that discriminator:

_comp_kind	transport	advection velocity
"field"	one solver FreeSurface builds ("slcn" or new "supg")	consistent.u − u_mesh for SUPG
"cls"	the MaterialCLS's own per-material LevelSetSolvers, via .solve(dt)	consistent.u − u_mesh
"swarm"	the swarm's own advection(V_fn, dt, ...)	consistent.u in full

The asymmetry in that last column is deliberate

The ALE correction u − u_mesh exists because an Eulerian field lives at nodes that move underneath it. A particle holds an absolute coordinate and has no such frame — subtracting the mesh velocity would drag every particle backwards by exactly the mesh motion, every step. _track_mesh_velocity is therefore False for a swarm, and those two MeshVariables plus the projection solver are never built at all.

It is still consistent.u, not free.u: the consistent solve is constrained so u·n̂ = Δh/Δt on the surface, which is exactly how far the mesh surface moves that step, so a particle on the lid travels with it. Using free.u would move near-surface particles at the un-relaxed rate and open a gap that compounds.

Ordering: the swarm is advected before _carry_and_deform, so the velocity is sampled on the geometry the Stokes solves actually ran on; _sync_swarm_after_deform repairs cell ownership and proxy fields afterwards.

advection parameter

advection=None (new default) selects the only scheme that can apply — "swarm" for a MaterialSwarm, "slcn" otherwise, so existing MeshVariable behaviour is unchanged. "swarm" with a mesh field raises; a mesh scheme with a MaterialSwarm warns and is ignored. "supg" is newly available for a single MeshVariable.

MaterialSwarm.advection's optional keywords are probed once with inspect at build time rather than caught at the call site — a TypeError raised inside advection would otherwise be mistaken for a bad keyword and the
particles advected a second time. New swarm_advection_order and swarm_restore_fn pass through; the latter defaults to the mesh's own return_coords_to_bounds when present.

Already-built MaterialCLS hand-over

A MaterialCLS that also supplies stokes.bodyforce/shear_viscosity_0 is necessarily already built by the time FreeSurface can see it, and its LevelSetSolvers have compiled their velocity into their weak forms. The new
MaterialCLS.sync_velocity_from() writes new data into the same fixed-identity MeshVariable every step, which works where reassigning .velocity does not. Constructing with velocity=v_mat (a MeshVariable, not v.sym) enables it; passing a bare sympy expression warns loudly that the materials will keep advecting with the velocity they were built with.

---

3. Bug fixes

3.1 Per-node relaxation rate manufactured node-scale artefacts (Need further discussion)

Was:

gamma = np.where(np.abs(displacement) > 1.0e-9, u_n / displacement, 0.0)
gamma = np.abs(gamma)

Now: one rate for the whole surface, by FE-trace-mass-weighted least squares

γ = ∮ w ḣ (h∞ − h) ds  /  ∮ w (h∞ − h)² ds

The exponential update is derived for a surface relaxing as a single mode, and with one mode the pointwise ratio is the same number at every node — so the least-squares rate recovers it exactly, and the update stays exact at any Δt. Confirmed to machine precision:

τ	recovered γ	1/τ
0.5	2	2
24.6	0.0406504065041	0.0406504065041
10000	0.0001	0.0001

Under a multi-wavelength load it is not one mode. ḣ comes out spatially broader than h∞, with different zero crossings, and the pointwise ratio then does not approximate a rate — it creates structure present in neither field:

• where ḣ = 0 but h∞ ≠ h, pointwise γ is zero, so e^{−γΔt} = 1 and that node is frozen while both neighbours relax past it — a one-node spike, sign-flipped against the local trend. This is the damaging case, and a guard on the denominator cannot see it.
• where h∞ → h the ratio has a pole, but that one is benign: the node is driven to h∞, which is where it already is.
• where the two disagree in sign — normal near a load's flanks, since neighbouring surface points are mechanically coupled — abs() silently reverses the motion.

The old 1e-9 guard was absolute and never fired against a field of scale 1e-4. Measured on Crameri Case 2: node-to-node curvature drops 10×, from 2.98e-6 to 2.97e-7 against a 1.1e-4 signal, and the frozen node at ḣ=0
goes from relaxing 0.43 % of the way to 49.61 % with its neighbours.

3.2 Second and later LevelSetSolvers on one mesh silently shared a field

phi_grad was named from a hardcoded literal:

self.phi_grad = discretisation.MeshVariable(r"|\nabla\psi|", ...)   # was
self.phi_grad = discretisation.MeshVariable(rf"|\nabla {self.phi.name}|", ...)  # now

With a fixed name, the second LevelSetSolver on the same mesh collides on this auxiliary variable. PETSc/UW3 prints Variable ... already exists - Skipping and hands back the first instance's phi_grad — so the second material's reinitialisation reads and writes the first material's gradient-magnitude field, with no exception anywhere. This is the "silently shared level sets" hazard MaterialDistribution._check_level_sets_are_new warns about, one level deeper, inside LevelSetSolver itself.

Every multi-material script with 2+ materials on one mesh hits this — both MultiMaterialLevelSet and MaterialCLS. It is a prerequisite for MaterialCLS working at all.

3.3 AdvDiffusionSUPG → AdvDiffusion

LevelSetSolver's SUPG path now uses systems.AdvDiffusion composed with the default EulerianSUPG transport manager, tracking the current solver API. Docstring updated.

---

Validation — Crameri et al. (2012)

Both benchmark cases run from a single script per case with a MATERIAL switch, so the swarm and level-set paths share every line except the materials block.

Case 1 — cosine relaxation (analytical solution, eq. 23: τ = 14.825 ka)

run	τ_fit [ka]	error	mode purity (RMS residual)
swarm, arithmetic, res 128	14.867	+0.28 %	0.00000
cls, arithmetic, res 128	17.244	+16.32 %	0.00000
cls, arithmetic, res 256	15.837	+6.83 %	0.00000
cls, geometric, res 128	14.524	−2.03 %	0.00000

All four stay exact single cosine modes to five decimal places — no node-scale structure survives in any run, which is the direct evidence that the global relaxation rate behaves as intended and 3.1 is fixed.

The CLS spread is the smeared-interface/mixing-rule sensitivity, not a free-surface defect: 18.3 percentage points from the averaging rule alone at fixed resolution, converging at order 1.26 under refinement.

Case 2 — plume beneath the lithosphere

t [Ma]	paper [m]	swarm [m]	diff
4	460	436.4	−5.1 %
8	660	633.6	−4.0 %
12	760	749.8	−1.3 %
16	810	791.3	−2.3 %
20	830	846.6	+2.0 %

Mean absolute error 2.95 % at res 128, where the plume radius is only 2.29 cells. Paper values are digitised from Figs 6(a)/7 (±20 m), so 12, 16 and 20 Ma are inside the reading uncertainty.

---

[lmoresi]

This one is not a conflict to fix — it is two PRs in a trench coat

Merging development gives 13 conflicted files and 57 hunks, 29 of them in
src/underworld3/systems/advection_diffusion_eulerian.py, and four of the
conflicts are add/add: both sides created the same file independently.

That is the signal. This branch and development developed Eulerian SUPG
transport in parallel, and development landed its version first:

this branch	development
SNES_AdvectionDiffusion_SUPG (754 lines)	SNES_AdvectionDiffusion_Composed (714 lines, has peclet_weight)
docs/developer/design/eulerian-supg-transport.md (343 lines)	same path, its own text
Ex_AdvectionDiffusionSUPG_RotationTest.py	same path
tests/test_1055_advdiff_supg_api.py, test_1100_advdiff_supg_rotating_gaussian.py, test_1077_advdiff_supg_parallel.py	same paths

So the SUPG half of this branch is superseded, and resolving those 29 hunks
by hand would mean re-deciding, line by line, which of two implementations of
one class survives — without either having been reviewed against the other.

What is actually still unique here

The level-set and free-surface work, which is what the PR title is about and
which nothing on development touches:

file	lines
src/underworld3/level_set.py	1307
src/underworld3/systems/free_surface.py	+674
Ex_Crameri2012Case1_uw3_FreeSurf.py, Case2	688
Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py	178
docs/advanced/level-set-transport.md	93
tests/test_1100_levelset_rotation.py, test_1058_dropped_meshvariable_aux_layout.py	205

Recommendation

Rebase this branch onto development keeping only the level-set and
free-surface files, and drop the SUPG half as landed elsewhere. That turns
57 hunks into a handful and leaves a PR that says what it does. It needs an
author decision, not a merge, so the merge has been left unattempted rather
than resolved into something nobody reviewed.

Note tests/test_1100_levelset_rotation.py collides on number with
development's test_1100_advdiff_supg_rotating_gaussian.py, and
test_1055_advdiff_supg_api.py with test_1055_yield_smoother.py.

Underworld development team with AI support from Claude Code

[lmoresi]

Review before the split: what the SUPG half holds that development does not

The decision is to keep the level-set and free-surface work and drop the
add/add SUPG half. Before dropping anything, every file was compared against
development to see what would be lost. Nothing is lost, and the reason
is stronger than "development's version is newer": most of this branch has
already landed on development by another route.

The SUPG solver: development's is a strict superset

	this branch	development
class	SNES_AdvectionDiffusion_SUPG	SNES_AdvectionDiffusion_Composed
exported as	systems.AdvDiffusionSUPG	systems.AdvDiffusion
public names	10	the same 10, plus peclet_weight
τ shapes	inverse_sum only	inverse_sum, brooks_hughes, doubly_asymptotic
τ lives in	the solver (_tau)	ddt.EulerianSUPG.tau()
tau_weights default	(2.0, 2.0, 4.0)	(2.0, 2.0, 4.0)
supg_weight default	1.0	1.0

The two files differ by 232 lines. The inverse_sum τ is algebraically
identical on both sides — same (c_t c_0/\Delta t)^2, (c_u|a|/h)^2,
(c_v\nu/h^2)^2 under the same square root, same BDF c_0, same weights.
The only substantive difference is where τ lives: development moved it into
the transport manager, which is the design the rest of ddt.py follows.

Already on development, so not at risk in the first place

Three of these were flagged as independent bug fixes that a naive drop would
lose. All three are already in:

change	state on development
custom_mg.py — the _pc_option_prefix is None guard (SEGV in _configure_pcmg on an adapt child)	present
discretisation_mesh.py — traverse DM fields BY NAME, not positionally (a collected variable leaves an orphan field; the cell-size field landed in a P2 slot as garbage)	present
_jitextension.py — _aux_component_offsets read from the DM	present
tests/test_1058_dropped_meshvariable_aux_layout.py	byte-identical
solvers.py — _advective_diffusive_dt + _global_max_diffusivity + _centroid_velocities_nd	present
analytic/transport.py	byte-identical
analytic/__init__.py, analytic/_base.py	identical
ddt.py — bdf_coefficient_expressions, am_coefficient_expressions	present
docs/.../semi-lagrangian-time-integration.md — "The Eulerian alternative"	present, and written against the correct name uw.systems.AdvDiffusion (this branch's copy says AdvDiffusionSUPG)

The split needs no rewire

level_set.py already calls systems.AdvDiffusion(self.mesh, self.phi, self.velocity, order=order, theta=theta) — development's name, not the one
being dropped — and that matches Composed.__init__(mesh, u_Field, V_fn, order=1, theta=None, ...) positionally. AdvDiffusionSUPG appears zero
times across level_set.py, free_surface.py, the Crameri and LeVeque
examples, level-set-transport.md and test_1100_levelset_rotation.py.

What is actually unique, and the only real merge work

	lines
level_set.py	1307
Crameri Case1 + Case2 free-surface examples	688
LeVeque SUPG-vs-SLCN comparison	178
level-set-transport.md	93
test_1100_levelset_rotation.py	92
free_surface.py — _advect_swarm, _build_mesh_velocity_tracking, _sync_swarm_after_deform, _relaxation_rate	+627/−47

free_surface.py is the only file with a genuine two-sided conflict left, and
it is lopsided: this branch +627/−47 against development's +86/−2.

Three problems to raise

1. peclet_weight defaults to 4.0 on development, and this branch has no
such term. It multiplies the SUPG weight by
|a|²h²/(|a|²h² + 4\,p^2\nu^2 + \epsilon), damping stabilisation where
diffusion dominates. For the level set it is inert: level_set.py pins
diffusivity = 0.0, and at ν = 0 the factor is |a|²h²/(|a|²h² + ε) ≈ 1
wherever the velocity is non-zero (and 0 where it vanishes, where there is no
advection to stabilise). Nothing else in the keep set sets a diffusivity. So
it does not affect this work — but it is a real behavioural difference for any
other caller, and it is development's default, chosen without this branch in
view. Worth confirming 4.0 is intended.

2. tests/test_1100_levelset_rotation.py collides on number with
development's test_1100_advdiff_supg_rotating_gaussian.py. Duplicate
numbers are endemic (70+ pairs already), so this is the house norm rather than
a new breach, but these two are both 1100 and both rotation tests, which is
needlessly confusing. Renumbering is cheap while the branch is being rebuilt.

3. This branch has never run CI — gh pr checks 677 reports "no checks
reported on the 'feature/levelset-supg' branch". So the 2358 lines being kept,
level_set.py included, are unvalidated by CI, and the rebase is the first
time they will be. That is the main risk in the plan, and it is a reason to do
the split rather than resolve 57 hunks: a smaller, CI-exercised branch.

Underworld development team with AI support from Claude Code

[lmoresi]

Split done as a merge, not a rebase — 5237 insertions become 2996

Resolved in 02389614. A merge rather than a rebase, so the branch is not
force-pushed and this review thread survives; the PR diff against
development is nonetheless only the level-set and free-surface work.

src/underworld3/level_set.py	1307
src/underworld3/systems/free_surface.py	+674/−47
Crameri Case1 + Case2	688
LeVeque SUPG-vs-SLCN	178
docs/advanced/level-set-transport.md	93
tests/test_1105_levelset_rotation.py	92
systems/__init__.py, __init__.py, docs/advanced/index.md	11

No rewire was needed. level_set.py already called
systems.AdvDiffusion — development's own name — and
AdvDiffusionSUPG appears nowhere in the kept files.

The two issues, fixed

1. import underworld3 raised ImportError on this branch.
systems/__init__.py had from . import level_set, but there is no
systems/level_set.py — the module is at the package top level. It has to
stay there: the Crameri example reaches uw.level_set.MaterialCLS and
free_surface.py imports underworld3.level_set. So the import path was
what was wrong, not the location. It is now from .. import level_set,
placed last in the file because level_set imports systems itself.

This is why nothing caught it: gh pr checks reported no checks at all
on this branch, so the package was never imported in CI.

Verified against the built library:

uw.level_set                        -> True
uw.level_set.MaterialCLS            -> True
uw.systems.level_set                -> True
uw.systems.LevelSetSolver           -> True
uw.systems.level_set.interface_thickness -> True
uw.systems.AdvDiffusion             -> SNES_AdvectionDiffusion_Composed
uw.systems.AdvDiffusionSUPG         -> gone

2. The test_1100 collision with
test_1100_advdiff_supg_rotating_gaussian.py — both rotation tests. Now
test_1105, which is free on development and on every open branch
(1100-1104, 1110, 1113-1121 are all taken somewhere).

Tests — the first this branch has ever had

• test_1105_levelset_rotation.py — 3 passed, 1 skipped
• test_1070 / test_1072 / test_1074 free surface — 14 passed
• check_test_coverage.py — all 338 files reachable from scripts/test.sh

Two things left for you

peclet_weight defaults to 4.0 on development's solver, a term this
branch's SUPG never had. It is inert here — level_set.py pins
diffusivity = 0.0 and at ν = 0 the factor reduces to ≈1 wherever the
velocity is non-zero — and nothing else in the kept set sets a diffusivity.
Raised because it is development's default for every other caller, chosen
without this branch in view, not because it affects this work.

shapely is not a declared dependency anywhere — not setup.cfg,
pyproject.toml or pixi.toml — so one case in test_1105 skips. It will
skip in CI too, and that case is the polygon initialisation. Either declare
it or accept that the polygon path is untested.

Underworld development team with AI support from Claude Code