Forward-from-nodes history; semi-Lagrangian schemes named by trace and launch; every history equals serial in parallel

[lmoresi]

Stacked on #795.

What changes

One entry point for the semi-Lagrangian histories. uw.systems.ddt.SemiLagrangian(mesh, psi_fn, V_fn, vtype, trace=, launch=) selects one of four schemes:

trace	launch	class (former name)
backward	nodes (default)	BackwardNodesSemiLagrangian (SemiLagrangian)
backward	integration_points	BackwardIntegrationPointsSemiLagrangian (IntegrationPointSemiLagrangian)
forward	integration_points	ForwardIntegrationPointsSemiLagrangian (ForwardSemiLagrangian)
forward	nodes	ForwardNodesSemiLagrangian (new)

A keyword the chosen scheme does not take is a TypeError. The former class names and the former stress_transport strings (semi_lagrangian, integration_point, forward) still work, with a FutureWarning. stress_transport takes backward_nodes (default), backward_integration_points, forward_integration_points, forward_nodes, lagrangian, eulerian; AdvDiffusionSLCN(transport=...) chooses its value history among the four semi-Lagrangian schemes.

Forward from nodes. The history launches its values from the field's nodes and from a lattice inside every element (the discontinuous basis two degrees up), carries them one step forward, fits the arrivals in each cell at the field's degree, and projects the per-cell fits onto the continuous store. A cell with too few arrivals takes a linear fit to its own arrivals, or their mean; an arrival on a shared face is fitted in every cell that contains it. A stress is committed by L2 projection onto the store and launched from that field. Rotating diffusing Gaussian, P2, SLCN solver, integral error:

	half turn, disc	half turn, box	full turn, disc
forward from nodes	1.78e-2	6.57e-2	8.14e-2
backward from nodes	2.52e-2	6.43e-2	8.94e-2

The box, where the flow crosses every wall, is slightly worse than the backward scheme. A serial-only version that let thin cells borrow arrivals from their neighbours reached 1.68e-2 and 6.32e-2; that borrowing depends on the partition and was removed. First order; a fixed mesh.

Every history gives the serial answer in parallel. tests/parallel/test_1066 holds the six stress histories (turned-over Maxwell box) and the four semi-Lagrangian value histories (rotating Gaussian with inflow) to 1e-6 of their serial values at np >= 3; verified at np 3, 4 and 6. Five sources of partition dependence were removed:

1. History projections were solved at the projection default tolerance (1e-4) with an aggregating multigrid; they are now CG with Jacobi to 1e-12 from a zero start (the warm start also made a restart differ from a straight run).
2. The nodal and integration-point histories recorded the current field by evaluating it at nudged nodes (serial) or by copying it (parallel, and only when handed the variable itself); a field that is a mesh variable is now always copied.
3. Behaviour change of the default scheme: the nodal trace-back took its start velocity at the node nudged 0.1% toward "the closest cell's" centroid. A seam node was nudged differently on each rank, and a node on a no-slip wall was given a velocity. The trace now starts at the node. Serial answers near walls change; all 388 tests in the transport, Navier-Stokes and viscoelastic set pass unchanged.
4. global_evaluate: a point that the swarm migration handed to a rank whose cells do not contain it was given the nearest-centroid rbf extrapolation (3e-2 at a steep stress, np 4). The rank whose cell contains it now evaluates it with the FE interpolant. This runs only where the cell hint is authoritative (no DMLocatePoints collective inside), and a finite rbf value is kept where the FE value is NaN.
5. Cell tie-breaks in the forward fits, described above.

Open

• Particle (Lagrangian) stress history differs from serial by 1.2e-5 at np >= 4 #797: the particle history differs from serial by 1.2e-5 at np 4 and 6 (np 3 matches); marked xfail.
• evaluate(..., evalf=True) flags every point as extrapolated #798, evaluate() is off by up to 4e-5 for a field the P2 space holds exactly, at points near a node #799: a false "not located" warning from evalf=True, and serial evaluate off by up to 4e-5 near nodes for an exactly represented P2 field.
• SLCN trace-back gives a partition-dependent answer (1.6% at np=8 on a 128x128 level set) #682 (SLCN level-set volume depends on the rank count) is probably the same set of causes; not rerun.
• TODO(DESIGN): face arrivals are offered to every rank, not only across seams; the two forward schemes use different rules for an arrival exactly on a face.
• Not yet offered for the Navier-Stokes velocity history.

Tests

tests/parallel/test_1066 (new), test_1062 tightened to 1e-6, test_1102 (new, forward from nodes against hard baselines), test_1063 covers forward_nodes, test_1059 covers selection, renamed strings and the Maxwell box for every history. Serial transport/NS/VE set: 388 passed. Parallel np 3/4/6: 10 passed, 1 xfail at np >= 4.

🤖 Generated with Claude Code

https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo

[lmoresi]

Adversarial review

Two review passes before the push (on c434d8a and on b20c992). Findings and disposition:

Fixed

• The SemiLagrangian entry point passed positional arguments to schemes whose positional order differs (varsymbol=2, order=True silently); now keyword-only after vtype, and an argument the scheme does not take is a TypeError.
• Forward from nodes kept the pre-solve fit for a cell with no arrivals, so an empty inlet cell read back the stress from before the solve; the fallback is now the field the cell launched, and nothing is carried between steps.
• Forward-from-nodes inflow detection used the per-rank points_in_domain, which treats partition faces as boundary: every seam node would have taken the inflow value. Now the global restore test the other flavours use.
• The global_evaluate containment round overwrote a finite rbf value with NaN for a located-but-NaN point (the case the fallback exists for), and could meet a DMLocatePoints collective on non-authoritative meshes through the interpolation cache. Now keeps finite values and runs only where the hint is authoritative.
• The owned-node mask relied on a write-then-read-back through the ghost update (wrong inside a deferred update); now read from the global section.
• A strictly-inside arrival could also be kept by a neighbour whose tolerance band it touched (partition-dependent); now only its own cell.
• containing_cells could miss the containing cell beyond the k nearest centroids; now falls back to the locator. Quadrilateral centroids used the simplex formula.
• History projections warm-started from a state a restart does not restore (restart test failed at 1e-12); now zero-start CG/Jacobi to 1e-12.
• test_1102 windows could not see a default change (+-10%); now 1e-4. test_1066 resets the model per case, asserts the seam exchange ran, runs the inflow path; test_1063 covers forward_nodes.

Deferred (TODO(DESIGN) in the code)

• Face arrivals are allgathered to every rank, including faces interior to a partition (no-slip wall nodes each step).
• The forward integration-point exchange gives a face arrival to one cell; forward from nodes to every containing cell.

Filed: #797 (particle history, np >= 4), #798 (evalf extrapolation flag), #799 (serial evaluate near nodes).

Underworld development team with AI support from Claude Code

[lmoresi]

Added (28314aa):

• SLCN trace-back gives a partition-dependent answer (1.6% at np=8 on a 128x128 level set) #682 fixed. global_evaluate now applies the monotone clamp on the rank that evaluates each point, not on the rank that asked. The level-set volume at np 1/2/4/8 is 0.0706771442/415/438/412; np 8 was 1.6% off. Serial results are unchanged.
• NavierStokesSLCN(velocity_transport=...). The velocity history can use any of the four semi-Lagrangian schemes, built by the same helper as AdvDiffusionSLCN(transport=...). test_1103 runs a lid-driven cavity at Re 100: the three schemes that accept P2 agree to 0.2%, and forward from integration points refuses a P2 velocity. test_1066 adds the cavity, which equals serial at np 3/4/6.
• Test runs. The whole tests/parallel suite at np 4: 168 passed, 1 xfailed (Particle (Lagrangian) stress history differs from serial by 1.2e-5 at np >= 4 #797), 1 failed. The failure is test_1063_constrained_freeslip_parallel[ti] (37.7525 against a recorded 37.0780); it fails identically on Viscoelastic stress past the conformation limit: a positive-definite convected step and a log-conformation history #795's 05d5248, so it predates this PR. The serial transport/NS/VE set: 393 passed.

Underworld development team with AI support from Claude Code

[lmoresi]

Added (2adf50d, c0704a8): the SUPG residual sees the carried elastic stress, and NavierStokesSLCN takes any stress history.

SUPG. uw.systems.NavierStokes weighted a momentum residual that did not include the stress divergence. At the exact solution that residual is the elastic stress divergence, the dominant term at high Wi, and SUPG turned it into an O(tau) error. The residual now includes the divergence of the carried stress. That term is the model's flux with the velocity's own derivatives set to zero, with spatially varying parameters expanded so their gradients appear. On a developing Oldroyd-B channel (Re 250, Wi 0.6), the difference from the Galerkin solution on the same mesh fell from 5.0e-5 to 3.2e-5 in velocity and from 3.4e-6 to 1.6e-6 in stress. Fully developed channel flow cannot show this, because the missing residual is constant along streamlines.

NavierStokesSLCN. It ignored stress_transport (the assignment just created an attribute) and always carried a backward nodal history, one step late (#742). It dropped the solvent stress from its momentum flux, forced its momentum order onto the viscoelastic model, and overwrote a supplied DFDt. It now shares one stress-history lifecycle with the Stokes family (_StressHistoryMixin), and both Navier-Stokes solvers form the viscoelastic momentum flux the same way, with the momentum scheme's weights. It reads the log-conformation store: N1 in upper-convected start-up is 1.193 against 1.188 exact, where the linear stress store gives 1.108.

Adversarial review of 2adf50d

Fixed in c0704a8

• NavierStokesSLCN with a viscoelastic model and velocity_transport="backward_integration_points" failed to compile: the solvent term differentiated an integration-point store. Stores without a derivative now give their nodal snapshots through _derivative_stand_ins.
• The theta rule took the stress history's Adams-Moulton weights, which by default gave a BDF2 time derivative against a Crank-Nicolson flux. It now uses the momentum scheme's weights, and a missing stored level is an error rather than a clamp.
• The SUPG memory divergence treated material parameters as constant containers, so the gradient of a spatially varying modulus or viscosity was lost. Non-constant containers are now expanded before differentiating.
• The memory stress was hard-coded to the isotropic form. It is now the model's flux with the velocity's own derivatives removed.
• The new NavierStokesSLCN tests ran a uniform-stress box, so every changed term had zero divergence and the tests could not fail. test_1104 adds a non-uniform developing channel against recorded values (solvent, log store, momentum order 2 against model order 1), the varying-modulus expansion, and an integration-point velocity history.
• The docs contradicted the code on integration-point stores; hasattr branches are replaced by a base _carried_stress_sym; the viscous-flux history is built when the model is assigned.

Open

• The objective-rate terms, the deformation-step terms and a yielding viscosity carry the velocity gradient, so their divergence needs second derivatives. They stay outside the residual, with the viscous term (TODO(DESIGN)).
• The integration-point snapshot stand-in differs from the carried stress by more than O(dt) at an inflow.

Serial transport/NS/VE set: 407 passed. Transport tests at np 3/6 pass; tests/parallel at np 4 is as before (#801 pre-existing, #797 xfail).

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review, two rounds (static reading; findings and what we did)

Round 1, the one-NavierStokes / one-AdvDiffusion commit (78e62a9). 14 findings; the HIGH ones were real defects: AdvDiffusion(transport="backward_integration_points") failed at its default θ = 0.5 (the diffusive flux differentiated the integration-point store; it now reads the snapshot stand-in as the Navier-Stokes solver does, and test_0066's "cn" case is a comparison, not a refusal); the composed AdvDiffusion had no CARRY remesh policy for a nodal trace-back (ported); the FutureWarning said the new class "is" the old one while the defaults, the meaning of order and the form of estimate_dt differ (the warning now says what differs). MED: the SUPG options were silently inert off the Eulerian path (refused, with the Navier-Stokes extrapolation state allocated on the Eulerian path only); V_fn did not reach the characteristic trace; estimate_dt returned inf for a uniform start on a semi-Lagrangian transport (defaults to the cell-crossing time there); a θ change across 1.0 did not rewire the integration-point form; a discontinuous field was accepted. Baselines for the composed solver with each transport (test_1102). Not done: sharing one trace between the velocity and stress histories (two traces, both correct).

Round 2, bf4e663..59366e9 (store smoothing on every history, ParticleL2Projector and reconstruction="global", the degree-2 projector, the deficit fill and bubble penalty, the smoothing-setter rewire fix, FENE-P, round 1's fixes). 22 findings. HIGH, all fixed in fda1286: the FENE-P fields were named by id(self), which differs across ranks, and the first parallel runs deadlocked at their first snapshot (named by a counter now); the projector's build raised on a rank with no cells; FENE-P with the default infinite extensibility was 0/0 (refused). MED, fixed: element= was a dead argument (derived when undeclared, enforced when declared); the forward integration-point global path averaged the previous field at a seam vertex over local cells only (summed across seams now); the forward nodal launch weights summed to 1.3 cell measures (each cell shares its measure among its points and dofs); no KSP convergence check; .coords instead of .coords_nd; fit_limiter settable after reconstruction="global"; the uniform-stress tests could not fail on the multiplicity or the weights (a face-arrival test and a weight-sum test added). LOW fixed: stale docstrings, the restart note, the smoothing check, the utilities export, .array writes, the per-cell overshoot reduced over ranks, a user's remesh policy left alone. Left as they are: _fene_exact_conformation (used by the tests), flux_smoothing (one caller left), the health check reading the exact f where the weak form reads the lagged f* (documented).

Checked and found correct by the reviewer: the P2 edge-to-vertex mapping and basis gradients, the Dunavant rule and bubble mass, the seam assembly with no-overlap distribution, the deficit fill against the inflow top-up, the multiplicity correction, the FENE-P lag's consistency, the rewire fix's structural compare.

Full serial suite after round 2: 2839 passed, 17 skipped, 20 xfailed (one test updated for the clearer ValueError); the stress-transport, projection, FENE-P and composed-solver suites re-run after the last fixes.

Underworld development team with AI support from Claude Code

[lmoresi]

Review of 3170929 (FENE-P: store log(f c), decode the conformation through the record's trace), fixes in 020b02f.

Cause of the FENE-P stall on the cylinder: the record was log c, the nodal projection that commits it overshoots at the wall (tr c 105 with L^2 100 from a stress whose own conformation had 70), the spring factor hit its floor and the velocity multigrid ground to the guard's hang limit. Saturation and the lagged rate were ruled out from the snapshots (f 4.65 at the last one, floor 100).

Findings and what we did:

• The trace decode f = 1 + (tr e^psi - d)/L^2 is positive only for L^2 > d; nothing checked it. The refresh now refuses an extensibility at or below the dimension (test).
• The nodal field the weak form reads was no longer tested after the round-trip test was rewritten. A test writes a record directly and checks the field, the carried strain and the carried stress against the closed form.
• The admissibility assertion was a band (9.7 to 10.0) around a closed-form value (502/51); it pins the value now.
• _peterlin_np had no caller; _fene_spring_factor_of_stress carried a floor and a docstring about the encode that the encode never used. Removed and rewritten as the record formula on sigma/G + I.
• Constructor comment described the old record; the diagnostic read L^2 with float(); docs now say which f* readers keep the lag and that a FENE-P snapshot from the previous record decodes short by f.

Verified correct: encode/decode algebra; every reader of the record (stress_star, E_eff, DEVSS projection, the solvers' decode paths) reads G(e^psi - I); the step algebra is unchanged; inflow encoding; level handling (FENE-P is order 1); units; parallel structure of the refresh.

Shear-box errors against the closed form: +0.0068/+0.0094 at dt 0.1, +0.0034/+0.0051 at 0.05, +0.0018/+0.0027 at 0.025. Cylinder Re 100: Wi 1 res 20 past step 300 (hung at 54 before), Wi 0.5 res 40 past step 100 (stalled at 71 before). All 8 tests in test_1068 pass.

[lmoresi]

Why this is red, and what it needs

The two CI failures are not this PR's. They are byte-identical across #785, #789, #795 and #800:

tests/test_1060_nitsche_freeslip.py::test_nitsche_normal_velocity_zero
    assert np.float64(0.0001008010278702676) < 0.0001
tests/test_1070_free_surface_plume.py::test_freesurface_strong_constraint_beats_penalty
    AssertionError: strong constraint not better: 1.73e-02 vs penalty 3.25e-02

Both tests were already fixed on development, for exactly these numbers, and neither name exists there any more. development's CI is green at 1391e4aa. The replacements name the failures in their own docstrings:

• test_1060 → test_nitsche_constrains_the_wall_normal_velocity, a relative bound of 0.1: "the previous version of this test asserted an ABSOLUTE 1e-4, which scaled with the buoyancy forcing rather than with anything about Nitsche, and sat a fraction of a percent from failing for unrelated reasons."
• test_1070 → test_freesurface_strong_constraint_tracks_the_prescribed_rate, an absolute bound of 5e-2 against the datum: "It replaces a comparison against the penalty path (errors["strong"] < 0.5 * errors["penalty"]) … measured 2026-09-14: 1.06e-2 on macOS, 1.73e-2 on the Linux CI runner, where the annulus triangulates differently."

So this branch is stale, not broken — 58 commits behind development.

The merge is not clean, and the conflicts are the same on all four

git merge origin/development conflicts in exactly two files on every one of the four branches (#800 adds a third, tests/test_0066_integration_point_slcn.py):

scripts/test.sh — one hunk, take THIS branch's side. development narrowed the glob to tests/test_1100*py; this side has tests/test_110*py. test_1101_advdiff_swarm_rotating_gaussian.py matches none of development's globs (test_1100*, test_1110*, test_1120*), so the narrow side would leave a test file dark — #721 recurring. scripts/check_test_coverage.py verifies the resolution.

src/underworld3/systems/ddt.py — five hunks, and four of them are mechanical. Four have a zero-line development side: they are this branch's own additions (applies_inflow_value, commits_flux_in_post_solve, the inflow_value setter — the #745/#783 inflow work), flagged only because surrounding context moved. None of those three names exists on development. Take this branch's side.

The fifth, at the top of _DDtBase, is a real design conflict and needs the author: this branch adds update_exp_coefficients / _exp_alpha / _exp_phi to the base class; development has since put update_exp_coefficients in three per-flavour classes (ddt.py:1299, 1734, 3612) with _exp_alpha/_exp_phi at :3623, 3628, and adds _note_history_shift in that same place. Keeping both gives a base-class definition under three overrides. development's per-flavour layout is the later design, but which one this branch's ETD work is written against is the author's call, not a merge-tool decision — and this is the DDt core, where CLAUDE.md asks for benchmarking rather than a judgement.

Order matters

All four touch ddt.py in the same region, so they conflict with each other as well as with development. They want merging one at a time, smallest first (#789, #785, #795, then #800), with a rebuild and the transport tests between each — the base moves for the rest after every merge.

Found in the backlog sweep (#819). No commits pushed to this branch.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.