Eulerian stress history blows up at a shear-wave reflection (Waters & King start-up)

[lmoresi]

What

The Eulerian (grid, streamline-upwind) stress history fails on the Waters & King 1970 start-up problem, where the two trace-back histories do not.

Channel between plates at y = ±h, fluid at rest, constant pressure gradient switched on at t = 0. Maxwell (β = 0), elasticity number El = λη/ρh² = 1, so the problem is hyperbolic and the shear-wave fronts cross the channel at t = 1, 3, 5.

stress_transport = "eulerian", BDF-1, resolution 1/16, Δt = 0.02. The run tracks the first overshoot correctly:

t	u centreline	reference
1.00	0.9526	0.9988 at t = 1.005
1.50	0.8774
2.00	0.7266

then oscillates with growing amplitude from t = 2.3:

2.300000,-0.84525938
2.320000, 1.04146532
2.340000, 2.55333802

and the linear solve stops converging (DIVERGED_LINEAR_SOLVE, then the SNES makes no progress and the run hangs). Reproduced twice at the same settings.

The failure is at the first reflection of the shear-wave front off the far wall. semi_lagrangian and integration_point run the same problem to t = 10 without this.

Where

SNES_Stokes.stress_transport == "eulerian" builds the history on the Eulerian SUPG DDt manager. The suspicion is the upwinding: at a reflection the stress discontinuity travels against the local flow direction rather than along it, and the streamline-upwind term stabilises the wrong direction.

Reproduce

~/+Simulations/stress_transport/waters_king/wk_uw.py (session scripts, not in the repo):

python wk_uw.py -uw_stress_transport eulerian -uw_resolution 16 -uw_dt 0.02 -uw_out out.csv

Records in ~/+Simulations/stress_transport/waters_king/wk_eulerian_r16_dt0.02.blowup.csv; the benchmark ledger has the full entry.

Why it matters

The grid history is the partition-independent one and the cheapest, and it is the flavour the cylinder drag tables have been run with. This says it cannot be trusted where the flow is elastically hyperbolic and a front reflects — which is exactly the regime the high-Wi runs are heading into.

[lmoresi]

The timestep sweep is in, and it removes the mechanism I proposed in the issue body.

At 1/16 the grid history diverges at every step tried:

Δt	Courant	diverges at
0.01 (1/32 mesh)	0.3	t = 1.41
0.02	0.3	t = 2.34
0.04	0.6	t = 2.56
0.08	1.3	t = 3.20

A larger step postpones the failure and does not prevent it, so this is not a Courant threshold.

I had proposed that the SUPG weight explains it: tau = weight / sqrt(transient + advective + viscous) with transient = (c_t*c0/dt)^2 and advective = (c_u*|a|/h)^2 (ddt.py, EulerianSUPG.tau), so as dt falls the transient term scales the upwinding away, with the crossover at exactly Courant one. That predicted dt = 0.08 would be stable. It is not. Withdrawn — the issue stands on the observation.

For contrast, on the same problem at 1/32 and Δt 0.01 the nodal trace-back gives peak 0.9783 at t = 1.060 and settles to 0.5054 against the reference 0.5055, with the residual ripple right to two digits; the integration-point history reaches the same peak and then rings permanently (#737). So it is the grid flavour specifically, not the timestep and not the resolution.