Viscoelastic flow in Underworld3: the methods library, the benchmarks, and what is open (October 2026)

[lmoresi]

This thread gathers the viscoelastic (VE) flow work of September and October 2026 in one place: what the library now contains, what the benchmarks say, and what is open. The Project items under Area solvers with "VE" in the title carry the individual tasks; this is the narrative they point back to.

The library

We treat VE flow as a library of methods rather than one scheme, the way a code carries several turbulence models. Three independent axes name a constitutive model (ViscoElasticPlasticFlowModel):

• element: the spring-dashpot arrangement, maxwell or jeffreys (a parallel solvent dashpot);
• relaxation: linear (Hookean spring) or fene_p (finitely extensible, Peterlin closure, Parameters.extensibility = L²);
• objective_rate: upper-convected and the others.

UCM is maxwell + linear + upper-convected; Oldroyd-B is jeffreys + linear + upper-convected; FENE-P is jeffreys + fene_p + upper-convected. The stress is carried by one of four semi-Lagrangian histories (stress_transport: backward or forward, from the nodes or from the integration points) or by the Eulerian SUPG/DEVSS path, in a log-conformation store with an exponential (ETD) step that keeps the conformation positive definite. The forward histories read their arrivals back either by a per-cell fit (reconstruction="cell") or by a global weighted least-squares projection (reconstruction="global", ParticleL2Projector), which is the stable choice and needs no smoothing or limiter. Regularisations exist as named, visible terms: store smoothing (α = c h²), a Barth–Jespersen limiter on the cell fit (kept, off, judged the wrong bound), a hierarchical bubble penalty for P2 projections, and a deficit fill where particles leave cells empty. The composed Navier–Stokes and advection–diffusion solvers take velocity_transport, transport and stress_transport so the momentum and stress histories are chosen independently. Docs: docs/developer/subsystems/stress-transport.md; PRs #789, #795, #800.

What the benchmarks say

Benchmark	Reference	Status
Steady simple shear, Oldroyd-B and FENE-P	closed forms (FENE-P: f²(f−1) = 2Wi²/L²)	Done. Oldroyd-B exact; FENE-P first order in dt/λ, errors halve cleanly (test_1068).
Waters and King start-up shear	analytic transient	Done: the nodal history converges, the integration-point history degrades as dt falls, the Eulerian history diverges at every dt.
DFG 2D-2 Newtonian cylinder, Re 100	cd 3.22–3.24, cl 0.99–1.01, St 0.295–0.305	Passed at res 40 for lift and St; drag 2% low (resolution).
Cross-slot instability, Oldroyd-B	pitchfork at De 0.35 (Poole)	Bifurcation found with the log store, SPD update and exp read, at De 0.71. Threshold off by two; parked.
VE cylinder, Re 100	Richter et al. 2010 (FENE-P, 3-D, β 0.9); no reference for the DFG geometry	Open. Shedding stops at Wi 0.5 with the global projection; the cell fit and the Eulerian stress keep shedding. Qualitative, scheme-dependent, no number to match.

The next quantitative target is the one classic VE cylinder benchmark we have not run: creeping flow past a confined cylinder, Oldroyd-B, β 0.59, 2:1 blockage, drag against Wi tabulated to four figures by Alves, Fan and Hulsen (132.36 at Wi 0 falling to about 117.8 by Wi 0.7). It needs a 2:1 channel and the Stokes limit of the existing driver, nothing new in the library, and it is the first place the four histories and two reconstructions can be ranked against a published number.

What the cylinder runs taught us

• At Re 200 the lift collapses with resolution in every scheme; the resolved answer is weak shedding with stress sheets on the separated shear layers. The forward family reaches a given answer about four times coarser than the Eulerian path.
• The Strouhal number follows the velocity history (SUPG 0.30, forward nodes 0.25 to 0.28), not the stress history.
• At Re 100, Wi 0.5, res 40, the stress reconstruction decides the regime: global projection goes steady, cell fit (cd 4.22, cl 0.24, growing) and Eulerian stress (cd 4.33, cl 0.52) keep shedding. Two readings remain: the projection damps the shear-layer perturbation, or the cell fit's overshoot at the rear stagnation point (Forward integration-point history: the per-cell fit is unstable where the flow empties a cell #811) is a perturbation source. A res-80 pair separates them.
• FENE-P at Wi 0.5 sits in the same near-steady regime with 12% less drag than Oldroyd-B (wall spring factor 2 to 3). Wi 2 is unresolved at res 20 for both laws.
• The FENE-P stall was the record: log c projected nodally overshoots L² at the wall. The record is now log(f c), the Oldroyd-B record, decoded through its trace (commits 3170929, 020b02f). Lesson: a bounded variable is not the one to project.
• A free-slip cylinder gives a steady birefringent strand and no shedding: a different, steady benchmark, not a sharper version of this one.

Open

Res-80 pair at Re 100 Wi 0.5; the creeping confined cylinder driver; Wi 2 at res 40; runs biased by the first deficit fill to rerun or discount (nodal_forward_global3, r40_fwd_Wi1, r40_fwd_Wi2, r20_fwd_Wi2); global projection as a swarm proxy option; second order for the forward methods; retiring the old SNES_NavierStokes and SNES_AdvectionDiffusion once the composed solvers settle; the cross-slot threshold; FENE-P at Re 100 against Richter 2010.