Zhong 2008 isoviscous response: constrained vs Nitsche free-slip behavior

[gthyagi]

Summary

While reproducing the Zhong et al. (2008) isoviscous Stokes response benchmark in UW3, constrained and Nitsche free-slip boundary conditions match different diagnostics:

• Constrained free slip provides useful built-in multiplier/topography diagnostics, but the official 8-rank 1/8 run overpredicts characteristic velocities.
• Nitsche free slip matches the Zhong Table 2 characteristic velocities well at 8 ranks, but raw boundary sigma_rr is not sufficient for dynamic-topography recovery.

This is related to the broader free-slip discussion in #104, but this issue is specifically about the Zhong 2008 spherical-shell response benchmark.

Benchmark script:

https://github.com/gthyagi/uw3-mantle-convection-benchmarks/blob/main/benchmarks/010_stokes_response_isoviscous.py

Benchmark setup

Default comparison case:

• Zhong et al. (2008) Table 2, first Stokes-flow set
• l = 2, m = 0, depth 0.5d
• r_b = 0.55, r_t = 1.0, r0 = 0.775
• uniform viscosity, eta = 1
• internal spherical-boundary load equivalent to delta(r-r0) P_l^0(theta) e_r
• response coefficients normalized by the measured angular norm of P_l^0
• serial validation: -uw_cellsize 1/4
• parallel validation: mpirun -np 8 ... -uw_cellsize 1/8

Zhong analytical references for this case:

Quantity	Reference
surface topography s	4.998e-01
CMB topography b	9.313e-01
surface characteristic velocity `	Ut
CMB characteristic velocity `	Ub

Results

Error summary

BC / run	Surface topo error	CMB topo error	Surface velocity error	CMB velocity error
constrained, serial 1/4	15.8%	16.4%	1.68%	3.88%
constrained, 8-rank 1/8	13.0%	7.66%	25.1%	63.1%
Nitsche gamma 10, serial 1/4	10.6%	15.5%	0.257%	5.87%
Nitsche gamma 10, 8-rank 1/8	48.0%	52.4%	0.654%	0.243%
Nitsche gamma 100, 8-rank 1/8	47.9%	51.5%	0.146%	0.038%

Serial actual response values

BC / serial 1/4	Surface topo	CMB topo	Surface velocity	CMB velocity
constrained	4.2090356241027566e-01	7.786127017538488e-01	1.0228577582537592e-02	1.2320478342852333e-02
Nitsche gamma 10	4.4705955177993373e-01	7.87091287177956e-01	1.0034120263339545e-02	1.2556363891469484e-02

8-rank actual response values

BC / 8-rank 1/8	Surface topo	CMB topo	Surface velocity	CMB velocity
constrained	5.648430331865263e-01	1.0026232680557212e+00	1.2584202252522738e-02	1.9340833150169143e-02
Nitsche gamma 10	2.5978694927744095e-01	4.43175114203523e-01	1.0125825484077849e-02	1.1888840166227975e-02
Nitsche gamma 100	2.6026572894429606e-01	4.517496256438572e-01	1.0074644204660118e-02	1.1855484341433173e-02

Interpretation

This appears to split into two issues:

1. Velocity response: Nitsche free slip is currently the better match to Zhong Table 2. The 8-rank 1/8 Nitsche run is below 1% error for both surface and CMB characteristic velocity. Constrained free slip overpredicts these velocities at 8 ranks.
2. Dynamic topography: constrained multiplier topography is useful, but constrained velocity response does not match the benchmark. Nitsche velocity is correct, but raw direct boundary sigma_rr is not enough for topography. This likely needs a CBF-style or projected boundary traction recovery path for Nitsche runs.

Current benchmark script behavior

010_stokes_response_isoviscous.py now defaults to:

• Zhong-style P_l^0 forcing
• measured angular-norm response normalization
• Nitsche free slip for the Zhong velocity comparison
• direct-stress topography diagnostics for non-constrained runs
• constrained free slip remains available with -uw_freeslip_type constrained for multiplier/topography diagnostics

Suggested next steps

• Investigate why constrained free slip overpredicts the 8-rank velocity response.
• Keep Nitsche as the default for Zhong Table 2 velocity-response comparison unless/until constrained velocity behavior is explained.
• Implement/test a Nitsche-compatible dynamic-topography recovery method, preferably CBF-style or projected boundary traction.
• After topography recovery is fixed, run the full Table 2 matrix for l = 2, 5, 8, 15 and depths 0.25d, 0.5d, 0.75d.

[lmoresi]

Investigation: the velocity catastrophe does not reproduce on current development (post-#242)

Built a faithful in-repo reproduction of the Zhong constrained setup — the exact test_1064 configuration (P₂⁰ radial traction via add_natural_bc on Internal, Stokes_Constrained free-slip on Upper/Lower, cellSize = 1/8, ksponly) — and ran it serial vs np=8 on current development (with #242 merged). I measured three quantities to separate the solve from the boundary extraction:

metric	serial	np=8	rel diff
volume L2 (∫v·v) — the solve	1.8260e-2	1.8195e-2	0.36%
|Ut| surface (BdIntegral on Upper) — extraction	1.0187e-2	1.0150e-2	0.37%
|Ub| CMB (BdIntegral on Lower) — extraction	1.2010e-2	1.1961e-2	0.41%

The serial values match this issue's reported serial constrained run (|Ut| 1.019e-2 vs 1.023e-2; |Ub| 1.201e-2 vs 1.232e-2), so the repro is faithful. But the np=8 catastrophe (25% / 63%) is absent — everything is partition-stable to ~0.4%.

Root cause: it was a boundary-label bug, fixed by #242

On pre-#242 development, SphericalShellInternalBoundary produced Lower (CMB) = 0 labeled points (and fragile internal labels). A free-slip BC applied to a 0-/under-populated Lower label leaves the inner surface unconstrained → velocity overpredicts there — which is exactly why #248's worst error was the CMB (63%). After #242:

• Lower is correctly labeled (160 pts at cellSize 0.3; full at 1/8), and
• the boundary areas are identical serial vs np=8: Lower 3.763621, Internal 7.510814, Upper 12.529644 (vs analytic 4πr² = 3.80 / 7.55 / 12.57, faceted-mesh approx). Partition-robust.

So the constrained solve and the characteristic-velocity extraction are both fine — the catastrophe was the internal-surface labeling, now fixed.

Recommended confirmation

Please re-run the external benchmark (010_stokes_response_isoviscous.py) at np=8 on current development. I expect the constrained velocity to now match Zhong Table 2 to ≲1% like the serial run. If it still overpredicts, we've ruled out the solve and the extraction in a faithful repro, so the remaining difference would be in the benchmark script's own mesh/normalization path.

Still open (separate from the velocity catastrophe)

Issue point 2 — dynamic-topography recovery for Nitsche (CBF-style / projected boundary traction) — is unaffected by this and remains a genuine TODO. The #242 residual diagnostics (compute_boundary_residual_fields) are now available to help with exactly that.

Repro scripts: ~/+Simulations/repro_248_internal_load.py, check_labels_np.py.

Underworld development team with AI support from Claude Code

[lmoresi]

@gthyagi — an update on the constrained-free-slip path for this benchmark, plus
a recipe for applying the internal delta(r-r0) P_l^0 e_r load so it is
partition-independent.

What changed

PR #265 (→ development) lands the parallel-correctness fixes tracked in #254.
Two of them are directly relevant to the velocities you compare here:

• False convergence is fixed. The grouped u | [p,h] Schur preconditioner
  uses inner iterative sub-solves (a variable preconditioner); the old default
  gmres is invalid for that and false-converged — it reported CONVERGED on the
  preconditioned norm while the true residual was large, landing on a wrong,
  partition-dependent answer. With fgmres + unpreconditioned stopping + a
  tightened Eisenstat–Walker it now genuinely converges and reproduces Zhong
  Table 2 (|Ut| ≈ 1%, |Ub| ≈ 0.1%).
• The rigid-rotation gauge is removed. On a fully free-slip velocity the
  rigid rotations are a true nullspace of the velocity block. The monolithic
  nullspace attached for the solve is not removed inside the fieldsplit/Schur
  iteration, so the inner velocity KSP left an unconstrained, partition-
  dependent rigid rotation in the velocity. The true residual was blind to it
  (still converged to machine precision), but it moved the tangential
  characteristic velocities |Ut|/|Ub| by ~0.1–1% serial-vs-parallel while the
  radial flow stayed clean. The solver now projects this gauge out post-solve
  (MatNullSpaceRemove, the velocity analogue of the pressure gauge). On the
  shell problem this drops the tangential serial-vs-np8 difference from 1.2% to
  ~2e-14, iteration count unchanged.

So the np=8 overprediction you were seeing should be largely accounted for by
these two — worth re-running 010_stokes_response_isoviscous.py against #265.

Applying the internal load partition-independently

There is one residual issue that #265 documents but does not fix: an
add_natural_bc on an internal surface assembles a slightly
partition-dependent load (details below). Two ways to drive the internal load:

(1) Hard fix — apply it as a body force (partition-clean, less efficient).
Realise the interface traction t(theta) e_r as buoyancy from a slightly
deflected density interface (or, equivalently, a thin radial body-force sheet of
strength t): a density jump exactly on the element-aligned interface plus a
small radial deflection delta h = t / (delta_rho g). This is a volume
integral, assembled over owned cells only, so it is partition-independent by
construction — the assembled load vector matches serial-vs-np8 to machine
precision (rel ~1e-16, vs ~3e-4 for the surface BC below). The cost is the extra
machinery to build the equivalent body force / mesh deflection (the solve
itself is no harder — same iteration count as the surface traction in our
tests); the payoff is a bit-level reproducible answer across decompositions.
Recommended when you need exact parallel reproducibility (e.g. regression
baselines).

(2) Convenient — add_natural_bc("Internal", t * e_r) (with a caveat).
This is the natural way to write the load and is fine for production accuracy,
but the assembled load carries the PETSc interior-facet partition-dependence
described below: the load vector differs serial-vs-np8 by ~3e-4. After the
elliptic Stokes solve smooths that localised seam error, the residual in the
velocity is only ~3e-6 (with the rotation gauge removed) — negligible for
most purposes, but not bit-level. Use (1) if you need the latter.

The underlying PETSc behaviour (cc @knepley)

The internal-surface natural-BC load is partition-dependent because of how the
boundary residual is assembled over interior facets. In
DMPlexComputeBdResidual_Single_Internal (src/dm/impls/plex/plexfem.c), each
labelled facet's element vector is integrated once and scattered into a single
support cell's closure:

DMPlexGetSupport(dm, point, &support);
...
DMPlexVecSetClosure(plex, NULL, locF, support[0], &elemVec[...], ADD_ALL_VALUES);

and the facet normal is oriented outward from support[0]
(DMFieldComputeFaceData_DS, dmfieldds.c). For an interior facet
(two support cells) on a non-overlapped DM, at a partition seam the locally
available support cell can be the opposite side from the one serial picks as
support[0]. The per-facet integral value is the same, but it is scattered
through a different cell closure; closure DOFs that are non-owned on that rank
are added to the local vector and then dropped at global assembly — so the
assembled internal-surface load under-counts by a partition-dependent amount
(~3e-4 for the l=2 shell load here; the force function integral / RMS are
identical to machine precision, so it is purely an assembly-side effect).

We keep the FE-assembly DM non-overlapped on purpose (a 1-cell overlap
double-counts volume assembly via LocalToGlobal + ADD_VALUES), and from the
application side we cannot impose a deterministic support[0] (support order
follows the DMPlex point partition). Things we tested that do not help:
owner-only relabelling of the facets (the ghost copies already contribute
nothing — PETSc integrates each facet once), and a 1-cell partition overlap
(verified it propagates but the assembled load is bit-identical, so it does not
recover the dropped DOFs). A clean fix really wants the interior-facet boundary
residual to scatter into the owned support cell (or a support-choice that
is partition-invariant) — DMSetNullSpaceConstructor-style, this is below the
application layer.

Flagging it here so it registers with the PETSc maintainers, since the
partition-dependent internal-facet support[0] scatter is a general
DMPlexComputeBdResidual property, not specific to UW3.

Underworld development team

[knepley]

Comment on "The rigid-rotation gauge is removed": The projection is fine. Another way to do this is to use Multigrid, and solve the coarse grid using PCSVD. This will project out the nullspace, and usually is fine. Sometimes you have to make sure that the prolongator does not put nullspace components back in.

[knepley]

@lmoresi I think I understand the problem (but that might be an illusion since I have not run it myself). However, how does this sound. I added some functionality for PyLith to orient an internal surface (this is necessary for faults). You can call DMPlexOrientLabel() where the label marks mesh point with their depth, and it will reorient the faces to make a continuous surface, changing the cell attachemnts as necessary. I think this will give you what you want.

[lmoresi]

Status after this session's merges: #242 fixed the SphericalShellInternalBoundary label bug (the main constrained velocity error), and #265 fixed the parallel-correctness residuals + added a guarded auto pressure gauge and gauge-invariant topography(reference='mean'). The constrained overprediction at 8 ranks should now be substantially reduced — needs an np=8 Zhong re-run to confirm (tracked with #254). The Nitsche raw sigma_rr → dynamic-topography recovery point is separate and still open.

Underworld development team with AI support from Claude Code

[lmoresi]

Rotated free-slip on this benchmark: the Schur blow-out is diagnosed and fixed (PR #306)

Running this benchmark's configuration with add_rotated_freeslip_bc (serial 1/8, tol 1e-7) exposed two independent defects in the rotated path:

| method | outer its | Ut err | Ub err | rotation content max|L| |
|---|---:|---:|---:|---:|
| Nitsche γ=10 | 1 | 0.48% | 0.63% | 1.6e-4 |
| constrained | 5 | 0.11% | 0.37% | 1.6e-4 |
| rotated (before) | 44 | 0.31% | 0.14% | 6.7e-4 |
| rotated (after, PR #306) | 1 | 0.23% | 0.15% | 1.6e-4 |

np=8 reproduces the serial values with 1 outer iteration.

1. Preconditioning: the rotated solve never assembled the Pmat, so the native 1/μ pressure-mass Schur preconditioner (schur_precondition=a11) was silently replaced by selfp+Jacobi — weak on curved boundaries and poisoned by the unit constraint diagonals along the whole spherical surface. Fixed by installing the native mass block as a USER Schur preconditioner (pressure is untouched by the rotation, so no remapping is needed).
2. Nullspace (the correctness item): the rotated nullspace handling was 2D-only. A spherical shell with free-slip on both boundaries admits all three rigid rotations; leaving them on the operator both grinds the Krylov solve and pollutes the velocity with arbitrary rotation content (visible as the elevated angular momentum above). All three e_k×r modes are now constraint-checked, attached, and projected out.

Also fixed on the way: the custom geometric-FMG velocity block never worked on a closed shell (Galerkin coarse operator inherits the rotation nullspace → LU zero pivot, outer reason -11) — its coarse solve is now SVD; and analytic normals given as mesh.CoordinateSystem.unit_e_0 now work directly (previously a NameError — they are unwrapped to Cartesian before compilation).

Underworld development team with AI support from Claude Code

[lmoresi]

Triage status (2026-07-23): every thread of this issue now has a landed fix — proposing we close after a confirming re-run

Summary of where the three original observations stand:

1. Constrained free-slip velocity overprediction at 8 ranks — root cause was the broken SphericalShellInternalBoundary Lower/Internal label (Fix spherical shell labels and add Stokes residual diagnostics #242, merged), with the parallel-correctness residuals (false convergence, gauge, knockout) fixed in Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 (merged). The remaining ~3e-6 partition-dependent residual on the 3D internal-boundary case is tracked in Stokes_Constrained: parallel-correctness residuals (partition-dependent [p,λ] gauge + constraint-assembly drift + lossy knockout) #254 as an upstream PETSc limitation.
2. Rotated free-slip on this benchmark — the Schur blow-out and missing 3D rotation nullspace were fixed in fix(rotated-bc): Schur preconditioning parity + 3D rotation nullspace — Zhong #248 blow-out (44 its -> 1) #306 (44 outer iterations → 1, rotation content back to the 1.6e-4 baseline).
3. Nitsche raw sigma_rr insufficient for dynamic topography — rather than repairing the raw-stress route, the recommended path is now rotated strong free-slip, where the constraint reaction gives the boundary normal traction directly (solver.boundary_normal_traction(boundary) / solver.dynamic_topography(...)), machine-precision v·n = 0 included.

@gthyagi — when you have a chance, could you re-run the official np=8 1/8 Zhong comparison on current development? If the constrained velocities now match Table 2 (and the rotated path gives you the topography diagnostics you wanted), I think this issue can close, with #254 carrying the last parallel residual.

Underworld development team with AI support from Claude Code

[gthyagi]

@lmoresi Confirming re-run completed against current development (e6923bbe) using the updated Zhong benchmark script. All runs below are 8-rank MPI; values are nondimensional. Errors are relative to Zhong et al. (2008).

No-self-gravity response (Table 4)

BC / resolution	Wall / outer its	U_t	U_b	surface topo	CMB topo	surface geoid	CMB geoid
Zhong	-	0.01006	0.01186	0.41920	0.77060	0.02579	0.03206
constrained, 1/8	302.25 s / 3	0.01007 (+0.12%)	0.01182 (-0.37%)	see note	see note	see note	see note
Nitsche, 1/8	25.70 s / 1	0.01011 (+0.48%)	0.01179 (-0.63%)	0.41407 (-1.22%)	0.70068 (-9.07%)	0.02349 (-8.93%)	0.02406 (-24.95%)
rotated, 1/8	18.43 s / 1	0.01008 (+0.23%)	0.01184 (-0.16%)	0.41291 (-1.50%)	0.79348 (+2.97%)	0.02495 (-3.24%)	0.03420 (+6.67%)
rotated, 1/16	67.33 s / 1	0.01008 (+0.16%)	0.01187 (+0.07%)	0.41740 (-0.43%)	0.76224 (-1.09%)	0.02528 (-1.97%)	0.03103 (-3.20%)

Self-gravity response (Table 2)

BC / resolution	surface topo	CMB topo	surface geoid	CMB geoid
Zhong	0.49980	0.93130	0.04486	0.05461
Nitsche, 1/8	0.48637 (-2.69%)	0.82437 (-11.48%)	0.04021 (-10.36%)	0.04204 (-23.02%)
rotated, 1/8	0.49168 (-1.63%)	0.96296 (+3.40%)	0.04381 (-2.34%)	0.05761 (+5.49%)
rotated, 1/16	0.49638 (-0.68%)	0.91801 (-1.43%)	0.04393 (-2.08%)	0.05295 (-3.04%)

Conclusions

• The constrained 8-rank velocity overprediction is fixed: U_t and U_b are now within 0.12% and 0.37% of Zhong, respectively. This confirms the Fix spherical shell labels and add Stokes residual diagnostics #242/Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 diagnosis.
• The rotated path confirms fix(rotated-bc): Schur preconditioning parity + 3D rotation nullspace — Zhong #248 blow-out (44 its -> 1) #306: one outer iteration, fastest 1/8 solve, and the best complete topography/geoid agreement. Refining from 1/8 to 1/16 reduces all reported topography/geoid errors to about 3.2% or less.
• Topography nuance: constrained topography(reference="mean") returned surface/CMB values -0.41781 / 0.76629, so its surface sign is not valid for this paper comparison. A direct rotated-reaction projection was also tested but did not recover the internal-load topography amplitude. The validated Nitsche/rotated values above therefore use spherical_shell_dynamic_response(..., topography_source="auto"), which selects the CBF-lumped residual path.
• Thus, the original constrained velocity issue is resolved. The benchmark validates rotated free slip for the solve, but it does not yet independently validate the direct rotated-reaction topography path for this internal-boundary load.

Full running notes and earlier comparisons are in the benchmark results document.

[gthyagi]

@lmoresi Follow-up on the direct rotated-free-slip topography path: the earlier near-zero result was a real 3D recovery bug, and the corrected reaction-based result now matches Zhong closely.

Root cause

stokes.boundary_normal_traction(...) uses the shared boundary-flux _desmear routine. That routine assembled a boundary mass in 2D, but its dim != 2 fallback only mean-removed the raw nodal reaction. In 3D it therefore returned a nodal load, not pointwise sigma_nn: no division by surface mass was performed. This explains both the near-zero amplitude and its mesh dependence.

The rotated Stokes solve and constraint reaction itself were correct. The defect was entirely in the 3D reaction-to-traction conversion.

Fix

The proposed fix:

• subtriangulates each P2 triangular boundary face into four positive-area triangles;
• globally deduplicates those triangles and assembles a lumped surface mass;
• computes sigma_nn = nodal reaction / nodal surface mass;
• removes the surface-mass-weighted mean;
• adds rotated_reaction to spherical-shell topography dispatch.

With topography_source="auto", postprocessing now uses:

• stokes.topography(...) for constrained multipliers;
• stokes.boundary_normal_traction(...) for rotated strong free slip;
• CBF residual recovery for the remaining BC paths.

The existing Zhong geoid and self-gravity formulas then consume the selected surface/CMB topographies; users do not need BC-specific geoid code.

Direct rotated-reaction results

Run	Surface topo	CMB topo	Surface geoid	CMB geoid	Surface topo SG	CMB topo SG	Surface geoid SG	CMB geoid SG
Zhong	0.41920	0.77060	0.02579	0.03206	0.49980	0.93130	0.04486	0.05461
serial 1/4	0.41541	0.75294	0.02471	0.02989	0.49250	0.90324	0.04288	0.05109
serial 1/8	0.42032	0.76698	0.02595	0.03173	0.50136	0.92634	0.04508	0.05416
np=8 1/8	0.42033	0.76694	0.02595	0.03173	0.50138	0.92628	0.04508	0.05416
np=8 1/16	0.41905	0.76894	0.02573	0.03187	0.49949	0.92873	0.04475	0.05432

For np=8 1/8, the no-self-gravity topography/geoid errors are +0.27%, -0.47%, +0.63%, and -1.03%. All four self-gravity errors are below 0.83%.

Refining to np=8 1/16 improves every response family. The no-self-gravity errors are -0.04%, -0.22%, -0.23%, and -0.59%; the corresponding self-gravity errors are -0.06%, -0.28%, -0.26%, and -0.54%. The run completed in 65.68 s with one outer iteration.

At the same 1/8 mesh, serial and np=8 differ by only 0.003% for surface topography, 0.006% for CMB topography, and at most 0.013% for geoid. The corrected 3D recovery is therefore partition-stable for this case.

Validation:

• level-2 postprocessing suite: 4 passed;
• existing rotated-topography and boundary-flux tests: 2 passed;
• serial 1/4, serial 1/8, np=8 1/8, and np=8 1/16 benchmark runs all passed with one outer iteration and topography_source=rotated_reaction.

Implementation commit: gthyagi@02b181f5

Detailed benchmark record: https://github.com/gthyagi/uw3-mantle-convection-benchmarks/blob/main/docs/results_bench_010_zhong2008_isoviscous_response.py.md

This supersedes the earlier conclusion that direct rotated reaction did not recover the Zhong-scale amplitude; that result exposed the missing 3D surface-mass de-smear.

[lmoresi]

Evidence relevant to closing this: the #438 pre-test measured the two free-slip implementations directly under anisotropy — the rotated strong constraint is exact (wall-normal flow ~1e-16, direct-LU control) while Nitsche-as-implemented leaks 13% at 100x anisotropy; in the isotropic limit the gap collapses to the known ~1e-3 Nitsche leak. That quantifies the constrained-vs-Nitsche behaviour difference this issue tracks: it is the weak-imposition leak, scaling with anisotropy, not a defect in either solve. The Zhong isoviscous response itself is exercised in test_1018 (SolCx topography via the constraint reaction). Suggest: close as answered by that evidence, or say what remains and it can ride the workstream-C benchmark pass (spherical Zhong is on that acceptance list).

Underworld development team with AI support from Claude Code