Repository navigation
Stress transport on the grid, forward from launch points, and on particles - #833
Merged
Merged
Conversation
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
…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
Contributor
There was a problem hiding this comment.
🟡 Changes recommended
Parallel boundary restoration, duplicated inflow terms, unreliable movement detection, and several contract violations remain unresolved.
20 open findings
Base timestep cap on the active objective operator · New Handle ranks with no physical boundary facets · New Avoid double-registering conditions on boundary labels · New Use a robust mesh geometry version check · New Track swarm movement beyond coordinate moments · New Validate that safety is positive · New Reject unsupported constructor options · New Check lag convergence before repeating the Stokes solve · New Document forward transport failure on sparse empty cells · New Document lagrangian inflow and viscoelastic failure limits · New Document DEVSS shear-wave reflection instability · New Clarify DEVSS refresh behavior across Stokes branches · New Preserve exact rational transport coefficients · New Move tier_b onto the test while retaining module marks · New Ensure release-gate tests have exclusive tier markers · New Apply tier_a to each test instead of only the module · New Apply tier_c to the test instead of only the module · New Ensure each test has exactly one tier marker · New Move tier_a onto the parametrized test · New Apply tier_b to each test instead of only the module · New
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 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
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
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.



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/-278across 16 files, of which 1,159 lines are tests.What it adds
A stress history that moves on the grid.
stress_transporton 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_smoothingputs a Laplacian in the store projection withalpha = c * cell_size^2— a field, so the dose follows the local cell;c = 0.07holds 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 + Iin 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 spandt*gammadot ~ 3. The viscoelastic model gainsmax_elastic_timestep(safety)andconformation_min_eigenvalue()— the health line that tells a lost preconditioner from a lost problem — and both trace-back histories exposecarried_tensors()for it.Issues in its territory it does NOT close
psi_star[0]holds the new stress on Stokes and the previous step's onNavierStokesSLCNnavier_stokes_eulerian.py, but no commit addresses the convention difference, and nothing here reconciles itLagrangian_Swarmgives a NaN first residual on the VE cylinderstress_transport = "lagrangian"is added here, and this is its open defect on that benchmarkfixed-in-PR"forward"transport this branch introduces.datawrites against Charter §7enhancement, it is a conformance choreWorth 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, andtests/parallel/test_1062_forward_stress_history_mpi.py.Hard baselines rather than method comparisons: one BDF-1 step from rest at
dt = lambdagives conformation eigenvalues1 -+ 1/2exactly; a wall four times faster gives-1everywhere. The ringing test at 1/16,dt 0.0125asserts 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 fromdevelopmentbetween — they touch the same region ofddt.py.Underworld development team with AI support from Claude Code