Request alternative free-slip BC implementations for shell geometries

[gthyagi]

Summary

The current penalty-based free-slip boundary condition path looks too fragile for spherical shell models and would benefit from alternative implementations.

In annulus benchmarks, the penalty method remains stable even for very large velocity penalties, although it can become very expensive. In spherical benchmarks, the same approach appears much more sensitive: free-slip runs can become unstable / memory-heavy at moderate resolution, and the solution quality is strongly penalty-dependent.

This suggests a scalability and/or formulation issue that is much more pronounced in spherical geometries than in annulus geometries.

Reproduction

Benchmark script:

• https://github.com/gthyagi/UW3_Annulus_Spherical_Benchmarks/blob/main/benchmarks/spherical/ex_stokes_kramer.py

Example spherical free-slip run that became problematic:

mpirun -np 8 python ex_stokes_kramer.py -uw_cellsize 1/16 -uw_case case2 -uw_vel_penalty 1e2

At coarser resolution (1/8), the same spherical free-slip case shows strong penalty sensitivity:

• vel_penalty=1e2: solves quickly, but accuracy degrades (v_l2 ~ 2.88e-1, p_l2 ~ 3.83e-2)
• vel_penalty=1e4: much better accuracy (v_l2 ~ 5.15e-2, p_l2 ~ 1.69e-2)
• vel_penalty=1e7: run became impractically slow / stalled

So for spherical shells, the penalty is doing double duty:

• too small: free-slip is under-enforced and errors increase
• too large: conditioning and runtime become problematic

Why this seems geometry-specific

In the annulus Kramer benchmarks, the penalty approach is much more robust. Even with a very large penalty (2.5e8), the free-slip run remains stable, although it is substantially more expensive than the no-slip run.

From the annulus timing comparison at the same resolution / MPI size:

• free-slip: case2_inv_lc_32_..._bc_paper
• no-slip: case4_inv_lc_32_..._bc_paper

Key PETSc timings:

• SNES_Stokes_SaddlePt.solve: 3.3674e+02 s (free-slip) vs 7.1149e+00 s (no-slip)
• MatMult: 3.2254e+02 s vs 1.6257e+00 s
• PCApply: 1.9113e+02 s vs 1.0640e+00 s
• PCSetUpOnBlocks: 1.3967e+02 s vs 7.3120e-01 s
• KSPSolve_FS_Schu: 1.6947e+02 s vs 9.5761e-01 s
• VecScatterEnd: 1.1411e+02 s vs 5.9967e-01 s

So the annulus penalty method is already expensive, but still usable/stable. In spherical models, the same general approach appears much less robust.

Request

It would be useful to add alternative free-slip boundary condition implementations to Underworld, especially for curved shell geometries.

Possible directions:

• strong enforcement in a local normal/tangential basis (equivalently: enforce zero normal velocity strongly while leaving tangential components free)
• Lagrange multiplier approaches for the normal constraint
• Nitsche-type slip enforcement or other weak methods that are less penalty-sensitive than a pure penalty term

Comparison with the benchmark paper

The original Kramer benchmark paper does not describe a pure penalty free-slip treatment for the Fluidity reference implementation.

Instead, it states that:

• strong Dirichlet boundary conditions are used
• for free-slip, the boundary velocity basis is locally rotated into normal/tangential components
• the normal component is constrained to zero strongly
• tangential components remain free
• rigid-body rotation modes are then projected out before computing errors

That is a materially different approach from the current penalty-based workaround, and it seems worth considering a UW implementation that is closer to that benchmark treatment for shell geometries.

Why this matters

This is relevant beyond the benchmark itself. If penalty-based free-slip is highly geometry- and penalty-sensitive in spherical shells, that can affect larger-scale spherical applications where free-slip is a natural boundary choice.

A more robust shell free-slip implementation would likely improve:

• reproducibility across resolutions
• solver conditioning / runtime
• pressure and velocity convergence behavior in curved geometries

[gthyagi]

Additional spherical free-slip sensitivity data for case2 (np=8, uw_cellsize=1/8):

• vel_penalty=1e4
  • v_l2 = 5.14773842247616e-02
  • p_l2 = 1.0866899839169922e-02
• vel_penalty=1e5
  • v_l2 = 1.8482878125726725e-02
  • p_l2 = 1.092159734140487e-02

The main point is that the velocity error is very sensitive to the penalty value, even when the pressure error changes only slightly. That is another sign that the current pure-penalty free-slip treatment is not a clean formulation for convergence studies in spherical shells: changing the penalty materially changes the discrete solution, not just the solver cost.

[lmoresi]

The local rotation formulation is tricky given the way PETSc is formulated. But the Nitsche approach is widely used. I don't have capacity to work on that but I can review the code you come up with. Please keep this on a feature branch as I expect it will create a lot of solver changes that we will need to manage carefully.

[lmoresi]

For the Lagrange multiplier approach, did you have in mind adding an additional term to the constraint equation ? The advantage of this is that it would effectively be a topography term obtained in a similar way to the pressure. Not sure about the implementation given the surfaces are not present on every node.

[gthyagi]

I’ll move the local rotation formulation to later this quarter or the next. I agree that the Nitsche approach is the most practical next step with minimal code changes. We can revisit the Lagrange multiplier approach in the future if we encounter a problem where it’s clearly needed.

For now, I’ll proceed with the penalty method and focus on getting the benchmarks running, so we can generate results and build some user traction for the spherical and annulus models.

[lmoresi]

Infrastructure for Nitsche BCs now in place

We've investigated the feasibility of Nitsche free-slip and confirmed it fits within PETSc's existing boundary callback API — no PETSc modifications needed. The key references are PETSc's own ex62.c (Nitsche Stokes) and Sime & Wilson (2020), "Automatic weak imposition of free slip boundary conditions via Nitsche's method" (arXiv:2001.10639).

What's been done

Commit 3b7ab6e on development enables the full boundary weak form infrastructure in the Stokes solver:

• f1_bd (gradient boundary residual): previously always NULL, now wired through when a natural BC provides fn_F. This is where the Nitsche symmetry term lives.
• g2_bd, g3_bd (boundary Jacobians): previously NULL, now computed via derive_by_array and passed to PetscDSSetBdJacobian when f1_bd is present.
• Pressure-velocity (pu) and pressure-pressure (pp) boundary Jacobian blocks: previously commented out, now enabled. Currently zeros — Nitsche will provide the u.n / n coupling terms.
• NaturalBC namedtuple gains fn_F field for the gradient residual.

No behavioral change for existing natural BCs (fn_F=None → same NULL pattern as before). All tests pass.

What remains for Nitsche free-slip

1. Symbolic construction: Build the Nitsche weak form terms (f0_bd, f1_bd, pressure coupling) from mesh.Gamma_N (normals) and the constitutive model viscosity. Auto-differentiation handles the Jacobians.
2. Mesh size h: Pass a per-cell size estimate to boundary functions (DG-0 auxiliary field with cell inradius).
3. User API: stokes.add_nitsche_freeslip_bc("Upper", gamma=10, theta=-1) — symmetric (theta=1) or skew-symmetric (theta=-1, unconditionally stable).

The skew-symmetric variant (theta=-1) works with UW3's existing FGMRES + fieldsplit solver stack — no solver changes needed.

Plan file: ~/.claude/plans/nitsche-freeslip-implementation.md

[lmoresi]

OK - turns out the Nitsche BC is not so tricky to implement once all of the infrastructure part is in place.

Main thing to note: symmetric version does behave better with current solver default config, so that's where we will leave it. There are some tests to show this.

[lmoresi]

Nitsche free-slip BC implementation — PR #106

@gthyagi — PR #106 adds stokes.add_nitsche_bc() as an alternative to the penalty free-slip. It implements the full Nitsche weak form:

• Penalty/stabilisation: (gamma * mu / h) * (u.n) * n
• Consistency: stress flux from the constitutive model projected onto the normal
• Symmetry: adjoint consistency term (theta parameter controls variant)
• Pressure coupling: p * n on velocity equation, u.n on pressure equation

Usage

# Replace this:
stokes.add_natural_bc(1e6 * Gamma.dot(v.sym) * Gamma, "Upper")

# With this:
stokes.add_nitsche_bc("Upper", gamma=10)

Validated on Cartesian box

Compared essential BC (strong free-slip), penalty (1e4), and Nitsche (gamma=10, theta=+1):

• Nitsche velocity is 0.08% from essential BC solution (penalty: 0.15%)
• Normal velocity constraint: Nitsche 2.3e-5, penalty 2.7e-5
• Solver cost: 1 Newton iteration, same as penalty

Next steps for spherical shell testing

The PR needs testing on the spherical shell cases from this issue. Could you try running your Kramer benchmark script with add_nitsche_bc("Upper", gamma=10) replacing the penalty free-slip on the outer boundary? The key question is whether Nitsche remains stable at moderate resolution where penalty became problematic.

The gamma parameter should be much less sensitive than the penalty — theory says gamma ~ 10 for P2 elements regardless of mesh size or viscosity.