Spherical Thieulot benchmark: remaining pressure mismatch at the inner boundary

[gthyagi]

This is a follow-up on the spherical Thieulot benchmark. Both m = -1 and m = 3 cases are now running without the earlier implementation bug. The earlier m = 3 failure was traced to an incorrect term in the manufactured body force in the benchmark script and is now fixed.

The remaining issue is the pressure mismatch at the inner boundary.

Reproduction scripts:

• ex_stokes_thieulot.py
• thieulot_field_plots.py

Tests performed:

• Used the essential BC solution as the baseline/reference.
• Compared against:
  • natural_full
  • natural_normal_petsc
  • natural_normal_analytic
  • natural_normal_projected
• Ran the comparison for m = -1 and m = 3 at cellsize = 1/8, Q2/Q1, p_bc = False.
• For m = 3, also retested the PETSc-normal case with a much larger BC penalty and a tighter Stokes tolerance.

Observed outcomes:

• natural_full stays very close to the essential baseline:
  • m = -1: p_inner_l2(base) = 0.023663
  • m = 3: p_inner_l2(base) = 0.00880887
• The normal-only BC does not improve the inner-boundary pressure:
  • m = -1:
    • PETSc normals: p_inner_l2(base) = 2.23654
    • Analytical normals: p_inner_l2(base) = 12.985
    • Projected normals: p_inner_l2(base) = 13.2424
  • m = 3:
    • PETSc normals: p_inner_l2(base) = 0.805627
    • Analytical normals: p_inner_l2(base) = 3.60611
    • Projected normals: p_inner_l2(base) = 3.87889
• For m = 3, increasing the PETSc-normal penalty to 1e10 or tightening the Stokes tolerance to 1e-8 did not materially change the PETSc-normal result.

Interpretation:

• Pressure appears highly sensitive at the inner boundary, likely because of the strong radial dependence there (roughly r^-n behavior).
• Swapping PETSc facet normals for analytical normals does not fix the remaining pressure mismatch at the inner boundary.
• In these tests, analytical/projected normals actually make the normal-only case worse relative to the essential baseline.
• This suggests the dominant effect is not just facet-normal accuracy. The bigger change is the BC class itself: natural_normal behaves like a free-slip-style normal-only condition, whereas essential / natural_full enforce the full benchmark velocity.

cc @lmoresi for visibility.

[gthyagi]

Spherical Thieulot BC-Normal Debug

Goal

Test whether switching the normal definition in a natural normal-only boundary
condition improves pressure on the inner spherical boundary, using the
essential solution as the reference baseline.

Driver

Standalone script:

benchmarks/spherical/ex_stokes_thieulot_bc_normals.py

The main benchmark driver
benchmarks/spherical/ex_stokes_thieulot.py
was not modified for this experiment.

For each parameter set, the script runs:

• essential
• natural_full
• natural_normal_petsc
• natural_normal_analytic
• natural_normal_projected when requested

and records:

• global pressure L2 vs analytical pressure
• inner / outer boundary pressure L2 vs analytical pressure
• inner / outer boundary pressure L2 vs the essential baseline
• inner / outer tangential velocity RMS

The key normal-only term is:

penalty * ((n . v_num) - (n . v_ana)) * n

Steps Taken

1. Added a dedicated standalone driver for boundary-normal experiments.
2. Verified the driver on m = -1 first.
3. Repeated the same matrix for m = 3.
4. Added a projected-normal variant as a third normal definition.
5. Ran targeted m = 3 PETSc-normal retests with:
   • higher penalty
   • tighter Stokes tolerance

Main Comparison Settings

• uw_cellsize = 1/8
• Q2/Q1
• continuous pressure
• uw_p_bc = False
• default pressure gauge: volume mean
• serial solve with the standalone driver
• baseline natural penalty: 1e8
• baseline Stokes tolerance: 1e-5

Results

m = -1

case                    p_vol_l2(ana)   p_inner_l2(ana)   p_inner_l2(base)   inner_tan_rms
essential               0.870782        2.00023           -                  4.1907
natural_full            0.870904        1.99414           0.023663           4.19116
natural_normal_petsc    0.903829        2.38951           2.23654            3.99403
natural_normal_analytic 1.27216         15.2982           12.985             0.144261
natural_normal_projected
                        1.27379         15.6381           13.2424            0.213833

m = 3

case                    p_vol_l2(ana)   p_inner_l2(ana)   p_inner_l2(base)   inner_tan_rms
essential               0.893876        1.83758           -                  15.3523
natural_full            0.894121        1.83568           0.00880887         15.3557
natural_normal_petsc    0.910169        3.06533           0.805627           14.3317
natural_normal_analytic 0.895693        10.4955           3.60611            0.659744
natural_normal_projected
                        0.89858         11.1348           3.87889            1.42643

Tuning Checks

Targeted m = 3 PETSc-normal retests:

Higher penalty

• uw_vel_penalty = 1e10
• uw_stokes_tol = 1e-5

natural_normal_petsc: p_inner_l2(base) = 0.805639
baseline PETSc-normal: p_inner_l2(base) = 0.805627

Tighter tolerance

• uw_vel_penalty = 1e8
• uw_stokes_tol = 1e-8

natural_normal_petsc: p_inner_l2(base) = 0.805627
baseline PETSc-normal: p_inner_l2(base) = 0.805627

Neither change materially improved the PETSc-normal result.

Positive Results

• The standalone driver isolates the experiment cleanly from the production
  benchmark path.
• natural_full matches the essential baseline very closely in both m = -1
  and m = 3.
• All three normal-only variants were run successfully: PETSc, analytical, and
  projected.

Negative Results

• Normal-only weak BCs do not improve inner-boundary pressure relative to the
  essential baseline.
• Analytical normals are substantially worse than PETSc normals for this
  benchmark objective.
• Projected normals track the analytical-normal behavior more than the
  PETSc-normal behavior.
• Increasing the penalty or tightening the solve tolerance does not recover the
  essential baseline for the PETSc-normal normal-only case.

Key Observation

The important distinction is not just PETSc normal versus analytical normal. The bigger change is full-vector benchmark BC versus normal-only BC.

natural_full still tries to match the full analytical boundary velocity, so it
stays very close to essential.

natural_normal_* changes the boundary-condition class. It only penalizes the
normal component, so it behaves like a free-slip style condition rather than
the original benchmark BC.

That is visible in the inner-boundary tangential RMS:

• natural_full stays almost identical to essential
• natural_normal_analytic suppresses tangential motion very strongly
• natural_normal_projected behaves similarly
• natural_normal_petsc leaks more tangential motion because the facet normals
  are imperfect

That leakage makes PETSc normals look less bad relative to the essential
baseline, but it is not evidence that PETSc normals are a more accurate version
of the benchmark BC.

Conclusion

Replacing PETSc facet normals with analytical or projected normals in a
natural normal-only BC does not reduce inner-boundary pressure error for
this spherical Thieulot benchmark.

If the goal is to match the benchmark / essential solution, the best choices
remain:

• essential
• natural_full

If the goal is a genuine free-slip style experiment, natural_normal_* is a
valid separate test, but it should be treated as a different boundary-condition
problem rather than as a pressure fix for the benchmark.

[gthyagi]

Recent pressure-anchor debugging results

Test titles:

• Test 1: Normal-vector boundary-condition comparison (essential, natural_full, natural_normal_petsc, natural_normal_analytic, natural_normal_projected)
• Test 2: Pressure-anchor comparison (none, lower, upper, both) with reporting gauges (none, volume_mean, inner_surface_mean, outer_surface_mean)
• Test 3: Serial direct solve check (np=1, m=-1, m=3, 1/8)
• Test 4: MPI iterative solve check (np=2, m=-1, m=3, 1/8)
• Test 5: Underworld feature inspection for pressure anchoring, nullspace handling, and boundary-mean diagnostics

Spherical Thieulot Pressure-Anchor Debugging

Goal

Investigate why pressure remains inaccurate at the inner spherical boundary, even after the m = 3 body-force bug was fixed, and determine which current Underworld features are viable for improving the pressure field.

Relevant drivers:

• benchmarks/spherical/ex_stokes_thieulot.py
• benchmarks/spherical/ex_stokes_thieulot_bc_normals.py
• benchmarks/spherical/ex_stokes_thieulot_pressure_anchor.py

Scope

The earlier m = 3 failure is resolved. The remaining issue is the inner-boundary pressure mismatch.

The working hypothesis tested here was:

1. pressure is highly sensitive near the inner boundary because of strong radial scaling
2. the global constant-pressure nullspace may not control the inner-boundary pressure trace well enough
3. explicit pressure anchoring or alternative gauge choices might reduce the error

Tests Run

1. Normal-vector boundary-condition tests

These were already completed in the separate normal-BC driver.

Cases compared:

• essential
• natural_full
• natural_normal_petsc
• natural_normal_analytic
• natural_normal_projected

Main result:

• natural_full stayed close to essential
• natural_normal_* did not improve the inner-boundary pressure
• analytical or projected normals were worse than PETSc normals in the normal-only BC

Conclusion:

• boundary-normal accuracy is not the dominant cause of the pressure error in the current benchmark setup

2. Pressure-anchor experiments

Implemented in ex_stokes_thieulot_pressure_anchor.py.

Pressure modes:

• none: PETSc constant-pressure nullspace, no pressure Dirichlet BC
• lower: pressure Dirichlet on Lower
• upper: pressure Dirichlet on Upper
• both: pressure Dirichlet on both Lower and Upper

Reporting-only gauges:

• none
• volume_mean
• inner_surface_mean
• outer_surface_mean

3. Underworld feature inspection

Checked what current UW exposes for pressure control.

Available now:

• named-surface pressure Dirichlet BCs
• PETSc constant-pressure nullspace
• volume and boundary integrals for post-solve diagnostics

Not first-class:

• boundary mean pressure constraints
• global pressure mean constraints as solver constraints
• point pressure pins

Results

Serial direct solve, np = 1, 1/8, m = -1

Summary file:

• pressure_anchor_stage m=-1 summary

Key result:

• none, lower, upper, and both were effectively identical
• gauge shifts changed p_inner_l2 only in the fourth decimal place

Representative numbers:

mode	gauge	v_vol_l2	p_vol_l2	p_inner_l2
none	volume_mean	2.29948e-2	8.70782e-1	2.00023
lower	volume_mean	2.29948e-2	8.70782e-1	2.00023
upper	volume_mean	2.29948e-2	8.70782e-1	2.00023
both	volume_mean	2.29948e-2	8.70782e-1	2.00023

Serial direct solve, np = 1, 1/8, m = 3

Summary file:

• pressure_anchor_stage m=3 summary

Key result:

• same pattern as m = -1
• pressure anchors did not materially change the solved pressure
• gauge shifts changed p_inner_l2 only slightly

Representative numbers:

mode	gauge	v_vol_l2	p_vol_l2	p_inner_l2
none	volume_mean	7.17376e-2	8.93876e-1	1.83758
lower	volume_mean	7.17376e-2	8.93876e-1	1.83759
upper	volume_mean	7.17376e-2	8.93876e-1	1.83759
both	volume_mean	7.17376e-2	8.93876e-1	1.83759

MPI iterative solve, np = 2, 1/8, m = -1

Summary file:

• pressure_anchor_mpi2 m=-1 summary

Key result:

• none converged and reproduced the known good low-error solve
• both diverged in the linear solve with DIVERGED_PC_FAILED

Representative numbers:

mode	gauge	v_vol_l2	p_vol_l2	p_inner_l2	status
none	volume_mean	5.90108e-3	2.61469e-1	11.7999	converged
both	volume_mean	8.20630e-1	1.0	1.0	diverged

This means the pressure-anchored solve path is not robust in the current parallel fieldsplit/preconditioner setup.

MPI iterative solve, np = 2, 1/8, m = 3

Summary file:

• pressure_anchor_mpi2_clean m=3 summary

Key result:

• none converged cleanly
• a mixed none,both run previously stalled on the anchored case and had to be killed

Representative none numbers:

mode	gauge	v_vol_l2	p_vol_l2	p_inner_l2
none	volume_mean	3.99436e-2	3.59465e-1	20.5627

What Worked

• fixing the earlier m = 3 body-force bug
• keeping bc_type = essential or natural_full for the actual benchmark comparison
• using PETSc pressure nullspace with post-solve diagnostics
• using BdIntegral-based diagnostics to inspect inner and outer pressure means separately

What Did Not Work

• replacing PETSc normals with analytical normals in a normal-only BC
• projected normals in the same normal-only BC
• using normal-only weak BCs as a way to improve pressure
• relying on a different constant gauge to fix the inner-boundary pressure mismatch
• using hard pressure anchors as a general fix
  • in serial/direct solves they barely changed the solution
  • in parallel/iterative solves they destabilized the solve path

Root Cause

The remaining inner-boundary pressure issue is not primarily a normal-vector problem and not primarily a constant pressure-gauge problem.

The evidence points to a boundary-local pressure-trace sensitivity near the inner spherical surface, where pressure has strong radial dependence and small algebraic errors are amplified. The global constant-pressure nullspace removes only the volume-wide additive mode; it does not provide a surgical control on the inner-boundary pressure trace. Hard pressure Dirichlet anchors change the algebraic problem and can alter or destabilize the iterative Stokes solve, but they do not represent a clean fix for the actual inner-boundary pressure shape.

In short:

• the pressure error is localized and boundary-sensitive
• the current nullspace/gauge treatment is too coarse to target it directly
• the available hard-anchor workaround is not robust in the iterative parallel path

Practical Recommendations Using Current UW Features

Best current approach for benchmark-quality solves

• keep bc_type = essential
• keep petsc_use_pressure_nullspace = True when there is no physical pressure BC
• do not use normal-only BCs as a pressure fix
• do not treat p_bc = true as a general solution improvement unless the physical model really has a pressure BC

Best current approach for debugging

• compare inner and outer pressure diagnostics separately using BdIntegral
• monitor:
  • volume pressure mean
  • inner-surface pressure mean
  • outer-surface pressure mean
  • inner/outer boundary relative pressure errors
• use reduced serial or small-MPI cases to separate algebraic issues from formulation issues

Cautious use of pressure anchoring

• named-surface pressure Dirichlet BCs are available now
• but they should only be used if the physics genuinely prescribes pressure there
• for this benchmark, they are better treated as a diagnostic perturbation than as a fix

Worthwhile UW Improvements

These would address the problem much more directly than current hard Dirichlet anchors.

1. First-class boundary-mean pressure constraints

   • example: enforce ∫_Lower p dS = 0
   • this is the most surgical control for the present problem

2. First-class point pressure pins

   • example: fix pressure at a single inner-boundary point
   • useful for gauge control without overconstraining the whole boundary

3. Better support for pressure constraints in parallel fieldsplit/preconditioned Stokes solves

   • current hard pressure anchors can destabilize the parallel solve path

4. Better boundary pressure / traction recovery utilities on curved spherical boundaries

   • especially for applications where tractions and surface/basal stresses matter

5. More explicit nullspace / gauge-control utilities in Python

   • boundary-local gauge tools are missing at present

Implications for Large-Scale Applications

This matters for global mantle convection, spherical-cap models, and global subduction models because those models also depend on accurate pressure or traction interpretation near curved shell boundaries.

Even when the velocity field looks acceptable:

• boundary-local pressure errors can remain significant
• global zero-mean pressure is not enough to control the boundary trace
• hard pressure anchors may change solver behavior rather than actually fixing the pressure field

For larger production models, the safest current practice is:

• preserve a physically consistent velocity BC treatment
• use PETSc nullspace when pressure is gauge-free
• diagnose boundary pressure separately from global pressure
• treat hard pressure BCs as model changes, not as benchmark-neutral fixes

[gthyagi]

Spherical Thieulot Radial Grading

Driver:

• ex_stokes_thieulot_radial_grading.py

Setup:

• radial mesh deformation only; benchmark fields and PDE are unchanged
• grading_ratio = h_outer / h_inner
• grading_ratio = 4 gives inner_scale ≈ 0.462

Short Summary

• 1/8, serial/direct:
  • m = -1: moderate inner grading helped; best tested grading_ratio ≈ 4
  • m = 3: inner grading did not help
• 1/8, np = 8:
  • m = 3: inner grading helped strongly; grading_ratio = 4 was the best tested compromise
• 1/16, np = 8, after the UW BdIntegral fix:
  • full inner/outer boundary diagnostics now complete for graded meshes
  • grading improves inner-boundary pressure but worsens outer-boundary pressure

Latest 1/16, np = 8 Results

Using the volume_mean gauge:

case	grading ratio	inner scale	velocity L2	pressure volume L2	pressure inner L2	pressure outer L2
m=-1	1	1.000	7.71041e-4	5.96441e-2	14.4094	3.88626
m=-1	4	0.462	4.65460e-4	2.41538e-2	5.51350	13.1031
m=3	1	1.000	5.57220e-3	4.05191e-2	16.4700	9.00753
m=3	4	0.462	1.29199e-3	2.45953e-2	4.56341	13.0849

Interpretation:

• m=-1: inner pressure 14.41 -> 5.51, outer pressure 3.89 -> 13.10
• m=3: inner pressure 16.47 -> 4.56, outer pressure 9.01 -> 13.08
• in both cases, grading helps the inner shell and global pressure norm, but pushes error toward the outer shell

Files

Recent full-diagnostic summaries:

• m=-1 grade 1
• m=-1 grade 4
• m=3 grade 1
• m=3 grade 4

Note: based on these results, it is also worth implementing curved elements (degree-2 mesh elements) in Underworld for spherical-shell benchmarks and related applications.