Fix spherical shell labels and add Stokes residual diagnostics

[gthyagi]

Summary

This PR adds the UW3 support needed to reproduce the spherical-shell Stokes response benchmarks used in Zhong et al. (2008):

• Fix uw.meshing.SphericalShellInternalBoundary() so Lower, Internal, and Upper labels are valid DMPlex boundary labels.
• Add low-level Stokes residual diagnostic hooks used for CBF-style dynamic-topography recovery.
• Add constrained/free-slip spherical-shell response regressions that document the current constrained grouped-Schur limitation.
• Keep the production Stokes/Nitsche solve path unchanged; the new residual hooks are diagnostics around the assembled problem.

Benchmark Theory

Zhong et al. (2008) compare Stokes response functions in a spherical shell with:

• inner radius r_b = 0.55, outer radius r_t = 1.0;
• radial delta forcing delta(r-r0) P_l^0(theta) e_r;
• free slip on the CMB and surface;
• response quantities at the surface and CMB: characteristic horizontal velocity, dynamic topography, and geoid.

In UW3 the delta forcing is represented as the equivalent radial natural traction on an embedded internal spherical boundary at r0. This requires SphericalShellInternalBoundary() to produce three usable surface labels:

• Lower: CMB boundary;
• Internal: load boundary;
• Upper: surface boundary.

The velocity response is a standard Stokes solve result. Dynamic topography is more delicate: direct pointwise/projected sigma_rr is not the same as the consistent-boundary-flux style stress recovery used by strong free-slip benchmark codes. The residual diagnostic hooks in this PR allow the benchmark script to recover boundary reaction/topography from selected Stokes residuals instead of relying only on direct stokes.stress post-processing.

Nitsche Free Slip: What Changed and Why It Matters

This PR does not change UW3's production Nitsche weak form.

The important update is diagnostic and validation support around the Nitsche path:

• The benchmark uses Nitsche free slip as the validated production path because it reproduces Zhong characteristic velocities robustly in serial and MPI runs.
• Direct Nitsche sigma_rr post-processing was not sufficient for topography/geoid comparison.
• The new residual hooks expose:
  • full FEM residual diagnostics;
  • boundary-specific weak residual diagnostics for registered natural/Nitsche boundary forms;
  • cloned-DM volume-only residual diagnostics for CBF-style recovery without appending registered boundary residuals.

This separation matters because a correct velocity solve and a correct topography recovery are related but not identical checks. The benchmark now confirms that Nitsche gives the correct velocity response, then uses residual-based recovery plus Zhong Appendix A post-processing for dynamic topography and geoid.

What Changed

Mesh

• Reworked SphericalShellInternalBoundary() geometry construction so the physical shell volume and the embedded internal spherical surface are preserved cleanly.
• Added radius-order validation for radiusInner < radiusInternal < radiusOuter.
• Added regression coverage for Lower, Internal, and Upper boundary areas.

Residual Diagnostics

• Added Stokes FEM residual field extraction.
• Added boundary weak residual extraction through the same boundary weak form registered by UW3.
• Added cloned-DM volume-only residual extraction. The cloned DM copies sections, fields, equations, constants, and auxiliary vectors but has no registered boundary objects, so PETSc does not append boundary residuals to the selected-cell volume residual.
• Fixed the residual diagnostic context signature to use void *, matching the Cython declaration and allowing all Cython extension modules to build.

Regression Tests

• Added tests/test_1064_constrained_spherical_shell_response.py.
• The test records:
  • Nitsche/default matches the Zhong velocity scale for the low-resolution response case;
  • direct-LU Nitsche and direct-LU constrained are useful same-path diagnostics but are not the benchmark reference;
  • constrained field-split variants that do not reproduce the validated Nitsche/default response are strict expected failures.

Benchmark Results

External benchmark repository:

gthyagi/uw3-mantle-convection-benchmarks

Benchmark script:

benchmarks/bench_010_zhong2008_isoviscous_response.py

Configuration:

• Nitsche free slip;
• Zhong P_l^0 internal-boundary forcing;
• volume_residual_cbf_lumped topography recovery;
• Zhong Appendix A self-gravity/geoid post-processing;
• Table 2 sweep over l = 2, 5, 8, 15 and depths 0.25d, 0.5d, 0.75d.

Default l=2, 0.5d Case

Zhong Table 2 analytical values:

Quantity	Analytical	UW3 8-rank 1/16	Error
surface velocity	1.006e-02	1.0076094909651055e-02	0.16%
CMB velocity	1.186e-02	1.186556831000749e-02	0.05%
surface topography, self-gravity	4.998e-01	4.9577940809588633e-01	-0.80%
CMB topography, self-gravity	9.313e-01	9.156930209567329e-01	-1.68%
surface geoid, self-gravity	4.486e-02	4.376413039396167e-02	-2.44%
CMB geoid, self-gravity	5.461e-02	5.265637036600948e-02	-3.58%

Full Table 2 Sweep, 8-rank 1/16

Errors relative to Zhong Table 2 analytical values:

depth	l	s err %	b err %	ht err %	hb err %	Ut err %	Ub err %
0.25	2	-0.75	-1.67	-3.31	-3.43	-0.01	-0.05
0.25	5	-0.71	-3.01	-1.70	-3.96	-0.06	-0.60
0.25	8	-1.28	-5.99	-2.43	-7.21	-0.19	-1.97
0.25	15	-4.18	-17.32	-6.30	-19.40	-0.90	-7.07
0.5	2	-0.80	-1.68	-2.44	-3.58	0.16	0.05
0.5	5	-1.08	-3.10	-1.87	-4.49	-0.13	-0.56
0.5	8	-2.05	-5.94	-2.94	-7.73	-0.38	-1.81
0.5	15	-6.28	-16.73	-7.75	-19.63	-1.74	-6.68
0.75	2	-0.76	-1.62	-2.42	-4.88	0.22	-0.02
0.75	5	-0.93	-2.34	-1.52	-4.33	-0.15	-0.27
0.75	8	-1.99	-4.48	-2.60	-7.06	-0.67	-1.13
0.75	15	-6.56	-15.49	-7.52	-20.54	3.19	-5.82

Interpretation:

• l <= 8 is reproduced well at 1/16: velocity errors are below about 2%, topography errors are mostly below 6%, and geoid errors are mostly below 8%.
• l = 15 remains resolution-sensitive, especially for CMB topography/geoid. This is now a convergence-resolution issue rather than a free-slip or forcing-definition issue.

Are All Recent Commits Essential?

For a single PR whose goal is “support Zhong-style spherical-shell benchmark reproduction”, yes, the commits are connected:

• SphericalShellInternalBoundary fixes are essential. Without them the benchmark cannot reliably integrate/apply loads on Lower, Internal, and Upper.
• Residual diagnostic hooks are essential for the validated topography/geoid path. Direct stress recovery was insufficient; CBF-style residual recovery is what made the Table 2 comparison work.
• The constrained response regression is not needed to run Nitsche itself, but it is important evidence. It documents why Nitsche/default is the current benchmark reference and why constrained grouped-Schur results should not be used as the paper comparison until that solver path is fixed.
• The final void * context-signature fix is essential because the shared compatibility header must build in all Cython extension modules.

If maintainers prefer smaller PRs, this branch can still be split into two review units later:

1. mesh-label fix for SphericalShellInternalBoundary();
2. residual diagnostics plus constrained-response regression.

Validation

Local validation in this branch:

• ./uw build passed.
• pixi run -e amr-dev pytest -q tests/test_0502_boundary_integrals.py::test_bd_integral_spherical_internal_boundary_areas passed: 1 passed in 4.70s.
• pixi run -e amr-dev pytest -q tests/test_1064_constrained_spherical_shell_response.py passed with expected failures: 2 passed, 2 xfailed in 180.38s.

Representative fixed serial boundary areas at ri=0.55, rint=0.775, ro=1.0, cellSize=0.25:

• Lower = 3.651738570092663;
• Internal = 7.402572621636933;
• Upper = 12.418895120898583.

Before the mesh fix, the same benchmark path reported Lower = 0.0 and failed the boundary-area check.

[lmoresi]

Review — strong work; one regression to resolve before merge

Reviewed in two parts (C/Cython diagnostics API, and the spherical meshing + tests). The boundary-label fix is correct and backward-compatible, and the new residual-diagnostics API is additive-only and memory-clean against the custom PETSc 3.25 build target. Two items should be resolved before this lands on development for the release, plus one C nit.

✅ What's solid

• Boundary labels: Lower=11 / Internal=12 / Upper=13 names and values unchanged; surface binning is done on analytic OCC geometry (bbox_radius), and the duplicate-inner-surface disambiguation via getBoundary(shell_vol) is sound. All three in-repo spherical benchmark examples keep working.
• Centre footgun: the new generator no longer materialises a Centre physical group, so the known "Centre pseudo-label hard-aborts PETSc" issue is not triggered (safer than before).
• Residual API (compute_volume_residual_fields / compute_boundary_residual_fields): purely additive — no existing solve/_build/assembly line changed. Cython resource lifecycle is correct (vecs restored in finally, every getSubVector/restoreSubVector paired, IS destroyed via created_is_field gate). C borrowed objects (sections, DS, aux) correctly not destroyed; signatures match custom PETSc 3.25. Boundary form-key matches the real BC assembly path (UW_Boundaries label + boundary_label_val).

⚠️ Should-fix before merge

1. Silent loss of Inner/Outer region extraction. The new single-volume cut+embed approach creates only the Elements volume group — the Inner=101/Outer=102 region groups are gone. But spherical.py still runs new_mesh.regions = regions, so the mesh advertises regions that don't exist on the DM. mesh.extract_region("Inner") would filter on a missing label → empty/failed submesh. Commit 49845c8 just added this capability "+ extract_region verified", so this silently undoes it. Nothing in-repo uses it on this mesh (only AnnulusInternalBoundary does), so the suite stays green — but advertising a broken capability is a footgun. Please either drop new_mesh.regions = regions (and the regions enum) for this generator, or document in the docstring + CHANGELOG that the internal-boundary spherical shell no longer supports Inner/Outer region extraction (inherent to the embed design).

2. Mis-tiered test. test_0502::test_bd_integral_spherical_internal_boundary_areas is a genuinely good area-integral validation, but it inherits the module-level pytestmark = [level_1, tier_a]. It builds a full 3D gmsh+embed mesh — neither "seconds/no-solve" (level_1) nor yet production-soaked (tier_a, CI/TDD-trusted). Please override to @pytest.mark.level_2 / tier_b until the new generator has soaked.

NIT

• petsc_compat.h UW_DMPlexComputeResidualByKeyVolumeOnly: the cloned vdm leaks one DM if DMPlexComputeResidualByKey errors — the PetscCall early-returns before DMDestroy(&vdm). It's a per-call diagnostic error path, but trivial to harden (destroy vdm before returning / cleanup goto).

Happy to push (1) and (2) myself if you'd prefer — just say the word. Leaving open pending the region-advertisement decision.

Underworld development team with AI support from Claude Code

[lmoresi]

Updated and ready — region regression fixed properly

Pushed three changes to address the review and bring this current (maintainer-edit), per Louis's direction:

1. Inner/Outer regions now materialised (not removed). extract_region is the sub-meshing entry point, so the fix is to make the regions real. Since the cut+embed design is a single OCC volume (gmsh can't emit two region groups), the Inner/Outer cell DM labels are created after import by classifying each cell's centroid radius vs radiusInternal — exact (the embedded internal surface is conformal, so no cell straddles it) and rank-local / partition-independent. Verified end-to-end: extract_region("Inner") → r∈[0.55, 0.775], extract_region("Outer") → r∈[0.775, 1.0].
2. Re-tiered the new 3D area test to level_2/tier_b (it builds a full 3D gmsh+embed mesh; was inheriting module level_1/tier_a).
3. Merged current development (13 commits incl. Use PETSc-native XDMF visualization topology #218/Per-iteration SNES update callbacks (pressure gauge + boundary-correct generic scatter) #250/fix(meshing): extract_region no longer raises (#197) #253). One .pyx conflict — your residual-diagnostics methods vs Per-iteration SNES update callbacks (pressure gauge + boundary-correct generic scatter) #250's callback methods, both added to SNES_Stokes_SaddlePt; resolved by keeping both (distinct methods).

Validated on the merged branch: extract_region Inner/Outer correct; Lower boundary now 160 pts (was 0 on development); test_1016 (callbacks), test_0771 (extract_region), test_0502 spherical-area, test_1010 (Stokes core) all green.

Remaining NIT from review (not blocking): the vdm error-path leak in UW_DMPlexComputeResidualByKeyVolumeOnly.

Underworld development team with AI support from Claude Code

[lmoresi]

Approved. Internal-surface labels fixed (Lower 0→160 pts), Inner/Outer region extraction now works end-to-end, residual diagnostics are additive + memory-clean. Region regression and test-tier addressed; merged current development (conflict resolved keeping both diagnostics + callback method sets). See detailed comment above.