Navier-Stokes with Eulerian SUPG momentum transport, and a partition-independent cell size (#687)

[lmoresi]

uw.systems.NavierStokesSUPG: the incompressible Navier-Stokes equations on the Stokes saddle-point solver with the momentum advection assembled implicitly and stabilised by the vector SUPG term, the counterpart of the scalar Eulerian solver in #673 (this PR stacks on that branch). Crank-Nicolson at order 1, BDF2 at order 2, the velocity history on the mesh; no stress history, the viscous stress at an earlier level is rebuilt from the stored velocity through the constitutive model, and pressure has no history. The advecting velocity is a choice: the second-order extrapolation 2uⁿ − uⁿ⁻¹ (one linear solve per step, the default), Picard passes on the latest iterate, or the unknown itself under Newton.

What the measurements say (design note, "Navier-Stokes with SUPG momentum transport"; study directory ~/+Simulations/navier_stokes_supg/):

• Kovasznay flow at Re 40: SUPG 6.6e-4, 2.6e-4, 6.7e-5 at h = 1/16, 1/32, 1/64; Galerkin (τ = 0) 1.1e-4, 1.6e-5; the semi-Lagrangian solver 5.8e-3, 2.9e-3 at seven times the cost per step. At this Reynolds number the stabilisation is not needed and costs accuracy; the missing viscous term in the first-derivative residual is the O(h²) consistency error. Newton, Picard and Courant 4 give the same steady state, as they must.
• Lid-driven cavity: Re 100 within 4% of Ghia on a 1/32 mesh at a fixed point; Re 400 at 1/48 within 6% (the mesh), where the extrapolated step at Courant 2 keeps a localised alternating mode of the lagged coefficient that one Picard pass, or Courant 1, removes; Re 1000 with geometric multigrid on the velocity block steps at one Krylov iteration per step but needs a longer run than was made here.
• Cylinder wake (DFG 2D-2, Re 100): Strouhal number 0.298 against the reference 0.295 to 0.305, the two fully implicit forms agreeing to three digits and the extrapolated step within 1% in frequency; the lift peak 11% low and the drag 23% low on the finer mesh, a deficit that does not close with resolution and points at the stabilisation's diffusion (the tau weights are the scalar solver's, untuned); the Galerkin control cannot run at this element Reynolds number, which is what the term is for. The semi-Lagrangian solver on the same mesh has the shedding 13% too slow at three times the cost.

Found and fixed on the way

• mesh.cell_size() was partition-dependent (mesh.cell_size() is partition-dependent: the per-cell radius comes from a kd-tree over the rank's centroids, so cells at partition boundaries get a different size on each rank count #687): the kd-tree radius picks the nearest centroid among the rank's cells, so the SUPG τ and the Nitsche penalty differed across a seam (two-rank Kovasznay error 5e-4 off serial, 1e-15 with a constant h), and after a deform it read stale vertex coordinates. The field now reports each cell's own radius from the DM's coordinates; the kd-tree radii still feed get_min_radius.
• The pressure gradient belongs in the SUPG residual: without it the term is O(1) at the exact solution and Kovasznay read three times worse.

Tests: the API contract (construction rules, one linear solve per step, Picard passes, the Stokes limit, runtime-constant timestep and theta), a two-rank Kovasznay error that matches serial to 1e-7, the scalar parallel reference re-recorded for the new cell size, and the Nitsche local-h tests reading the field's definition.

Deferred, recorded in the note: a recovered Laplacian to complete the residual for P2 velocity; the time-dependent (Cahouet-Chabard) pressure Schur approximation for large ρ|a|Δt/η; vector and tensor unknowns for the scalar solver and the viscoelastic stress equation, which is where the DDt V_fn contract gets generalised.

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

[lmoresi]

Adversarial review (self), findings and evidence.

1. The SUPG residual is incomplete for P2 velocity. It carries the time derivative, advection, pressure gradient and body force, not the viscous term. On Kovasznay this is the whole gap to Galerkin (6.6e-4 against 1.1e-4 at 1/16, converging at 1.4 to 2.0 against 3) and it will bias any smooth, resolved flow the same way. We ship it because the term needs second derivatives the kernels do not provide and the stabilisation is meant for convection-dominated cells; a recovered Laplacian is the fix and is on the list. Users at low element Reynolds number should set supg_weight = 0.
2. The theta rule's level-n stress is exact for constant viscosity only. It is formed as 2η ε̇(uⁿ) with the current effective viscosity, so a strain-rate dependent viscosity gets the wrong level-n stress at order 1. Documented; BDF2 avoids it. Not measured.
3. The extrapolated default has a failure mode, and we have seen it. Cavity Re 400 at Courant 2: a localised alternating mode of the lagged coefficient at the lid singularity, saturating at 0.1 in the max norm while the centreline values sit still. One Picard pass or Courant 1 removes it. The default stays extrapolated because Courant 1 is where the scheme is accurate anyway, but a user at Courant 2 and element Reynolds number 8 will meet this.
4. Every linearisation gives the same steady state, so Kovasznay and the cavity cannot rank them. Only the cylinder wake did: fully implicit forms agree to three digits, the extrapolated step is 1% off in frequency and 8% on the lift peak at Courant 1. One time-dependent case, one mesh; more would be better.
5. The cylinder drag is a quarter low and does not close with the mesh (28% at resolution 20, 23% at 40 with the cylinder cells at 1/160), while the Strouhal number and pressure difference sit on the reference. The likely cause is the streamline diffusion at the scalar solver's tau weights, untuned for momentum; the Galerkin control that would isolate it cannot run at this element Reynolds number (four Newton iterations and 22 s for the first step, the second never finishing). The drag against tau is the open measurement and the first thing to do after this lands.
6. The Re 1000 cavity is not a converged comparison: 1200 steps to t = 19 at 89 to 91% of Ghia's extrema and the secondary vortex not yet formed. It shows the solver stepping at one Krylov iteration per step with geometric multigrid at element Reynolds number 16, nothing more. Two earlier attempts failed for driver reasons (a collective inside a rank-0 branch; the hang watchdog on silent ranks), not solver reasons, and the note says so.
7. cell_size() changed meaning for every consumer: the own-centroid radius is 5 to 15% larger than the kd-tree radius on irregular meshes. The Nitsche penalty scaling and the scalar SUPG τ move by that much; test_1077's serial reference moved by 4e-7 relative and was re-recorded, test_1065 now reads the new definition. get_min_radius is untouched. Anyone with a tuned Nitsche γ on an irregular mesh will see a small change.
8. Parallel coverage is two ranks on a small mesh. Four ranks ran the cavity and cylinder without incident but no four-rank number was checked against serial.

[lmoresi]

Added 9d3aa5f: Swarm.advection failed on any rank holding no particles (#693, the reshape of an empty velocity array in estimate_dt), found when releasing passive tracers at the inlet of the cylinder benchmark on four ranks. One-line fix, no behaviour change for non-empty ranks.