fix: swarm empty-rank save deadlock (MPIO dtype) and estimate_dt reshape crash by bknight1 · Pull Request #680 · underworldcode/underworld3

bknight1 · 2026-09-03T07:43:08Z

Summary

Two fixes making passive-swarm operations safe when one or more MPI ranks hold zero particles.

estimate_dt() empty-rank guard — zero-particle velocity data is reshaped to an explicit (0, dim) array instead of crashing on reshape(n, -1) with ValueError: cannot reshape array of size 0 into shape (0,newaxis).
Empty-rank field-dtype preservation — fixes a parallel-HDF5 (MPIO) collective-close deadlock in SwarmVariable.save / write_timestep. An int swarm variable (e.g. add_variable("uid", dtype=int)) held int32 on a non-empty rank but float64 on an empty rank (the empty-rank unpack_raw_data_from_petsc fallback returned a bare float64 np.zeros((0, n))). save() passed the per-rank dtype to the collective create_dataset, so HDF5 metadata diverged and the collective close deadlocked. The field PETSc dtype is now stored (_petsc_dtype) and used for every zero-length fallback so dtypes agree across ranks.

Intended outcome

write_timestep / swarm save no longer hang with empty ranks (was reproducible 6/6 at np=4, now clean 6/6).
estimate_dt() no longer crashes on empty ranks.

The evaluate / global_evaluate / advection empty-rank deadlock was already fixed upstream by #611 / PR #656; coverage here guards that against regression.

Tests

New regression tests (both run with --with-mpi):

tests/parallel/test_0795_swarm_empty_rank_evaluate_save.py — evaluate + write_timestep on empty ranks. test_passive_swarm_save_empty_ranks reproduces the user's crust-only-tracer pattern and hung at np=4 before this change, now passes.
tests/parallel/test_0796_swarm_advection_empty_rank.py — advection on the default non-evalf path with empty ranks.

Result: 5 passed at np=2 and 5 passed at np=4 (4 from test_0795 + 1 from test_0796).

Underworld development team with AI support from Claude Code

On an MPI rank holding zero particles (e.g. crust-only tracers confined to a subset of ranks), the velocity array entering the estimate_dt() max-speed reduction is empty; reshape(n, -1) then raises 'ValueError: cannot reshape array of size 0 into shape (0,newaxis)' because NumPy cannot infer the implied dimension from zero elements. Guard the empty case explicitly so the rank contributes zero to the global max and estimate_dt() returns cleanly. The empty-rank collective point-location deadlock in advection's global_evaluate is fixed upstream (issue #611 / PR #656) and covered by tests/parallel/test_1076_global_evaluate_empty_rank.py, so no local workaround is needed here. Regression: test_0796 exercises the DEFAULT (non-evalf) advection path on empty ranks, which also runs estimate_dt via order=2. test_0795 covers evaluate + write_timestep on empty ranks. Underworld development team with AI support from Claude Code

An int swarm variable (e.g. add_variable('uid', size=1, dtype=int)) holds particles on some ranks but not others. On an empty rank (local_size == 0), SwarmVariable.unpack_raw_data_from_petsc() returned np.zeros((0, n)) which is float64, so an empty rank's data array had dtype float64 while a non-empty rank's was the field's PETSc type (int32). In parallel HDF5 every rank must create an identical dataset, but SwarmVariable.save() passed the per-rank local_data.dtype to the collective create_dataset(); the divergent dtype left the collective close(s) unsynchronised and the save deadlocked — the same silent-hang class as issue #151. Fix: store the field's PETSc dtype (_petsc_dtype) and use it for every zero-length fallback array, so empty-rank data agrees with the field's true type across all ranks. Regression covered by tests/parallel/test_0795_swarm_empty_rank_evaluate_save.py::test_passive_swarm_save_empty_ranks, which hung at np=4 before this change and now passes. Underworld development team with AI support from Claude Code

Copilot

🟡 Changes recommended

The new advection regression test can be flaky/incorrect due to add_particles_with_coordinates() migrating particles across ranks, and the new tests introduce disallowed data-access patterns (mesh.access/.data) per the style charter.

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

Pull request overview

Fixes two empty-MPI-rank failure modes in swarm workflows (advection timestep estimation and parallel-HDF5 saving), and adds MPI regression tests to prevent deadlocks/crashes when some ranks hold zero particles.

Changes:

Preserve PETSc field dtype on empty ranks for SwarmVariable unpack/initialisation, preventing dtype divergence in collective parallel-HDF5 writes.
Guard Swarm.estimate_dt() against zero-length velocity arrays on empty ranks.
Add MPI regression tests covering evaluate/save and advection with empty ranks.

File summaries

File	Description
`src/underworld3/swarm.py`	Stores PETSc dtype on swarm variables for consistent empty-rank fallbacks; adds empty-rank handling in `estimate_dt()`.
`tests/parallel/test_0795_swarm_empty_rank_evaluate_save.py`	New MPI regression coverage for `evaluate()` + `write_timestep()` under empty-rank distributions.
`tests/parallel/test_0796_swarm_advection_empty_rank.py`	New MPI regression coverage for `swarm.advection(..., order=2)` exercising `global_evaluate` and `estimate_dt()` on empty ranks.

Review details

Suppressed comments (2)

tests/parallel/test_0796_swarm_advection_empty_rank.py:47

New code should not use with mesh.access(...) (style charter forbids it); use the variable's .array property directly for data assignment.

    v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=1)
    with mesh.access(v):
        v.data[:, 0] = mesh.X.coords[:, 1]
        v.data[:, 1] = 0.0

tests/parallel/test_0795_swarm_empty_rank_evaluate_save.py:348

New code should not use with mesh.access(...) or the flat .data compatibility layer (see the style charter’s Data Access rules). Use the MeshVariable .array property directly for assignment, and read coordinates from mesh.X.coords.

    with mesh.access(T_mesh):
        T_mesh.data[:, 0] = 300.0 + 100.0 * T_mesh.coords[:, 0] * T_mesh.coords[:, 1]

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

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

        if vel.ndim == 3:
-            vel = vel.reshape(vel.shape[0], -1)
+            # Guard against empty ranks: an array of size 0 cannot be
+            # reshaped with a `-1` axis (NumPy cannot infer the implied
+            # dimension from zero elements) — e.g. (0, 1, dim) -> (0, -1)
+            # raises ValueError. A zero-particle rank legitimately has no
+            # velocities and contributes 0 to the global max below.
+            if vel.size == 0:
+                vel = np.zeros((0, vel.shape[2]) if vel.ndim >= 3 else (0,))
+            else:
+                vel = vel.reshape(vel.shape[0], -1)


+    # Particles only on rank 0 (collective call)
+    if uw.mpi.rank == 0:
+        coords = (np.random.rand(100, mesh.dim) * 0.8 + 0.1)
+    else:
+        coords = np.empty((0, mesh.dim))
+    swarm.add_particles_with_coordinates(coords)
+
+    uw.mpi.comm.barrier()
+
+    # Sanity: the empty-rank precondition must actually hold (all 100 points
+    # are added on rank 0 and no migration happens here, so the other ranks
+    # hold zero particles).
+    sizes = uw.mpi.comm.allgather(swarm.local_size)


+import os
+import numpy as np
+import pytest
+import sympy as sp


+]
+
+
+def test_advection_empty_rank_default(tmp_path_factory):


…rd supersedes the branch's reshape fix Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

…s run with #680 Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

…independent cell size (#687) (#688) * Expose the BDF and Adams-Moulton coefficient symbols on the DDt managers A solver that assembles its own weighted sum of history terms (an Eulerian scheme applying a multistep rule to a spatial operator) needs the constants-routed coefficient expressions, not just their current values. Read-only accessors; no behaviour change. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Pack and index auxiliary fields by DM field, not by position in mesh.vars A MeshVariable that is dropped and garbage-collected (the default Model holds the only strong reference; uw.reset_default_model() releases it, and the statistics helpers delete temporaries deliberately) leaves its PETSc field in the DM. Mesh.update_lvec zipped mesh.vars.values() against the field decomposition by position, and the JIT's petsc_a[] offsets were a running count over the live variables, so every later variable was packed into, and read from, the wrong slots. Measured: 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. update_lvec now packs by field name and zeroes an orphaned field; the JIT reads component offsets from the DM's own field list and patches each variable from its field_id. Regression test: 2 of its 3 checks fail without the fix. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Add the RotatingGaussian transport oracle; fix the integral-norm error for scalar variables A Gaussian carried round the origin by rigid rotation while diffusing is exact at every time (rotation commutes with the Laplacian), so a transport scheme's error can be measured directly and the round trip after one revolution is an absolute check. AnalyticSolution.error(norm='integral') added a 1x1 Matrix symbol to a scalar expression and had never been exercised on a scalar variable. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Extract the per-element timestep estimate shared by the advection-diffusion solvers The cell-crossing / diffusion-time reduction (isotropic or direction-aware, minimum or percentile) becomes a module-level helper so the Eulerian solver can call it rather than carrying a copy. SLCN behaviour unchanged. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Skip the mesh-owned multigrid pickup for a solver that owns its preconditioner A solver with no managed option block (_pc_option_prefix is None) sets its own PC; installing the adapt child's PCMG hierarchy on it segfaulted inside PETSc (additive-Schwarz PC, PCMG calls). The gate now treats that state as the explicit choice it is, alongside preconditioner='gamg' and the user override latch. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Eulerian advection-diffusion with SUPG: BDF and Adams-Moulton orders from the symbolic history uw.systems.AdvDiffusionSUPG(mesh, T, V_fn, order=N, integrator='bdf'|'am') assembles the implicit weak form from the Eulerian DDt history: the BDF stencil or the Adams-Moulton weights on the advective and diffusive terms at every stored time level, plus the SUPG flux tau R u with the strong residual of the same scheme. Timestep, multistep coefficients and the tau weights are runtime constants of the compiled kernels, so a change of dt costs nothing (the issue #657 prototype recompiled on every change). Diffusivity comes from the constitutive model like every scalar solver. Measured on the rotating Gaussian: stable at any cell Courant number, error set by u dt against the feature width (dt^2 for the second-order schemes), unchanged to three digits by a band refined to h/9 at local Courant 13; Crank-Nicolson reproduces the prototype's numbers to four digits. Tests: API and no-recompile contract, temporal convergence (slopes 0.8/0.9 for BDF1, 1.9 for BDF2), band invariance, round trip, np=2 = serial. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note for the Eulerian SUPG solver; BDF2 becomes the default from the integrator study Rotating-Gaussian study at res 32, Courant 0.25 to 8, pure advection and kappa 1e-3: 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, backward Euler carries 20-40% error at any practical timestep. Cost per step is the same for every scheme. BDF2 is the robust default; the note records the alternatives and when to pick them. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Integrator study, res 64: BDF3 grows slowly on pure advection at any Courant number BDF2 and Crank-Nicolson track their res-32 errors at the same u dt. BDF3's stability region misses the imaginary axis near the origin, so the low-frequency modes of a finer mesh grow: 31x the exact field after 590 steps at Courant 1. Safe only with diffusion, below Courant 2. Note and docstring updated; the BDF2 default stands. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * AdvDiffusionSUPG takes the semi-Lagrangian solver's interface: a drop-in replacement The constructor, order, theta, f, V_fn, constitutive_model, delta_t, estimate_dt and solve keep the meaning they have for AdvDiffusionSLCN, so a script changes the class name and 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 documentation gives. The trace-back-only arguments (restore_points_func, monotone_mode, old_frame_traceback, DFDt) are accepted and ignored with a warning. integrator is inferred and only needs setting to reach the higher Adams-Moulton rules. delta_t is settable and solve() reuses it; the notebook viewer reports the scheme. User page docs/advanced/eulerian-advection-diffusion.md with the swap table and the when-to-use-which guidance. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Drop the integrator argument: order and theta already reach every safe scheme The only schemes the argument added were Adams-Moulton at orders 2 and 3, which the integrator study shows blowing up on advection from Courant 1. The multistep family now follows the order (the theta rule at order 1, BDF above); the higher Adams-Moulton assembly stays in the code, reachable only by switching the family on the instance, which is how the study measured it. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * An accuracy-based timestep for the Eulerian solver; credit NengLu in the module and note estimate_dt now returns the step at which the field changes by a fraction (0.02) of its range: from the advective rate |u . grad phi| at the vertices before the first solve, and from the rate the last step actually produced after it. The cell-crossing time the semi-Lagrangian solver reports is not a stability limit for this scheme and says nothing about its accuracy; it stays available as basis='resolution'. The estimate is mesh-independent, which the band test now checks (the resolution estimate collapses 3x on the refined child, the accuracy estimate moves under 25%), and at the default fraction Crank-Nicolson completes the rotating-Gaussian round trip under one per cent. The advective rate uses the vertex Clement gradient rather than a point evaluation of a derivative expression, which fails on a mesh carrying many variables. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Match the Krylov tolerance to the SNES tolerance, and make preconditioner="fmg" a real switch on the SUPG solver The Eulerian SUPG step took two Newton iterations on a linear operator: the Krylov default (rtol 1e-5) does not reach the SNES tolerance (1e-8), and the second Jacobian assembly cost more than every linear solve of the step. The Krylov tolerance is now 1e-9 and a step is one Newton iteration: 1.54 s to 0.91 s per step at 256^2 in serial. Measured against geometric multigrid at matched tolerances (design note, "Preconditioner"), GMRES with additive-Schwarz ILU is the cheaper linear solve at every Courant number from 1/2 to 32 and its iteration count is the same on one and eight ranks; the multigrid's cycle count grows with the Courant number nearly as fast, and a cycle costs about three Schwarz iterations. Schwarz stays the default on every mesh. preconditioner = "fmg" now hands the block to the managed multigrid route (custom-P transfers over the refinement hierarchy or an adapt child's coarse tail, flexible GMRES outside) for the rank count where a one-level method runs out of coarse space. The solver's solve() builds through the base _build, where a preconditioner choice is resolved; the pre-run of the three setup stages marked the solver set up first, so the request was silently inert. The semi-Lagrangian solvers share that pattern and the defect (#683). Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the 512^2 rows at matched tolerance * Let theta be set after construction, as the semi-Lagrangian solver allows The shipped convection examples set adv_diff.theta = 0.5 after building the solver; the Eulerian drop-in refused it. The blend is a runtime constant refreshed from the history manager before every solve, so the setter updates it without a recompile (order 1 only, the constructor's rule). Vector and tensor unknowns join the design note's deferred list: the solver is scalar, where the semi-Lagrangian trace-back carries them. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Navier-Stokes with Eulerian SUPG momentum transport, and a partition-independent cell size uw.systems.NavierStokesSUPG: the incompressible Navier-Stokes equations on the Stokes saddle-point solver with the momentum advection assembled implicitly and stabilised by the vector SUPG term F1 = tau R (x) a, the counterpart of the scalar Eulerian solver. Crank-Nicolson at order 1, BDF2 at order 2, with the velocity history on the mesh; no stress history, the viscous stress at an earlier level is rebuilt from the stored velocity through the constitutive model. The advecting velocity is a choice: the second-order extrapolation 2u^n - u^{n-1} (one linear solve per step, the default), Picard passes on the latest iterate, or the unknown itself under Newton. The strong residual the SUPG term sees carries the pressure gradient; without it the term is O(1) at the exact solution and costs fifty times the Galerkin error on Kovasznay flow. mesh.cell_size() now reports each cell's own radius, the RMS distance of its vertices from its own centroid, taken from the DM's coordinates. The kd-tree radius it used to copy picks the nearest centroid among the rank's cells, so the field differed with the partition (#687, found because the two-rank Navier-Stokes answer differed from serial by 5e-4 and matched to 1e-15 with a constant h); after a deform it also read stale vertex coordinates against fresh centroids. get_min_radius and the other consumers of the kd-tree radii are unchanged. Tests: the solver's API contract (construction rules, one linear solve per step, Picard passes, the Stokes limit, runtime-constant timestep and theta), a two-rank Kovasznay error that matches serial to 1e-7, the scalar parallel reference re-recorded for the new cell size, and the Nitsche local-h tests reading the field's definition. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Document the SUPG Navier-Stokes solver: user page and the design-note section with Kovasznay and cavity results Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: cavity Courant 1 row, cylinder wake rows, and the corrected Re 1000 status * Design note: the semi-Lagrangian cylinder row * Design note: Re 1000 cavity rows, the finer cylinder mesh, and the Galerkin control that cannot run * Swarm.advection: let estimate_dt see a rank that holds no particles (#693) The velocity evaluated for the timestep estimate has shape (0, 1, dim) on an empty rank, and reshape(0, -1) cannot infer the trailing size; the empty-rank handling a few lines below never ran. Give reshape the size explicitly. Found with passive tracers released at the inlet of the DFG cylinder on four ranks, where every rank but the inlet's is empty at the first step. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Set the PETSc constants on the DS the integrals use, so expression values reach the kernels (#695) uw.maths.Integral, BdIntegral and CellWiseIntegral compile their integrands through the same JIT as the solvers, which routes every uw.function.expression to PETSc's constants array, but none of them ever called PetscDSSetConstants: the kernels read zeros, so any integrand with a viscosity, a time or another expression in it integrated to nothing, and a fresh Integral returned the same zero from the cache. Found on the DFG cylinder drag, where the viscous traction (eta is an expression) vanished and the drag read 23 to 28% low on two meshes without moving with the SUPG weights. Each class now packs the manifest and sets the constants right after the objective; the boundary integral sets them on its sandbox DS, which has its own discrete system. Regression test test_0503 covers the three classes, a changed value without recompilation, and the constitutive-flux traction that found it. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the cylinder drag was the missing viscous traction (#695), the vortex-decay benchmark, and #696 The DFG cylinder section is rewritten around what the tau sweep found: the stabilisation moves the drag by 2.6% and the deficit was the boundary integral dropping the viscous part (#695). With the integrals fixed and only the cylinder cells refined through gmsh at a fixed time step, drag, pressure difference and Strouhal number converge onto the reference bands on the 1/20 channel mesh, the traction and reaction measurements close on each other, and the whole-mesh 1/40 run buys less than the 1/320 cylinder cells do. The Galerkin form that "could not run" was the GAMG fallback; the refinement callback gives FMG on the gmsh mesh. New Taylor-Green vortex-decay subsection (dt and h sweeps for CN and BDF2, Galerkin against SUPG, the viscosity range, the advecting-velocity choices), and the two defects it found: #695 and the zero-valued expression folding (#696, raised, not patched). Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the FMG rows of the cylinder table (base mesh refined through the circle callback) Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the Picard row at 1/640 cylinder cells settles the lift overshoot as the extrapolation lag Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Take swarm.py from development (#680): the empty-rank estimate_dt guard supersedes the branch's reshape fix Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: LU is serial-only on the velocity block; parallel tracers run with #680 Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the FMG cylinder-refinement table (Picard, Newton, BDF2, four ranks) Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * NavierStokesSUPG: an opt-in recovered viscous term in the SUPG residual recovered_viscous=True projects the deviatoric stress of the advecting velocity onto a continuous symmetric tensor before each solve pass and puts its divergence in the strong residual the SUPG term sees. Without it the residual lacks the viscous term (the kernels see first derivatives only), an O(h^2) inconsistency for P2 velocity that shows on resolved viscous flow: four times the Galerkin error on the 1/64 vortex-decay mesh, sixteen times on Kovasznay at 1/32. The projection's function is set on first use (it needs the constitutive model) and its default is a zero matrix, not None. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * NavierStokesSUPG: the recovered viscous term is the previous level's momentum balance The differentiated projection of the stress was unstable (1/64 vortex decay and Kovasznay at 1/32 blew up) and did nothing at 1/32. Louis's form: the momentum balance of the stored level gives div sigma^n = rho (Du/Dt)^n + grad p^n - f from first derivatives of stored fields, so the residual the SUPG term weights becomes the increment of the out-of-balance force between levels, at the cost of one stored pressure level and no extra solve. At a discrete steady state that residual vanishes and the stabilisation switches off, which is a property to measure, not assume. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the recovered viscous term measured (balance form = Galerkin accuracy on resolved flow, unstable on the cylinder) Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * NavierStokesSUPG: recovered_smoothing projects the balance term with a screened-Poisson length The plain balance term is unstable where advection dominates because it carries the grid-scale residual of the previous step. With a smoothing length the term is projected onto a continuous vector field (one vector projection per step), keeping the smooth viscous divergence and filtering the rest; zero keeps the plain form. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Revert "NavierStokesSUPG: recovered_smoothing projects the balance term with a screened-Poisson length" This reverts commit 19edda1. * Revert "NavierStokesSUPG: the recovered viscous term is the previous level's momentum balance" This reverts commit c5c72ec. * Revert "NavierStokesSUPG: an opt-in recovered viscous term in the SUPG residual" This reverts commit 42aaedf. * Design note: the recovered viscous term measured three ways and withdrawn Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * NavierStokesSUPG: tau_shape selects the Brooks-Hughes or doubly asymptotic parameter The inverse-sum tau (Shakib-Tezduyar) is above the optimal 1-D curve at cell Peclet numbers of order 1 to 10, where the resolved benchmarks sit. The optimal shape tau = (h/2|a|)(coth Pe - 1/Pe) and its two-limit approximation (h/2|a|) min(Pe/3, 1) are now selectable, each combined with the transient term so the time step still caps them. coth is written through tanh: the C printer rewrites coth through exp and drags the square root in |a| into exp(log(.)), which brings arg() into the kernel. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the shape of tau measured (Brooks-Hughes, doubly asymptotic) Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the Re 1000 cavity converged (94 to 96% of Ghia at 1/64); the rank-local v_max explained Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * NavierStokesSUPG: peclet_weight turns the SUPG term off where the cell is diffusion-dominated The term is multiplied by Pe^2 / (Pe^2 + Pe_c^2) with Pe the cell Peclet number of the advecting velocity, so it is absent where it is not needed (where it costs a fixed multiple of the Galerkin error) and full where advection dominates. Zero (default) leaves the weight uniform. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the weight by cell Peclet number measured (Galerkin accuracy where resolved, stabilisation kept on the cylinder) Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * NavierStokesSUPG: the cell-Peclet weight is the default (Pe_c = 4) Louis's ruling, the code being unreleased: the SUPG term is weighted by Pe^2 / (Pe^2 + 16) by default, off where a cell is diffusion-dominated and full where advection dominates. Kovasznay at 1/8 (the parallel test's reference) goes from 3.83e-3 to 1.42e-3; the design note's earlier tables were made at the uniform weight and say so. The scalar transport solver keeps the uniform weight until its convection benchmarks are re-measured. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * AdvDiffusionSUPG: the cell-Peclet weight, as for the Navier-Stokes solver (Pe_c = 4) Written without dividing by kappa, so pure advection (the default kappa = 0) keeps the uniform weight and its tests do not move. The convection benchmarks are re-measured with it in the design note. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Docs: the cell-Peclet weight on the scalar solver's user page Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Design note: the cell-Peclet weight is the default of both solvers; the convection rows with it Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Examples: the SUPG Navier-Stokes solver on the lid-driven cavity and the Taylor-Green vortex Two runnable examples in the repository's format: the cavity at Re 100 against Ghia (about four minutes) and the Taylor-Green vortex decay with its exact error and energy decay (about a minute). Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * The DDt history manager is the transport plugin: EulerianSUPG assembles advection and SUPG, the solvers compose A solver that owns an unknown now composes its residual from three contributions of its DuDt (time_derivative, advection, stabilisation_flux) plus the levels and weights of the scheme (states, spatial_weights), and never asks which flavour it holds. The new ddt.EulerianSUPG assembles the implicit advection component-wise for a scalar, vector or tensor unknown and the SUPG flux tau R (x) a of the solver's strong residual; the history-carrying flavours answer zero for both. V_fn is data on the manager (V_fn_history names the carrier of the stored levels, the stored velocity for momentum), the timestep is a runtime constant every flavour writes (delta_t), and the stabilisation knobs live on the manager with the solvers' properties passing through. AdvDiffusionSUPG and NavierStokesSUPG lose their own residual code and compose the same way. A SemiLagrangian manager dropped into the scalar solver reproduces AdvDiffusionSLCN on pure advection; a flattened symmetric tensor is transported through SNES_MultiComponent with a residual that is only the manager's terms (test_1057). The plain Eulerian manager keeps its explicit splitting correction behind an _advection_mode gate and gains num_components for MATRIX histories. Regression: Kovasznay 1/16 and 1/32, the vortex decay at 1/32, the Blankenbach box and both examples reproduce their recorded numbers to every printed digit; the cylinder keeps its mean drag, lift extrema and Strouhal number, with the drag peak moving 3.0797 -> 3.0802 (evaluation order in a shedding wake). Two-rank tests keep their serial constants. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * The composing solvers take the generic names: AdvDiffusion and NavierStokes; the semi-Lagrangian classes keep their SLCN names A solver that composes its transport from its DDt manager is not an SUPG solver: SUPG is a property of the EulerianSUPG manager it holds by default, and a SemiLagrangian manager makes the same solver a semi-Lagrangian scheme. So uw.systems.AdvDiffusion and uw.systems.NavierStokes now name the composing solvers (SNES_AdvectionDiffusion_Composed, SNES_NavierStokes_Composed) and the SUPG class names are gone. The semi-Lagrangian classes stay reachable as AdvDiffusionSLCN, NavierStokesSLCN and NavierStokesSwarm; every existing use of the generic names with the semi-Lagrangian meaning in docs, notebooks, examples and tests is moved to the explicit SLCN name, so nothing changes scheme. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Restore the branch's mesh changes the whole-file merge resolution dropped; Picard reductions on every pass The merge of development took discretisation_mesh.py and test_1065 wholesale from development, losing the orphaned-field packing by name (test_1058) and the rest of the branch's non-conflicting edits; this is the hunk-by-hunk resolution with the landed cell_size (#692). The Picard loop of the Navier-Stokes solver now takes its two reductions on every pass, and the break predicate is recorded as rank-uniform in the collective-guard scan. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Parallel tests: a refinement hierarchy for the Navier-Stokes reference, platform-tolerant comparison to the serial error The GAMG fallback on a mesh without a hierarchy gave a platform-dependent answer (7% on the Linux CI); the test now refines a 1/4 mesh once so the velocity block runs geometric multigrid, and the serial reference (0.00132279) is met by two and four ranks to 3e-10. Both tests compare to the serial error at 1e-6 relative: the partition effect they guard against was 5e-4 (#687), platforms differ at 1e-7. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Adversarial review of the plugin and the rename: six fixes From three reviews of the branch head (findings posted on #688): - the base contract's shape helper collided with Symbolic's `_shape` attribute, so Symbolic.advection() raised instead of answering zero; renamed; - a user-supplied EulerianSUPG on NavierStokes advected the stored level with the new velocity: the solver now sets V_fn and V_fn_history whoever built the manager, and its advection setter only steers such a manager; - the change-rate bookkeeping read the manager's history `.array`, which fails for a SemiLagrangian history under units and for a swarm-backed history; it now diffs a copy of the unknown's data; - a supplied manager silently overrode `order`/`theta`; a mismatch is an error, and the theta setter refuses a manager without theta; - the 1-D tau shapes divided by the diffusivity (zoo at the manager's default); - the timestep and SUPG knobs are created with unique names like the BDF coefficients, so they do not accumulate in the persistent registry; - a bare scalar residual is accepted by stabilisation_flux. Rename loose ends: an example that imported the bare NavierStokes name now uses NavierStokesSLCN explicitly; tutorial 9 prose; the solver-unification design table; API entries for the composing classes and the manager. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * DDt: a quantity timestep is non-dimensionalised before it reaches the kernels (#701) _as_float took the magnitude of a Pint or UW quantity, so a semi-Lagrangian solver stepped with 100 kyr under a 1 Myr reference time wrote 100 (not 0.1) into the manager's runtime timestep and into the variable-step BDF bookkeeping. It now goes through uw.non_dimensionalise, which handles both quantity types; without reference scales the magnitude is what comes back. Test with a negative control in test_1057. Closes #701. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Swarm.advection: a clear error for a swarm that was never populated (#702) DMSwarm reports a local size of -1 until particles are added on some rank, and the advection then failed inside numpy with 'negative dimensions are not allowed'. The empty rank of a populated swarm (size 0) is unchanged. Closes #702. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * UWexpression never reports is_zero, is_positive or is_negative from its value (#696) A runtime constant's value can change after construction, so sympy must not fold on its current sign or on it being zero: exp(c) with c created at 0 evaluated to 1 at construction and a time ramp that started at t = 0 stayed frozen (found on the Taylor-Green Dirichlet case). The three assumptions now answer None, as for a plain Symbol; the value is read when the expression is unwrapped for compilation. Control in test_0503: exp(c) survives, integrates to 1 at c = 0 and to e at c = 1. Level-1 suite: 1704 passed. Closes #696. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL * Copilot review of #688: NavierStokes.estimate_dt dimensionalises the accuracy estimate; EulerianSUPG takes no mutable default bcs The accuracy basis returned a bare non-dimensional number while the resolution fallback returns a quantity under a scaling model; both now come back through _dimensionalise_dt (test under reference scales). The manager's bcs default is None -> a fresh list; a caller's list is still kept by reference on purpose, so a solver's live essential_bcs reach the projections. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --------- Co-authored-by: Claude Fable 5.1 <noreply@anthropic.com>

bknight1 added 2 commits September 3, 2026 15:39

bknight1 requested a review from lmoresi as a code owner September 3, 2026 07:43

Copilot AI lite review requested due to automatic review settings September 3, 2026 07:43

Copilot started reviewing on behalf of bknight1 September 3, 2026 07:44 View session

Copilot AI reviewed Sep 3, 2026

View reviewed changes

Merge branch 'development' into bugfix/swarm-empty-rank-save-reshape

6d3a358

lmoresi merged commit 37308de into development Sep 4, 2026
2 checks passed

This was referenced Sep 6, 2026

Navier-Stokes with Eulerian SUPG momentum transport, and a partition-independent cell size (#687) #688

Merged

Advecting a swarm that was never populated fails with 'negative dimensions are not allowed' (DMSwarm local size -1) #702

Closed

Navigation Menu

Uh oh!

fix: swarm empty-rank save deadlock (MPIO dtype) and estimate_dt reshape crash - #680

fix: swarm empty-rank save deadlock (MPIO dtype) and estimate_dt reshape crash#680
lmoresi merged 3 commits into
developmentfrom
bugfix/swarm-empty-rank-save-reshape

bknight1 commented Sep 3, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

		]


		def test_advection_empty_rank_default(tmp_path_factory):

Uh oh!

Conversation

bknight1 commented Sep 3, 2026

Summary

Intended outcome

Tests

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🟡 Changes recommended

Suppressed comments (2)

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants