Document Stokes constant-penalty compatibility after viscosity-scaled penalty change

[gthyagi]

Summary

Commit a388902b80754654e20cab572f63b72f4483f6dd (2026-06-10, Stokes: viscosity-scale the augmented-Lagrangian (grad-div) penalty) changed the Stokes augmented-Lagrangian incompressibility penalty from a bare constant penalty to a viscosity-scaled penalty.

Conceptually this changed the Stokes operator from:

stress + penalty * div_u * I

to:

stress + penalty * viscosity * div_u * I

This is physically sensible for spatially variable viscosity because the penalty then scales with the local stress scale. However, user scripts and benchmark scripts that assumed the old constant-penalty semantics can silently change behaviour and performance unless they are updated.

Why this needs stronger user documentation

The current docstrings describe the new viscosity-scaled penalty, but users also need a clear migration note explaining:

1. stokes.penalty is now a dimensionless viscosity-scaled multiplier.
2. stokes.penalty = C no longer means the old constant penalty C; it means C * viscosity in the operator.
3. Benchmarks/scripts that require old constant-penalty behaviour should use the compatibility form:

# old constant penalty lambda_const under current UW3
stokes.penalty = lambda_const / viscosity
stokes.saddle_preconditioner = 1.0 / (viscosity + lambda_const)

4. For the current viscosity-scaled formulation, users should usually start with the default stokes.penalty = 0; if a penalty is used, saddle_preconditioner should usually be left unset or kept consistent with the viscosity-scaled operator.
5. Examples and benchmarks that use patterns like stokes.penalty = 1.0 together with 1 / (viscosity + stokes.penalty) should be reviewed, because that expression mixes the old constant-penalty interpretation with the new viscosity-scaled operator.

Benchmark evidence

This affected a 2D BH92-style fault benchmark with strong viscosity contrast and localized weak/fault-related structure.

Case tested:

• n = 3
• lc = 0.02
• faultDist_fact = 2.5
• viscosity contrast = 1e5
• 8 MPI ranks
• transverse-flow BH92 no-fault benchmark variant

Observed timings:

case	total time	GMRES orthogonalizations	notes
archived old UW3 / old constant penalty	50.49 s	96	reference behaviour
current UW3 using unmodified old script settings	574.8 s	723	much slower; changed penalty semantics
current UW3 with constant-penalty compatibility form	45.83 s	100	matches old timing and fields
current UW3 viscosity-scaled penalty + default pyx solver, reduced lc=0.04 diagnostic	172.6 s	230408	much slower diagnostic case

Field comparison for the fixed current-UW3 constant-penalty compatibility run against the archived old run was very close:

• velocity relative L2: 1.746e-05
• pressure relative L2: 1.458e-04
• viscosity relative L2: 7.744e-05
• strain-rate invariant relative L2: 1.131e-04
• stress invariant relative L2: 1.871e-04

So for this benchmark, the compatibility formula restores both the old timing behaviour and the old solution fields.

Suggested action

Please add a user-facing documentation section or release-note entry explaining the difference between:

• constant penalty: historical lambda_const * div(u) behaviour
• viscosity-scaled penalty: current lambda * viscosity * div(u) behaviour

and provide explicit recipes for both.

Also please review/update examples and benchmark scripts that set stokes.penalty and stokes.saddle_preconditioner manually, especially where the preconditioner still assumes 1 / (viscosity + penalty) under the old semantics.