Skip to content

Stress transport on the grid, forward from launch points, and on particles - #833

Merged
lmoresi merged 33 commits into
developmentfrom
feature/forward-parallel
Oct 8, 2026
Merged

lmoresi merged 33 commits into
developmentfrom
feature/forward-parallel

Conversation

@lmoresi

@lmoresi lmoresi commented Oct 7, 2026

Copy link
Copy Markdown
Member

The branch four open PRs are stacked on, and the reason none of them can land: it has had no PR of its own. #785 → #789 → #795 → #800 all base here, so they have nothing to merge into.

32 commits, +3294/-278 across 16 files, of which 1,159 lines are tests.

What it adds

A stress history that moves on the grid. stress_transport on a Stokes or Navier-Stokes solve selects where the viscoelastic stress is carried: the Eulerian SUPG manager transports its own history on the grid, the integration-point history carries the stress instead of rebuilding it, and the theta rule reads the history for its stored flux level. The history manager commits its own flux and shifts its own levels, so the solver no longer has to know which flavour it holds.

Two more transports. stress_transport = "forward" carries the stress from launch points inside the cells, and "lagrangian" makes the particle history an equal citizen. Both are snapshotted and restored, and the forward one runs in parallel.

An exponential integrator on every history. integrator='etd' was semi-Lagrangian-only (#739); the coefficients are now refreshed before the first carry (#740, #741), and the integration-point flavour carries a strain-rate history along the characteristic for ETD-2.

An inflow condition. Applied by the sign of the normal velocity, living on the base history so no flavour drops it silently (#733), with a true inflow value at restored departure points (#745).

DEVSS, opt-in, on the Stokes momentum flux for a discontinuous elastic stress, with the pair cancelling on a plain Stokes solve too (#754).

AdvDiffusionSwarm, a swarm (Lagrangian) advection-diffusion solver, its history proxy defaulting to degree 1 so it holds a boundary-crossing flow.

Issues this closes

Eight of the eleven issues its commits name are already closed — #727, #732, #733, #739, #740, #741, #745, #754. Two more close with it:

Closes #737. The integration-point history stores the stress by a global L2 projection and samples it back at the quadrature points, and that cycle has no dissipation at the cell scale: below Courant one a cell-scale mode grows from round-off, measured at 2.4 per unit time on the Maxwell Waters-King start-up. Removing the element null space does nothing, private nodes per cell diverge, DEVSS does not touch it. store_smoothing puts a Laplacian in the store projection with alpha = c * cell_size^2 — a field, so the dose follows the local cell; c = 0.07 holds the 1/32 case unconditionally at 0.5% on the peak, and the irregular mesh needs that value.

Closes #768. On the confined cylinder every history lost positive-definiteness of tau*/G + I in the first step at Courant one on the far-field mesh, because the wall shear rate is ten times the far-field one and the explicit stretching term cannot span dt*gammadot ~ 3. The viscoelastic model gains max_elastic_timestep(safety) and conformation_min_eigenvalue() — the health line that tells a lost preconditioner from a lost problem — and both trace-back histories expose carried_tensors() for it.

Issues in its territory it does NOT close

issue verdict
#742 psi_star[0] holds the new stress on Stokes and the previous step's on NavierStokesSLCN still live — the branch touches navier_stokes_eulerian.py, but no commit addresses the convention difference, and nothing here reconciles it
#749 Eulerian stress history blows up at a shear-wave reflection (Waters & King) still live — the grid history is what this branch adds, and this is its known failure; not fixed here
#735 nodal semi-Lagrangian carries ~50% too much stress at a no-slip wall still live
#738 Lagrangian_Swarm gives a NaN first residual on the VE cylinder still live — stress_transport = "lagrangian" is added here, and this is its open defect on that benchmark
#797, #783, #784, #788 carried by the stack above (#785, #789); labelled fixed-in-PR
#811 forward per-cell fit unstable where the flow empties a cell still live — against the "forward" transport this branch introduces
#722 flat .data writes against Charter §7 still live, and partly from this work — relabelled enhancement, it is a conformance chore
#682, #683, #685, #712, #402 still live, not this branch's territory

Worth saying plainly: this branch adds three transports and leaves each one's known defect open behind it (#749 grid, #811 forward, #738 lagrangian). That is deliberate rather than an oversight — the transports are usable and measured, and each defect is a bounded case recorded against the scheme it belongs to.

Tests

tests/test_1059_stress_transport.py (740 lines), test_1060_stress_store_smoothing.py, test_1061_stress_forward_history.py, test_1063_stress_history_restart.py, test_0074_return_to_bounds_on_a_file_mesh.py, test_1101_advdiff_swarm_rotating_gaussian.py, and tests/parallel/test_1062_forward_stress_history_mpi.py.

Hard baselines rather than method comparisons: one BDF-1 step from rest at dt = lambda gives conformation eigenvalues 1 -+ 1/2 exactly; a wall four times faster gives -1 everywhere. The ringing test at 1/16, dt 0.0125 asserts the measured growth and its suppression, and reads the content after the trace-back, since after the store it is zero by construction.

A branch-level adversarial review was applied before this PR, as its own commit.

For the merge

Four PRs are stacked on this branch. Merging it lets #785 retarget to development, and the rest follow in order (#789, #795, #800), each wanting a re-merge from development between — they touch the same region of ddt.py.

Underworld development team with AI support from Claude Code

lmoresi and others added 30 commits September 9, 2026 11:21
The viscoelastic Stokes solve reached into its DFDt to project the new stress
into level 0, fan the flat multi-component result back into the tensor, and
shift the levels. That is the manager's business, and doing it through
commit_flux_to_history() is what lets a solver hold a different flavour of
history without knowing which. Behaviour unchanged: tests/test_1051_VE_shear_box
(the analytic Maxwell shear box) and the constitutive and snapshot regressions pass.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
A viscoelastic solve carries a stress that is not its own unknown, so the
manager must move it: transport_on_update makes update_pre_solve carry every
stored level forward by one implicit SUPG step, in place of a semi-Lagrangian
trace-back. The independent components are flattened onto one matrix variable
and marched together through a multi-component solve, backward Euler in time,
with the manager's own advecting velocity and stabilisation parameter.

tests/test_1059_stress_transport.py is the transport test the existing
viscoelastic benchmarks cannot give us: they are all spatially uniform and so
transport nothing. Uniform translation of a stress blob against the exact
answer and against the semi-Lagrangian history on the same case; components
carried round a rigid rotation; and the transport staying off for a manager
that assembles its advection in a solver's residual instead.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…o managers in the test

Backward Euler on the transport step made it first order and eight times worse
than the semi-Lagrangian trace-back on uniform translation; the theta rule (the
manager's, Crank-Nicolson by default) brings it to 0.0139 against 0.0017. The
test shared one stress variable between the two managers, which coupled the
transports it was comparing: each now has its own.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
The transport step's theta was a hidden constant whose docstring claimed it
followed the manager's; it is its own knob, Crank-Nicolson by default, and the
docstring now records why lowering it to backward Euler is not a real option.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…transport

The solver names the flavour that carries its stress history, defaulting to the
semi-Lagrangian trace-back. To let the Eulerian manager stand in its place the
snapshot machinery that keeps the flux projection one-shot moved from the
semi-Lagrangian class onto the base, the Eulerian projection now goes through
the same source builder and follows psi_fn, and projection setup is idempotent.

A history committed this step must not be shifted again by the post-solve: the
Eulerian post-solve shifted the levels a second time, pushing the new stress
straight into level 1. Invisible at order 1 and a 150-fold error at order 2.
With that fixed the two flavours agree to six digits on the analytic Maxwell
shear box at both orders (1.54e-2 at order 1, 9.33e-4 at order 2, both against
the exact solution), which is the check the uniform benchmarks can give: the
stress does not move there, so any difference would be the plumbing.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
A shear box whose shear modulus varies in x: the stress is non-uniform and the
shear carries it, so the transport term is active. No closed form, so the
schemes are judged against each other and under refinement. They agree to
1.0e-3 at order 2 and 5.5e-3 at order 1, in both cases the size of the shift
from halving the timestep, so the gap between them is discretisation and not a
defect. The grid history is the cheaper of the two here.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
Assigning a viscoelastic constitutive model to the Eulerian Navier-Stokes
solver now creates a stress history, as it does on Stokes, and stress_transport
chooses the flavour. The solver drives the viscoelastic sequence itself, once
per step around its Picard passes rather than once per pass: the sequence is
extracted from the Stokes solve into _stress_history_pre_solve and
_stress_history_post_solve, and the momentum passes ask the Stokes solve to
leave the history alone. A stress history needs order 2 here, because the theta
rule weights the viscous flux at the stored velocity levels and would rebuild
them blind to elasticity.

The two flavours agree to 3e-5 on a sheared box with inertia. That case has no
closed form (at this effective viscosity the velocity is far from steady simple
shear within the run), so it is a consistency check, not an accuracy one.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
Crank-Nicolson weights the momentum flux across time levels. For a viscous
fluid the stored level is rebuilt from the stored velocity as 2 eta strain(u);
for a viscoelastic one that rebuild is wrong, and the object the theta rule
wants is already held by the stress history. The two histories meet there: the
flux history of the time scheme IS the elastic stress history.

That removes the order=2 restriction I put on a viscoelastic Navier-Stokes
step, which was a consequence of the viscous rebuild and not of the physics.
Order 1 is the scheme the Navier-Stokes benchmarks use and the default. Both
stress transports agree to 5e-6 at Crank-Nicolson with a first-order stress
stencil, and to 4.4e-5 with a second-order one.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…ves the order ramp (#727)

The stress history's effective order ramps over the opening steps, and when it
changes the compiled functions must be rewired. Setting that flag after the
build had already decided whether to set up left the managed multigrid block
asking a preconditioner for sub-solvers it had not created, and PETSc refused:
the first step was healthy and the run died around step 13. The preparation
(elastic timestep, order check, flux expression) now runs before the build and
the advance (carrying the history, refreshing coefficients) after it, for both
the Stokes and the Navier-Stokes paths.

Verified on the DFG cylinder at Weissenberg 0.3 with a coarse base refined once
to the same resolution the direct solve used: no state error, two Krylov
iterations per step.

Closes #727.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…ormal velocity

A transported quantity needs data wherever the flow enters and nowhere else, and
on a shedding wake those places move: the outflow boundary of the DFG cylinder
carries reversed flow over 20 to 40 per cent of its nodes, migrating as the wake
sheds. Left unconstrained, whatever the solve produces there is injected and
carried upstream; measured, the stress maximum leaves the cylinder for the outlet
at step 15 and the run is destroyed by step 33.

 adds the condition weakly on every boundary through
the negative part of u.n, so it is active exactly where the flow enters, vanishes
where it leaves, and contributes nothing at a wall. With it the same case runs to
completion with the traction outlet in place: drag 2.370 against 2.368 for a
prescribed outflow, peak stress 0.326, and no reversed flow at the outlet at all
-- the injected stress was driving the reversal that admitted it.

The sign is not obvious and is recorded with its measurement: the other way round
the term is anti-dissipative and the stress reaches 16 within ten steps.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
development added a guard so a particle-carried stress history skips the nodal
projection and shift, which it does itself in its post-solve. This branch had
moved that block onto the history manager; the two are reconciled by a
commits_flux_in_post_solve flag the manager declares, so the solver still does
not ask which flavour it holds.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…g it (#732)

IntegrationPointSemiLagrangian.update_pre_solve swallowed store_result into
**_ignored. The viscoelastic path passes store_result=False so the history is
not re-recorded from psi_fn, which for a viscoelastic solver is the
constitutive flux and so a function of the history itself; SemiLagrangian
honours it and this flavour did not. Every step it shifted its snapshots and
re-recorded the flux over a field that was already the stress, applying the
constitutive update a second time and discarding what commit_flux_to_history
had placed. On the analytic Maxwell shear box at order 1 that overshot the
relaxed stress by 12.8% where the nodal and grid flavours sit at 1.5%.

Honour store_result, and commit the flux to the nodal snapshot as well as to
the integration points. The two ladders have to move together: the assembler
reads psi_star, but the next trace-back samples psi_snap, so a flux committed
only to psi_star does not survive a step. The commit projects rather than
evaluates, because the trace-back samples the snapshot between its nodes and
needs the L2 fit of the flux on the history space.

Maxwell shear box, sigma_xy at the origin, all three stress_transport options:

  order 1   0.851356   0.851356   0.851356   (exact 0.864665, 1.54%)
  order 2   0.863858   0.863858   0.863858   (exact 0.864665, 0.09%)

identical to six digits, where the integration-point flavour read 0.975696 at
order 1 before. Order 2 also exercises the snapshot ladder: a level shifted
twice shows up there as a 150-fold error, as it did for the Eulerian flavour.

On the varying-modulus shear box, where the stress is non-uniform and the shear
genuinely carries it, the three norms are 0.873247 / 0.872329 / 0.872064, a
0.14% spread against a 0.4% self-drift when the step is halved.

The test over the Maxwell box now runs all three flavours and judges each
against the analytic relaxation curve; the tolerance it already carried is what
catches this defect.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…ence (#733)

The property was defined only on EulerianSUPG. On any other flavour the
assignment landed as a plain instance attribute and nothing read it: a driver
that sets ns.DFDt.inflow_value whatever stress_transport names got the
condition on one flavour and nothing on the other two, with no error and no
warning. On the viscoelastic cylinder the inflow condition is the difference
between a run that completes and one that dies at step 33.

The property now sits on _DDtBase with the shape check. EulerianSUPG declares
applies_inflow_value = True and keeps a setter override that drops its compiled
transport solver. Every other flavour accepts and stores the value and warns
once that it does not apply it: a trace-back or particle flavour restores an
out-of-bounds departure point to the boundary and reads the transported field
there, which constrains the inflow but is not the value the caller set. Wiring
it in as a true out-of-bounds value is a separate 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
…c stress

With a transported elastic stress the momentum flux is 2 eta_eff edot(u) +
c sigma*, and as the Weissenberg number rises the elastic term dominates: the
velocity is driven by the divergence of a field coarser than edot(u),
discontinuous across elements, and not derived from the velocity. The
discrete operator loses its viscous character and the stress's inter-element
jumps force the velocity at mesh scale with nothing to damp them. Measured on
the viscoelastic cylinder with an integration-point stress history: 14% of the
stress by rms is inter-element jump content, every step, and the near wake
scallops at cell scale before shedding begins.

DEVSS (Guenette & Fortin 1995) adds and subtracts a viscous term,
F1 += 2 eta_a (edot(u) - D), with D the strain rate projected onto the stress
history's space and lagged one step. At convergence the pair cancels to
projection error, so nothing physical is added; but the first term is implicit
in the velocity and the second is data, so a mesh-scale velocity response --
which the projection does not carry -- sees the full viscosity eta_a while
smooth modes see none. It acts on the response to the discontinuous stress,
not on the stress: the history keeps everything it carries.

Stokes.devss_viscosity = eta_a switches it on (None, the default, is off and
the flux is unchanged). D lives in the stress history's space (degree u-1,
continuous), is set from the initial velocity at first setup, and refreshed
after every solve. The Navier-Stokes SUPG flux takes the same term.

Test: on the uniform Maxwell box D equals edot exactly and the answer does not
move for eta_a = 1 (1e-8 of the exact value); on the varying-modulus box the
term is live and moves the answer by projection error only.

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='etd' failed with AttributeError on the eulerian and
integration_point stress histories: the _exp_alpha / _exp_phi accessors
and update_exp_coefficients lived on SemiLagrangian only, and the
integration-point flavour did not allocate the coefficients at all.
They move to _DDtBase, the base class gains forcing_star = None (the
second-order forcing history that only the nodal flavour allocates), and
IntegrationPointSemiLagrangian allocates the coefficients like the other
two. The Maxwell shear-box test over all three flavours now takes the
integrator as a parameter; ETD-1 is exact for a constant strain rate, so
its tolerance is 1e-4.

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
…arries (#740, #741)

A trace-back stress history initialises its first level from the
constitutive flux of the velocity it finds at the first carry. In the
Stokes family that carry ran before the constitutive model's coefficient
update, so with integrator='etd' the recorded flux used the viscous-limit
coefficients (alpha = phi = 0): on the viscoelastic cylinder, which starts
from a velocity field already in place, that gave 6.5x the viscous drag of
BDF-1 from the first step, with pressure drag and stress norms identical
(#740). _stress_history_prepare now refreshes the coefficients as soon as
the elastic step is set; the post-carry refresh stays for the BDF ones.

The trace-back Navier-Stokes solver never set dt_elastic or refreshed the
coefficients at all, so a viscoelastic model there ran without its memory
term and ETD in its viscous limit (#741). Its solve now does both around
the carry, as the Stokes family does.

Tests: the Maxwell shear box on the trace-back Navier-Stokes solver with
both integrators, and a preset shear velocity giving BDF-1 and ETD-1 the
same first stress to O(dt/t_r). The trace-back solver's level 0 holds the
previous step after a solve (#742), so that test reads the constitutive
flux there.

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
… integration-point flavour

The second-order exponential integrator needs the forcing (strain rate)
of the previous step alongside the stress. The nodal flavour holds it as
a nodal field re-evaluated in place; the integration-point flavour had no
slot and the factory refused it. It now carries one: a point variable and
a continuous snapshot, committed from the solved strain rate by the
constitutive model's post-solve hook (an L2 fit read back at the points,
as the stress commit does) and sampled at the same departure point as the
stress before the next solve. So the forcing the parcel saw travels with
it, which is what the exponential update along the characteristic wants.

Why: on the cylinder the two-level BDF-2 stress extrapolation blew up on
the nodal history at Wi 0.5 and raised the integration-point Courant
floor at Wi 5, while ETD-2 on the nodal history stayed clean. ETD-2 keeps
one stress level and extrapolates only the forcing; this lets the
integration-point history use it.

Maxwell shear box: ETD-2 within 1% of the analytic curve on the nodal and
integration-point flavours, the two agreeing to 1e-6 (the grid flavour has
no forcing slot yet and is skipped for that case).

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
…parture points (#745)

On a box channel the integration-point stress history grew a mode in the
inlet cell column from t ~ 0.8 at Wi 1 and Wi 5, at either time step, with
DEVSS on or off and with either outlet; the nodal and grid histories were
stable, and the same channel from a gmsh rectangle was stable for all
three. The difference was the box mesh's coordinate restore: a departure
point that leaves through the inlet is clamped onto the edge and sampled
there, and the L2 commit over the inlet column fed that sample back.

The flavour now applies inflow_value: for points the trace restored (the
unclamped end point differs from the clamped one) the history takes the
inflow expression at the restored position instead of the edge sample.
Meshes without a restore function, the gmsh ones, are unchanged. With the
value set to the relaxed stress of the incoming flow the box channel is
clean at Wi 5 and Wi 1, stress peak 0.020 against 0.021 on the grid.

Test: uniform flow into a box carrying zero stress; after k steps the
inflow tensor has entered k * speed * dt and no further, and without the
value the edge sample brings in nothing.

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 Maxwell element carried its stress as a passive tensor: no
upper-convected or co-rotational term, so it was a transported linear
Maxwell fluid, and there was no solvent viscosity, so not Oldroyd-B. The
quantitative benchmarks (confined cylinder drag, 4:1 contraction,
cross-slot) all rest on the upper-convected derivative.

ViscoElasticPlasticFlowModel(..., objective_rate="none" | "upper_convected"
| "jaumann") adds L.sigma + sigma.L^T (or W.sigma - sigma.W) on the carried
stress with the current velocity gradient, linear in the unknown and first
order in time, weighted by the newest level's coefficient in the BDF form
and by alpha dt in the exponential form. Parameters.solvent_viscosity adds
2 eta_s E to the momentum flux; the history now commits history_flux, the
polymer stress alone, so the solvent part is not fed back through the
memory (the Stokes family and the trace-back Navier-Stokes solver both
read it).

Tests: start-up shear of the UCM fluid builds N1 = 2 eta lambda gdot^2
(1 - e^{-t/lambda}(1 + t/lambda)) to within 10% at dt/lambda = 0.1 on the
nodal and integration-point histories, with the shear stress unchanged
and the passive element making no N1; a solvent of eta_s adds exactly
eta_s gdot to the total shear stress.

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
DEVSS adds 2 eta_a (E - D) to the momentum flux, with D a lagged continuous
projection of the strain rate. The pair is meant to cancel, leaving only the
part of E the projected space cannot represent -- the cell-scale roughness the
term exists to damp -- so the rheology is unchanged.

The refresh that lets D track E lived only in the stress-history post-solve. A
plain Stokes solve takes the other branch, so it never ran: D stayed at its
initial zero and the term was a bare 2 eta_a E. The run then silently used
eta + eta_a, and repeated solves did not heal it, because the branch that
refreshes is not re-entered.

Measured on a fully developed channel, eta = 1, eta_a = 0.25, where the exact
gradient is -12:

  before: dp/dx = -15.00004 on every one of four successive solves (ratio 1.2500)
  after:  dp/dx = -12.00000

Refreshing after the solve alone is not enough: D starts at zero, so the FIRST
solve is still wrong and only a second call would be right -- and a steady
problem is solved once. So the fix lag-iterates within the call, refreshing D
from the velocity just found and re-solving until the pair has settled, bounded
at four passes. A time-stepping viscoelastic run already gets its catch-up
across steps and exits after the check solve.

The viscoelastic path was never affected, which is why this hid: each timestep
re-enters the stress-history branch and refreshes. The existing DEVSS test
covers only that path; the new one pins the plain-Stokes path against the exact
fully developed gradient, so it fails the moment the cancellation stops.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
A mesh read from a file (every gmsh geometry) has no analytic closure for
return_coords_to_bounds, so the property returned None and a trace-back foot
that left through an inlet was never restored: it fell through to the
evaluator's distance-weighted fallback, slow and wrong. The general facet
restore already handles the case and leaves interior points untouched, so it
is the fallback whenever there is no analytic closure, not only after a
deformation. Tested on a box written to a file and read back.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…h the conformation (#737, #768)

The integration-point history stores the stress after every solve by a global
L2 projection onto the continuous space and samples it back at the quadrature
points. That cycle has no dissipation at the cell scale, so below Courant one
a cell-scale mode of the carried stress grows from round-off (2.4 per unit time
on the Maxwell Waters-King start-up) and rings. Removing the element null
space does nothing, private nodes per cell diverge, DEVSS does not touch it;
a Laplacian term in the store projection holds it. `store_smoothing` on
IntegrationPointSemiLagrangian sets alpha = c * cell_size^2, a field so the
dose follows the local cell; c = 0.07 holds the 1/32 case unconditionally at
0.5% on the peak and the irregular mesh needs that value.

On the confined cylinder every history lost the positive-definiteness of the
conformation tau*/G + I in the first step at Courant one on the far-field
mesh: the wall shear rate is ten times the far-field one and the explicit
stretching term cannot span a step of dt*gammadot ~ 3. The viscoelastic model
gains `max_elastic_timestep(safety)`, the safety factor over the largest
strain rate, and `conformation_min_eigenvalue()`, the health line that tells
a lost preconditioner from a lost problem. Both trace-back histories expose
`carried_tensors()` for that. Hard baselines on the shear box: one BDF-1 step
from rest at dt = lambda gives eigenvalues 1 -+ 1/2 exactly; a wall four
times faster gives -1 everywhere. The ringing test at 1/16, dt 0.0125 asserts
the measured growth and its suppression; it reads the content after the
trace-back, since after the store it is zero by construction.

Documented in docs/developer/subsystems/stress-transport.md.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…transport = "forward"

A fourth stress history. The launch set is the integration points with their
weights; each step every point moves forward one step along the velocity, the
arrivals in each cell are fitted by weighted least squares to a linear
polynomial (centred and scaled by the cell, so the eigenvalue ratio of the
moment matrix is a conditioning test), and that discontinuous P1 field is
what the weak form reads. After the solve the flux is read back at the launch
points through a continuous P1 projection. A boundary cell whose face has
fluid entering and that received less than it launched has the missing share
filled with the inflow value. Nothing persists at the arrivals: one fixed
point set, one velocity-dependent map, one fit, no particle state.

Why: no launch point sits on a no-slip wall with zero velocity (the nodal
history's wall-layer defect), the per-cell fit is cheaper than sampling a
global store at every foot (17 s a step against 28 on the confined cylinder),
and the wall stress it carries is the more self-consistent (drag by stress
integral and by reaction within 1.4% where the integration-point history has
them 12% apart). Cylinder Wi 0.4, dt 0.04: drag 118.76 (reference 120.6),
conformation positive throughout.

What it is not: a scheme without the cell-scale mode. It transports its memory
consistently and rings below Courant one on a Maxwell element exactly as the
integration-point history does (Waters-King 1/16, dt 0.0125: diverges at t 2.4
unsmoothed). `flux_smoothing`, a number or a field, is the Laplacian
coefficient of the read-back projection; at c = 0.023 in cell-size units the
same case runs clean to t 8 (0.9543 at t 1, 0.5185 at t 6.5, nodal 0.9622 /
0.5171), now a level-3 test with those baselines. A swarm prototype of the
scheme that held without smoothing turned out not to be transporting its
memory at all (the proxy was rebuilt from the particles at their launch
positions before the flux was read); that is recorded so it is not rediscovered.

Serial and first order only; a periodic seam is not crossed; the launch set
does not follow a moving mesh. The setter refuses the flavour in parallel.
Documented in docs/developer/subsystems/stress-transport.md with the
recommended configuration.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
An arrival that leaves its rank's partition travels with its value and
weight to every rank, and each rank keeps the offered points that land in
its own cells and locates them exactly; only the seam layer travels, not the
set. Domain boundary faces come from the mesh's boundary label rather than
the local support count, which in parallel would have made every partition
face an inflow candidate; the inflow value is read at every boundary-cell
dof on every rank before masking, since that read is collective; a rank with
no cells keeps nothing. The serial refusals in the class and the solver
setter are gone.

Tested: np 2 and np 4 give the serial BDF-1 shear-box value 0.85136 with
arrivals crossing the seams (asserted), and the Waters-King start-up at
np 2 matches the serial trace to every printed digit (0.954273 at t 1,
peak 0.97151 at t 1.088).

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
… and restored

Neither registered with the model as a state bearer, and the forward flavour
kept its launch values in a plain array the snapshot never saw, so a restart
of either would have re-carried the wrong stress. Both now expose a snapshot
state (the slots and snapshots by name, the step bookkeeping, the forward
flavour's launch variable and smoothing) and register at construction; the
launch values live in an integration-point variable, captured with every
other variable.

Test: six steps of the Maxwell shear box, save, six more, restore, the same
six again, equal to round-off, for the nodal, integration-point and forward
histories, and the result is the twelve-step value, not a re-initialised
six-step one. Passes in serial and at np 2.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…nd run the restart test in CI

flux_smoothing is a user parameter, not evolution state: capturing it made a
load_state() revert any smoothing tuned after the snapshot, and the documented
field form (a coefficient times the mesh cell size) does not survive the
deepcopy on the in-memory path nor the disk path, which skips it. The
integration-point flavour never captured its store_smoothing; the two now
agree. The launch values are written in one assignment, so a commit costs one
PETSc flush rather than one per component. The stale docstring saying the
integration-point flavour has no checkpoint state is corrected.

scripts/test.sh pulls test_1063 forward out of the unbatched test_106* group,
as test_1072 is: it is the only guard on a stress history surviving a
snapshot and restore (level 1, tier a, 16 s). test_1060 (15 min) and
test_1061 stay out of the CI batch.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
Three reviewers over the histories (ddt.py), the constitutive and solver code,
and the tests and docs. The findings that mattered were parallel correctness and
two physics errors.

Parallel:
- The integration-point inflow write called a collective (evaluate) under a
  per-rank `if left.any()`: a run with the inlet on one rank deadlocked. The
  read is made on every rank now, and only the restored rows are kept.
- The forward fit called a collective on the populated ranks that a rank with no
  cells skipped, and then divided by zero. Ownership at a seam went through the
  evaluation slab, which could keep a point a hair across the seam on the wrong
  rank; it is strict containment now, the same rule serial and parallel.
- The file-mesh return-to-bounds fallback took every locally-single-cell facet as
  a domain boundary, so in parallel it snapped a foot crossing a partition seam to
  the seam, and read get_min_radius (collective) under a rank-local branch. It
  reads the mesh's boundary label now and the collective is unconditional.
- The forward boundary "normal" was the cell-centroid-to-face vector, not the face
  normal, so half the wall cells of an unstructured mesh were marked inflow; a
  free-slip wall flipped sign from round-off. The face normal is used, oriented
  outward, with a round-off band on the sign.

Physics:
- The BDF-2 objective-rate source was weighted by the level-0 BDF coefficient,
  which is 2 at order 2: the stretching term was doubled. The source carries
  weight one at every order.
- The theta rule read only the memory part of the stored stress, dropping the
  solvent viscosity's contribution at the old level in the composed solver.
- The DEVSS first-step refresh had been placed in SNES_TransientDarcy.solve (which
  has no DEVSS); it is in the Stokes setup path where it belongs, and the lag loop
  exits on a reduced criterion and is skipped inside a composed solver's passes.

Also: max_elastic_timestep returns a physical time when reference scales are
active, as estimate_dt does, and raises on a non-finite rate rather than
returning a silent nan; two class docstrings were shadowed by attributes placed
before them; a history committed before the first carry is kept, not overwritten;
the forward launch set refuses a moved mesh (by a coordinate stamp, since the
topology version bumps on every variable creation); the integration-point history
refuses ignored bcs; dead state removed.

Tests: the parallel forward test runs on a gmsh mesh with a non-uniform stress and
a flow that crosses a seam of either orientation, and matches serial to 1e-4 (a
per-cell fit sums its arrivals in a partition-dependent order); the three
flavour-vs-flavour tests carry hard serial baselines; the solvent test constrains
the assembled momentum flux through a channel's first-step speed; the ringing test
pins the growth rate rather than a round-off-seeded amplitude; markers corrected;
the file-mesh test renumbered off its clash with test_0073. Docs numbers
reconciled with the tests. test_1063 wired into CI (16 s); test_1060 and test_1061
stay out (15 min and longer).

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…w ones tier C

development now runs tests/test_106* and test_107* as a whole band and excludes
the slow cases by being tier_c (reported, not gating), the way test_1064 is,
rather than by number. So the line that pulled test_1063 forward is redundant:
the band runs it, and the coverage checker (check_test_coverage.py, also from
development) verifies every file is reachable. test_1063 (16 s, tier_a) gates;
test_1060 (15 min) and test_1061 (20 min) move to tier_c so they report their
hard baselines without gating, as test_1064 does.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…transport = "lagrangian"

The internal-swarm Lagrangian DDt carries the stress on a swarm the solver owns
and advects, reading it back through a discontinuous cells proxy and never
projecting it to the mesh. It is now selectable the same way as the other
histories, alongside semi_lagrangian, integration_point, forward and eulerian;
the same scheme on a user-supplied swarm stays available as Lagrangian_Swarm
passed through DFDt= (test_0070).

Two defects had to be fixed for it to work through the solver factory:

- update_post_solve and initialise_history evaluated the constitutive flux and
  wrote each tensor component in the same loop. The flux reads psi_star itself,
  and writing one component marks the swarm proxy stale, so the next component
  read a half-updated history and the relaxation advanced twice a step (recorded
  values ran sigma_1, sigma_3, sigma_5...). Every component is evaluated before
  any is written now, the pattern Lagrangian_Swarm already used (audit SWARM-06).
  The flavour then matches the integration-point history to 2e-8 on the uniform
  Maxwell shear box.
- the owned swarm emptied downstream on an open boundary and the cells proxy had
  nothing to interpolate. It now sets population control that both fills starved
  cells and caps over-full ones, with bounds from the initial occupancy so the
  swarm neither starves nor grows without bound.

Order 1 BDF only for now: the particle flavour has no exponential coefficients
and order 2 is not yet validated, so the factory refuses the exponential
integrator and order 2 cleanly rather than crashing in the first solve. The
conformation check refuses it, as it does any history without a per-point tensor.

Tests: lagrangian joins the test_1059 Maxwell-box matrix (held to the same 1e-6
spread as the mesh flavours) and the test_1063 restart round-trip (the raw
per-particle array is not compared, since population control need not reproduce
the particle set). Documented as the fifth history.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
The particle counterpart of AdvDiffusionSLCN (semi-Lagrangian) and AdvDiffusion
(streamline-upwind). It carries the advected scalar's history on a swarm the
caller supplies through a Lagrangian_Swarm manager, so it rides on the material
swarm a coupled model already has rather than a private one. The solver declares
its history variable on the swarm (so build it before populating) and does not
advect the swarm itself: the flow solver advects the shared swarm in a coupled
run, and a standalone loop advects it before each solve. It warns once if the
swarm has not moved between solves, so a forgotten advection is not a silent,
un-transported answer.

Validated against the diffusing rotating Gaussian (uw.analytic.RotatingGaussian,
exact at every time) on a disc: over a half revolution it gives the integral-L2
error 1.7e-2, matching SLCN (6.4e-2) and SUPG (6.5e-2) and diffusing to the exact
peak, held to a hard baseline in tests/test_1101. The particle-scheme trap, a
half-blend (step_averaging=2) that keeps the particle's old sharp value and
under-diffuses, is caught by a second test.

Two things the demonstrator turned up and the docstring now records. First, the
mesh solution must return to the particles in full each step (PIC with
step_averaging=1) or the field under-diffuses; a diffusion-matched FLIP retention
keeps sub-cell sharpness where that matters. Second, a particle scheme needs its
cells kept populated: a square box under rigid rotation loses its corners (they
sit at radius > side/2 and rotate out of the domain), which starves the corner
cells and destabilises any particle scheme though the mesh schemes are
indifferent, so rotation demonstrators use a disc or annulus.

scripts/test.sh widens the adv-diff batch glob to test_110* so the new test runs.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
lmoresi and others added 2 commits September 23, 2026 10:40
…rossing flow

The particle advection-diffusion solver diverged over a full rotation on a square
box (peak 13.6 by a half turn) while holding a disc. The cause was not
depopulation: tracking the per-cell particle count, the particle values and the
nodal projection showed the cells never starve (min 15 particles each). Where the
flow crosses a wall, out-flowing particles are clamped back just inside it and
pile up in a thin layer along the boundary (249 in one cell). The history proxy
was degree 2, and a six-coefficient per-cell least-squares fit of particles strung
along a line is only mildly ill-conditioned, under the projector's guard, so it
overshoots the nodal value about fivefold; that overshoot enters the solve,
corrupts the particles a few steps later, and diverges.

A degree-1 fit needs three spanning points, is well conditioned on the clamped
layer, and holds. It defaults to degree 1 now (`proxy_degree`, exposed for the
rare case where the flow keeps every cell's particles well spread and the extra
accuracy is wanted). Measured on the rotating diffusing Gaussian: the square box
now runs to a full revolution bounded and diffusing to the exact peak (half-turn
error 6.7e-2 where degree 2 blew up), and the disc is 2.0e-2 over a half
revolution, still between SUPG (1.7e-2) and SLCN (2.5e-2). test_1101 gains the
box case; the disc baseline moves to the degree-1 value.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
This branch is the base of a four-PR stack (#785 -> #789 -> #795 -> #800) and
was 71 commits behind development, which is why all four showed the same two
CI failures: tests development had already renamed
(test_1060_nitsche_freeslip's absolute 1e-4 bound and
test_1070_free_surface_plume's strong-vs-penalty ratio). The refresh belongs
here, at the root, so each PR's diff against its base stays its own work rather
than growing development's history.

Two conflicts, both this branch's own additions. scripts/test.sh keeps the
broader `tests/test_110*py` glob -- development narrowed it to `test_1100*py`,
which matches none of the other test_110* files here (#721 recurring);
scripts/check_test_coverage.py verifies it. Four of the five ddt.py hunks have
an empty development side, and the fifth takes both: this branch's base-class
`update_exp_coefficients` / `_exp_alpha` / `_exp_phi` (#739, which replaced
three per-flavour copies dating from b6b0e7a in April) and development's
`_note_history_shift`. Verified on the merged tree: one definition of each.

Underworld development team with AI support from Claude Code
Copilot AI balanced review requested due to automatic review settings October 7, 2026 23:13

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

Parallel boundary restoration, duplicated inflow terms, unreliable movement detection, and several contract violations remain unresolved.

20 open findings
What changed in this PR

This PR generalizes viscoelastic stress-history transport and adds swarm-based advection-diffusion support.

Changes:

  • Adds Eulerian, forward, integration-point, and particle stress-history capabilities.
  • Adds objective-rate, DEVSS, inflow, restart, and conformation diagnostics.
  • Adds extensive serial/MPI tests and stress-transport documentation.
File Description
src/​underworld3/​systems/​ddt.py Expands history managers and transport contracts.
src/​underworld3/​systems/​solvers.py Integrates stress histories, DEVSS, and swarm advection-diffusion.
src/​underworld3/​systems/​navier_stokes_eulerian.py Connects transported stress to Eulerian Navier–Stokes.
src/​underworld3/​systems/​__init__.py Exports new solver and history APIs.
src/​underworld3/​constitutive_models.py Adds objective rates, solvent stress, and health diagnostics.
src/​underworld3/​discretisation/​discretisation_mesh.py Restores out-of-bounds points on file meshes.
docs/​developer/​subsystems/​stress-transport.md Documents transport choices and limitations.
docs/​developer/​index.md Registers the new subsystem guide.
scripts/​test.sh Includes the expanded advection-diffusion test band.
tests/​test_0074_return_to_bounds_on_a_file_mesh.py Tests file-mesh boundary restoration.
tests/​test_1059_stress_transport.py Covers stress transports and constitutive behavior.
tests/​test_1060_stress_store_smoothing.py Tests integration-point smoothing.
tests/​test_1061_stress_forward_history.py Tests stabilized forward transport.
tests/​parallel/​test_1062_forward_stress_history_mpi.py Checks forward transport under MPI.
tests/​test_1063_stress_history_restart.py Tests history snapshot restoration.
tests/​test_1101_advdiff_swarm_rotating_gaussian.py Validates swarm advection-diffusion.

🧠 Review effort: Balanced


Give feedback about Copilot approvals in this survey to enter a drawing for a $150 gift card.

Comment on lines +2323 to +2324
E = sympy.Matrix(self.Unknowns.E)
rate = sympy.sqrt(2 * (E.T * E).trace())
Comment on lines +3615 to +3616
return (numpy.array(facets, dtype=int).reshape(-1, d),
numpy.array(opp, dtype=int).reshape(-1))
Comment on lines +2349 to +2350
for boundary in self.mesh.boundaries:
solver.add_natural_bc(condition, boundary.name)
Comment thread src/underworld3/systems/ddt.py Outdated
Comment on lines +5554 to +5555
coords = np.asarray(self.mesh.data)
return (coords.shape, float(coords.sum()))
Comment on lines +5041 to +5043
def _swarm_position_key(self):
c = np.asarray(self.swarm._particle_coordinates.data)
return (c.shape, float(c.sum()), float((c * c).sum()))

import itertools

pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] # several hundred solves: minutes

import underworld3 as uw

pytestmark = [pytest.mark.level_2, pytest.mark.tier_c] # 15 min: a slow characterisation with a hard baseline, reported not gated (as test_1064 is)
Comment on lines +15 to +18
pytestmark = [pytest.mark.level_3, pytest.mark.tier_c] # 20 min on Waters-King: reported not gated, as test_1064 is


def test_the_forward_history_with_its_read_back_smoothing_holds_the_maxwell_start_up():

import underworld3 as uw

pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] # solves: not level 1

import underworld3 as uw

pytestmark = [pytest.mark.level_2, pytest.mark.tier_b]
The deprecated-pattern scan fails this branch on one line: `ddt.py:5554`
reached the vertex coordinates through `self.mesh.data`, which S7 reserves for
old code. `mesh.data` itself emits a DeprecationWarning naming the replacement.

Checked rather than swapped blind, because this stamp decides whether a
snapshot restore is refused: on a unit box at cellSize 0.25 both return shape
(30, 2) with coordinate sum 30.113992809675707, identical. So the stamp means
exactly what its docstring says it means -- the vertex count and the
coordinate sum -- before and after.

tests/test_1063_stress_history_restart.py and
tests/test_0074_return_to_bounds_on_a_file_mesh.py: 5 passed.

Underworld development team with AI support from Claude Code
@lmoresi
lmoresi merged commit f4733d4 into development Oct 8, 2026
2 checks passed
lmoresi added a commit that referenced this pull request Oct 8, 2026
The rename lands in four places the branch had not reached, and one anchor was
stale.

`SemiLagrangian` is a factory over (trace, launch), not a class, so it carried
no class-level description -- and `tests/test_0017_describe_and_render.py`
requires every family to answer at that level, while
`docs/developer/subsystems/describe-and-view.md` calls
`SemiLagrangian.describe_class()` directly. The factory now delegates
`describe_class` and `view` to the default scheme, which is what it builds when
neither axis is given, and exposes `schemes` for a caller that wants another.

Three expectations follow the rename: the describe contract and the step
transcript now name `BackwardNodesSemiLagrangian`, and the transport-schemes
guide's front matter lists the four trace-and-launch families instead of two
names no class carries (`tests/test_0030_capability_guides.py`).

`docs/advanced/eulerian-advection-diffusion.md` described three managers; the
section now says how the two axes select a scheme, and records that a forward
trace's per-cell fit is unstable where the flow empties a cell (#811).

The eulerian cavity anchor is re-recorded. It moved 4e-05 relative on merging
development, and the new value is development's own: the same cavity with the
default EulerianSUPG history gives [-0.129019, -0.05467728] there, before and
after #833, unchanged to nine digits from tolerance 1e-8 through 1e-12. The
original was recorded 2026-09-29, before this branch last merged development.

That leaves something worth knowing: the Eulerian SUPG velocity answer moved
on development and nothing there anchors it. This file is the only test that
can see it.

level_1 tier_a: 1428 passed, and the one failure it had (the capability guide)
is fixed here. test_1103: 7 passed.

Underworld development team with AI support from Claude Code
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants