Repository navigation
Integration-point stress history grows an inlet instability when dt is halved (Wi 5 cylinder) #737
Description
Activity
Discriminators run (1/20, dt = 1/240, DEVSS eta_a = nu/4, 240 steps to t = 1). Stress peak by step:
run 100 140 180 220 240 IP, Wi 5 0.087 0.146 0.493 1.602 2.281 nodal, Wi 5 0.178 0.268 0.330 0.365 0.371 IP, Wi 1 0.175 0.242 0.252 0.254 0.254 grid, Wi 5 0.088 0.099 0.112 0.116 0.114 The nodal trace-back at the same dt and Wi does not grow, so the boundary-carry candidate in the report above is not the mechanism: the clamp at the inlet is shared by both trace-back flavours and only the integration-point one grows. The IP history at Wi 1 does not grow either. The growth needs all three of: the integration-point sampling, the halved time step (Courant about 0.3 at the inlet centre instead of 0.6), and the long relaxation time (per-step retention 0.992 at Wi 5 against 0.958 at Wi 1).
That is the low-Courant growing mode of the integration-point scheme (#703 analysis: inter-cell aliasing of the shifted sampled field, a creep of ~1e-3 per step on the scalar problem) with the viscoelastic memory term as its amplifier: whatever the transport gains at grid scale is re-injected almost in full at high Wi, which turns a creep into a doubling every ~25 steps. It shows at the inlet first, where the flow is unidirectional and the stress signal is smallest, so the mode is not masked.
The damping the #703 analysis calls for, kappa_stab = eps h^2 / (pi^2 dt), has no diffusion term to attach to in a stress history; the analogue here is a smoothing weight on the per-step commit projection (
SNES_MultiComponent_Projection, currently smoothing 0). A run with that weight set from the formula is the confirming test, and if it holds the option belongs with the stabilisation contributions being collected on the planning list, off by default.Ruling (Louis): the integration-point history is operated at Courant ~1, as settled for the scalar case, and the damping term is not the remedy — it was not liked then and any use of it needs a careful look first. Below Courant 1 the scheme is to be treated with care, more so with a memory term. For convergence checks on this history, refine h and dt together at fixed Courant rather than halving dt alone. The grid history is unaffected by dt and remains the reference for a dt sweep.
Integrator check. ETD-1 on the integration-point history reproduces BDF-1 step for step: flat at dt 1/120 for 240 steps at Wi 2, 5 and 10, and at Wi 5, dt 1/240 the same growth with the same doubling time (stress peak 0.14 to 2.10 over steps 140 to 240, against 0.15 to 2.28 for BDF-1). Consistent with the mechanism as stated: the amplifier is the per-step retention of the carried stress, and e^{-dt/lambda} and lambda/(lambda+dt) agree to four digits there. The integrator is an accuracy lever on this history, not a stability one.
Order 2 on the integration-point history, revisited. BDF-2 on IP at Wi 5, dt 1/120 doubles every ~10 steps from step 60 (peak 5.06 by step 120) where BDF-1 and ETD-1 are flat. ETD-2 with a traced forcing history (5e1266e: the strain rate sampled at the same departure point as the stress) is flat for 240 steps at the same settings, peak 0.09 throughout. So the second level itself is not the problem; extrapolating the stress across two levels with opposite-sign weights is. The Courant floor of this history is unchanged by the integrator (ETD-1 = BDF-1 step for step at dt 1/240), and order 2 is now available on it without lowering that floor.
Cross-reference #745: on a plain box channel with no cylinder the integration-point history grows at the inlet from t ~ 0.8 at Wi 1 and Wi 5, at dt 1/60 and 1/120 alike and with DEVSS off, while nodal and grid are stable. That inlet growth is dt-independent there, so it is a different mechanism from the low-Courant mode discussed here, though the dt 1/240 cylinder case above also grew at the inlet first.
A benchmark with an exact reference (Waters & King 1970 start-up, El = 1, beta = 0) sharpens this from "do not halve dt alone" to a monotone degradation with a mechanism.
All at 1/32, the same mesh, only the step changing:
dt Courant peak settles to ripple (t 7-10) 0.04 1.28 0.95757 0.50493 0.0203 0.02 0.64 0.96949 0.49038 0.0165 0.01 0.32 0.97830 0.37367 0.2786 0.005 0.16 diverged at t 3.02 reference 0.99880 0.50550 0.0230 Clean, clean, bounded ringing, divergence. The mesh is not involved: 1/16 and 1/32 at dt = 0.02 agree to four digits (0.96955 / 0.96949).
Integrator-independent, as this issue says. At 1/32, dt = 0.01, ETD-1 rings identically to BDF-1 (peak 0.97598 vs 0.97830, settles to 0.37644 vs 0.37367, ripple 0.277 vs 0.279). At dt/lambda = 0.01 the two retention factors differ by 5e-5, so there is nothing for the integrator to change.
Proposed mechanism. The integration-point history stores at quadrature points and must rebuild its field from a per-cell fit to be sampled at a departure point. That fit error is spatial, so it is injected once per step and is roughly independent of dt, while relaxation only damps what is already present, by exp(-dt/lambda) per step. For dt << lambda the noise floor is then ~ eps*lambda/dt: halving the step doubles it. Beyond some point injection outruns damping and bounded ringing becomes growth, which is the 0.01 -> 0.005 transition above.
The control that makes the case. The nodal trace-back on the same problem, same mesh, same steps:
dt peak settles to ripple 0.01 0.97835 0.50537 0.0221 0.005 0.98468 0.50543 0.0223 It converges first order in dt toward a spatial floor near 0.8% of the peak, with the settled value and the physical ripple unmoved. Waters & King is unidirectional, so the stress is x-independent while the trace-back displacement is purely in x: exact transport is the identity, and a nodal interpolation reproduces a field constant along the displacement direction exactly — zero injection. The integration-point reconstruction does not, which is the whole difference.
If that is right, the fix is on the reconstruction (a representation that is exact for the fields it is asked to carry, or an update that does not round-trip through a fit every step), not on the integrator and not on the timestep heuristic.
Diagnosis update (2026-09-20), Waters & King 1/32, dt 0.01, and the confined cylinder.
- Both trace-back flavours store the stress after every solve by the same global L2 projection onto continuous P1; the IP flavour only samples that field at the quadrature points instead of the vertices. Nothing is carried at the points between steps.
- The per-cell element null space of the point values is real and grows with the blow-up, but removing it every step changes nothing to six digits over 500 steps (regular and irregular meshes). It is a symptom, not the driver.
- Private nodes per cell (
continuous=False, degree 1 and 2) diverge by t 2 at 1/16, dt 0.0125. Shared vertices are what damp the inter-cell jump mode. - DEVSS at eta_a = eta_p does not affect the ringing.
- A Laplacian term in the store projection does: alpha ≈ gamma·dt·(h/pi)^2, gamma the measured growth rate (2.4 per unit time on the Maxwell WK case, 1.8 at beta = 0.2). Measured window at 1/32, dt 0.01: 3e-6 rings, 1e-5 stable over the run (0.1% on the peak), 3e-5 unconditionally (0.5%), 1e-4 too strong (3%). As a field, alpha = c·h(x)^2 with c = 0.07 on the irregular mesh: peak 0.9747 / settled 0.4882 against nodal 0.9784 / 0.4864.
- On the cylinder (beta 0.59) the mode does not appear over the run lengths used; the failures there were the timestep against the wall strain rate (conformation lost in the first step at dt 0.4; admissible everywhere at dt 0.04 to Wi 0.6).
Not yet in the code: the smoothing is applied by wrapping the commit projection's
solvein the run scripts (~/+Simulations/stress_transport/waters_king/wk_smooth_field.py,cylinder_ob/ob_cylinder.py -uw_store_smooth c). The natural home is an option onIntegrationPointSemiLagrangiansetting_commit_projection.smoothing = c * mesh.cell_size()**2.Underworld development team with AI support from Claude Code
- added a commit that references this issue
on Sep 20, 2026 Implemented on feature/stress-transport-ddt (71bb82d):
IntegrationPointSemiLagrangian.store_smoothing(alpha = c·h², c ≈ 0.03–0.07),ViscoElasticPlasticFlowModel.max_elastic_timestep(safety)andconformation_min_eigenvalue(), hard-baseline tests, and docs/developer/subsystems/stress-transport.md. Closes with the branch PR.Underworld development team with AI support from Claude Code
Fixed by #833, merged to
development(f4733d4f).store_smoothingonIntegrationPointSemiLagrangianputs 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.The mechanism, for the record: the history stores the stress by a global L2 projection onto the continuous space and samples it back at the quadrature points, and that cycle has no dissipation at the cell scale — so 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. Three things that did not work, so nobody retries them: removing the element null space (inert), private nodes per cell (diverge), DEVSS (does not touch it).
Guarded by the ringing test at 1/16,
dt 0.0125, which asserts the measured growth and its suppression — and reads the content after the trace-back, because after the store it is zero by construction.Closing manually: this repo's
Closeslines fire whendevelopmentreachesmain, and merges tomainare infrequent.
On the viscoelastic DFG cylinder (pure Maxwell, Wi 5, 1/20 mesh, DEVSS eta_a = nu/4, EulerianSUPG momentum) the
integration_pointstress history is stable at dt = 1/120 for 240 steps and grows an instability at dt = 1/240 that starts at the inlet, not at the cylinder. Theeulerian(grid) history at the same two time steps is unchanged to 0.7% in drag.Wi 5, 1/20, drag window t = 0.67 to 1, near-cylinder stress L2, max stress peak:
Stress peak on the IP dt = 1/240 run by step (dt = 1/240, so step 240 is t = 1): 140: 0.15, 160: 0.27, 180: 0.49, 200: 1.04, 220: 1.60, 240: 2.28. Doubling every ~25 steps, still growing at the end of the run. The stress invariant frames show it as a faceted, cell-scale stress structure attached to the inlet boundary and spreading downstream to the cylinder; the grid run is clean there.
Candidate mechanism, not yet tested. The trace-back flavours do not use
inflow_value; a departure point that lands outside the domain is restored to the boundary and takes the transported field's value there (docstring oninflow_value). At the inlet that is an exact carry of the boundary value from the committed flux, with no relaxation applied on the way, so a boundary node whose stress is above the relaxed inflow stress feeds itself through the memory term. The same exact-carry-at-a-boundary mechanism is the candidate for #735 at the no-slip wall. The grid history is clean because it has a real inflow condition (-uw_stress_inflow 1in the benchmark script, the relaxed stress of the incoming flow).Predictions being tested: the nodal trace-back at the same dt and Wi should show the same inlet growth (same boundary carry); the IP history at Wi 1 and dt = 1/240 should not, or much more slowly (relaxation per step is five times stronger).
Runs:
~/+Simulations/stress_transport/cylinder/{w5_grid,w5dt_grid,w5_ip,w5dt_ip}, scriptrun_wi5_dt.sh; framesframes_w5dt_tauII(steps 120, 168, 232). Worktree branchfeature/levelset-supgat b5bfc29.