Schur preconditioner scaling omits the solvent and DEVSS viscosities, and the pressure solve runs away at high Wi

[lmoresi]

Symptom

Creeping Oldroyd-B flow past a confined cylinder stops making progress above about Wi 0.2. The step cost is flat and then cliffs, and the solution degrades on the same step:

step	wall time	drag (reference 120.59)
13	70.8 s	114.37
14	77.2 s	113.48
15	1435.7 s	86.07

No iteration cap is hit. Newton converges in one iteration on the absolute residual and the outer Krylov solve converges on relative tolerance in one to four iterations, every step, including the bad one.

Where the time goes

-log_view, 15 steps, 2118 s total:

event	time	calls	per call
KSPSolve_FS_Schu (pressure)	1990.7 s	26	76.6 s
KSPSolve_FS_0 (velocity, FMG)	10.1 s	26	0.39 s
SNESJacobianEval	67.3 s	164	0.41 s
DEVSS projection	66.8 s	32	2.09 s
PCSetUp	0.8 s	164	0.005 s

94% in the Schur (pressure) solve; the velocity multigrid is 0.5% and is not the problem. Underneath: 202,871 MatMult and 709,321 KSPGMRESOrthog calls, so the pressure Krylov solve is running to enormous iteration counts.

Tolerances are as designed (outer 1e-6, pressure 1e-7, velocity 3.3e-8), one order apart, so this is not a tolerance problem.

Likely cause

The pressure block's Schur preconditioner is the viscosity-scaled pressure mass, _pp_G0 = 1 / constitutive_model.K (petsc_generic_snes_solvers.pyx:8044). For the viscoelastic model,

K = dt * eta_p * mu / (dt * mu + eta_p * c0)   ==   eta_p * dt / (lambda + dt)

which is the polymer effective viscosity only. The momentum operator it is meant to precondition actually carries

eta_operator = eta_p * dt / (lambda + dt)  +  eta_s  +  eta_a(DEVSS)

Neither the solvent viscosity nor the DEVSS artificial viscosity appears in K. The mismatch is the ratio eta_operator / K, which grows with Wi and shrinks with dt. For the confined-cylinder benchmark (beta = 0.59, so eta_p = 0.41, eta_s = 0.59):

Wi	dt	K (what the pc uses)	true operator viscosity	ratio
0.2	0.4	0.273	0.863	3.2
0.4	0.4	0.205	0.795	3.9
5	0.0125	0.00102	0.591	578

An Oldroyd-B fluid with a significant solvent fraction is exactly the case where this bites, and the cylinder and contraction benchmarks are both in it. DEVSS adds a further eta_a the preconditioner also does not see.

Suggested direction

1. Include the solvent and DEVSS viscosities in the Schur scaling, so _pp_G0 uses the viscosity the momentum operator actually carries rather than the polymer part alone. Cheap to try: set stokes.saddle_preconditioner = 1/(K + eta_s + eta_a) on a failing case and see whether the pressure iteration count collapses.
2. Grad-div augmentation (stokes.penalty) as the other standard lever on Schur conditioning.
3. If neither is enough, the Schur approximation itself may need revisiting for the viscoelastic operator.

Reproducer

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

UW_OB_DEBUG=1 python ob_cylinder.py -uw_wi 0.4 -uw_dt 0.4 -uw_t_end 8 -uw_out run

Half domain, blockage 0.5, beta 0.59, Re 0, 11744 cells, three multigrid levels. Fails at step 15. -uw_penalty and -uw_p_rtol are already wired for experiments. The Newtonian case of the same rig reproduces the benchmark drag to 0.04%, so the rig itself is sound.

Related: #754 (DEVSS not subtracting its own viscosity on a plain Stokes solve).

[lmoresi]

Correction to the issue body: it is the velocity block, not the pressure solve.

-log_view books the inner velocity solves to KSPSolve_FS_Schu, because they happen inside the Schur complement application. Counting the convergence reports directly instead:

inner velocity solve	count	mean iterations
hit the 200 cap	2496	200.0
converged	2336	86.3

More than half fail outright and the rest take 86 multigrid-preconditioned iterations. The 26 outer velocity solves at 0.39 s each, which the body quotes as evidence the multigrid is fine, are a tiny minority.

The K hypothesis in the body is also withdrawn. The material parameters are uniform, so getting that scalar wrong is a uniform rescaling of the Schur preconditioner, and a Krylov method with a relative tolerance is invariant to it. It would not have cost an iteration.

What is actually wrong is structural: the viscoelastic velocity operator is not a scaled Laplacian. The upper-convected terms add velocity-gradient-dependent stretching, the operator loses symmetry, and the elastic stress layer is thinner than a cell near the wall. SOR has no purchase on that at any sweep count.

Partial remedy, now on development via #757. solver.strategy = "bulletproof" (gmres/8 + additive Schwarz smoothing):

smoother	the failing step	drag
sor	26947 s	−4642
asm	115 s	86.05

Capped fraction 90% → 77%. It survives, and it does not fix the underlying problem: the same wrong drag, 86.05, appears at the same step under both smoothers. Two different iterative paths reaching the same wrong number says the fault is outside the linear algebra.

Refinement does not rescue it either. Halving the cylinder cell (0.04 → 0.02) makes it fail earlier, at step 14, with the velocity solve hitting 200 iterations, the pressure residual then growing, and an overflow to not-a-number inside the failed solve. Every stored field is finite and below 1e8 at the end of every completed step, so nothing upstream is injecting the not-a-number; it is generated by the divergence.

So the open question for whoever picks this up is narrower than the body suggests: what makes the viscoelastic velocity operator unsolvable at Wi ≳ 0.3 on this geometry, given that the elastic layer is already sub-cell (half-height 0.0318R against a 0.04R cell at Wi 0.2) and refining it makes the linear algebra worse rather than better. Patch-based/Vanka smoothing inside a monolithic (u, p, tau) multigrid is the direction the literature takes; augmented Lagrangian is the cheaper thing to try first.