TI-VEP viscosity reports the yield-limited weak-plane value; Nitsche scales on K (#463)

[lmoresi]

TI-VEP .viscosity now reports the yield-limited fault-plane viscosity the stress actually uses

Fixes#463

The bug

TransverseIsotropicVEPFlowModel.viscosity computed the yield-limited
fault-plane viscosity eta_1_eff and then returned the un-yielded bulk
shear_viscosity_0 — the computed value was a dead local. The residual was
never affected (the flux path builds the rank-4 tensor through
_eta_for_tensor, which applies the yield correctly), but every reporting
surface that reads .viscosity — diagnostics, viscosity renders, projections —
saw no yielding at all.

Measured on a yielding shear state (tau_y = 0.5, fault-plane shear rate 50,
eta_0 = 1, eta_1 = 0.01): projecting .viscosity gave a constant 1.0
(bulk) while the tensor coefficient the stress uses projects to 0.005
(the yield limit tau_y / (2 |gamma_dot|)).

The fix

• viscosity now delegates to _eta_for_tensor with the ACTIVE integrator
  mode, so the property returns exactly the coefficient the stress tensor is
  built from. Under yield that is the persistent _eta1_yield_eff container
  introduced in Stokes uu_G3 transposed tangent + TI un-frozen Picard: fix issue #457 at source #493 — a wrapped UWexpression atom, so the Picard freezing
  contract holds for anything that bakes the property. For the hybrid
  integrator it reports the BDF (yield-clipped) branch, matching the default
  self._c.
• TransverseIsotropicVEPFlowModel.plastic_fraction had the same drift
  squared: it compared the bulk eta_0 against eta_1_ve (identically zero
  after the Max(0, ...)), and a misplaced ternary made the hasattr guard
  a no-op. It now reports 1 - eta_1_eff / eta_1_ve on the weak-plane
  channel (the bulk is structurally non-yieldable).
• K (the Schur/stiffness scale) stays shear_viscosity_0 — deliberately,
  per the issue.
• The two Nitsche boundary kernels (SNES_Vector.add_nitsche_bc,
  SNES_Stokes_SaddlePt.add_nitsche_bc) used .viscosity as their PENALTY
  scale. That consumer wants the stiffness scale, not the (now weak-plane,
  possibly yield-collapsed) effective viscosity — they now read .K, which
  is bit-identical to the old behaviour for every model (isotropic models:
  K == viscosity; TI models: K == eta_0 == the old .viscosity return).

The deferred .sym tidy from #493 (isotropic VEP)

PR #493 noted that ViscoElasticPlasticFlowModel._build_c_tensor still baked
the UNWRAPPED .sym contents of the effective viscosity into the c-tensor
(tangents were symmetric either way because the VEP solve flux goes through
stress()). That tidy is done here for contract uniformity:

• ViscoElasticPlasticFlowModel.viscosity now stores the combined
  (yield-limited, floored) coefficient inside a persistent
  _vep_eff_viscosity container and returns it — the
  ViscoPlasticFlowModel._plastic_eff_viscosity pattern. (Bonus: the
  property now always returns a UWexpression, so _object_viewer no longer
  breaks on a yielding VEP model.)
• The _build_c_tensor override is deleted; the base ViscousFlowModel
  build bakes the wrapped container atom element-wise (and the stray debug
  prints in the override go with it).

Proof nothing changes:

• The isotropic Picard tangent was ALREADY frozen before the tidy (the
  strain-rate content sat inside wrapped atoms one level down) — the new
  frozenness test test_vep_c_tensor_coefficients_are_frozen passes on both
  sides of the change, with a Newton-unwrap positive control.
• Solve identity: a 6-step shear-box loading sequence through yield
  (_yield_mode="min", tau_y = 0.5, the test_1052 configuration) was run on
  the pre-change and post-change builds for BOTH the isotropic VEP and the
  TI-VEP models. TI-VEP: v, p and the tau_xy trace are bitwise identical
  (its flux path is untouched). Isotropic VEP: identical to 3.2e-10
  max-abs (v: 9.4e-11, p: 2.3e-10, tau_xy: 3.2e-10) against a 1e-6 solver
  tolerance — not bitwise, because routing the coefficient through the
  container changes the compiled kernel's expression form, not its value.

Tests

tests/test_1066_stokes_jacobian_layout.py (alongside the #457/#493
freezing-contract tests):

• test_ti_vep_viscosity_property_reports_yield_limited_eta1 — the property
  must return the same container _eta_for_tensor bakes; frozen under
  d/dL, Newton unwrap reveals the yield law (positive control), and with
  the yield limit disabled it reports the unlimited VE fault-plane viscosity
  (negative control). Validated fail-before: on pre-fix code it fails at
  assert visc is not cm.Parameters.shear_viscosity_0.
• test_vep_c_tensor_coefficients_are_frozen — isotropic VEP pin for the
  tidy (passes pre- and post-change, as intended for a no-behaviour-change
  refactor).

Stakeholder suites: test_1066 (7 passed), test_0104 / test_0610 / test_1055 /
test_1057 / test_1059 / test_1060 / test_1065 / test_1067 (69 passed),
test_1052 VEP stability regression, level_2 (5 passed).

Gate: pytest tests -m "level_1 and tier_a" -q --ignore=tests/test_0050_utils.py
→ 592 passed, 0 failed, 17 skipped (MPI opt-in), 1 xfailed, in 7:09.

Underworld development team with AI support from Claude Code

[lmoresi]

Review

We verified the branch independently of the implementing agent:

• The defect is real and was measured 200×: pre-fix, cm.viscosity projected 1.0 while the flux coefficient projected 0.005 (τ_y/(2|γ̇|)) on an actively-yielding TI-VEP state. The property now delegates to _eta_for_tensor(...) — it returns the same wrapped container the tensor bakes, so reporting and stress cannot drift again and the freezing contract is intact (d/dL frozen assertion + Newton-unwrap positive control in the new test, fail-before validated).
• The Nitsche penalty-scale switch (.viscosity → .K) is the right call and is scaling-identical on every existing model (isotropic: K == viscosity; TI: K == η₀ == the old return) — verified K is defined on the base class with family overrides, so no model reaching the kernels lacks it. Post-TransverseIsotropicVEPFlowModel.viscosity discards the yield-limited eta_1 #463-semantics, .viscosity would have under-scaled the penalty by the yield ratio on TI-VEP; this change is protective, not behavioural.
• The Stokes uu_G3 transposed tangent + TI un-frozen Picard: fix issue #457 at source #493 deferred tidy (isotropic VEP wrapped-atom container, .sym-baking override deleted): frozenness test passes on both sides — the tangent was already frozen one atom-level down, so the tidy is structural; solve identity 3.2e-10 max diff vs 1e-6 tolerance, TI-VEP arm bitwise identical.
• plastic_fraction's double defect (η₀-vs-η₁ comparison + the no-op hasattr ternary) is fixed on the weak-plane channel.
• Gate 592 passed / 0 failed; eight constitutive/Nitsche/yield stakeholder files 69/69; test_1052 VEP stability 5/5.

Clean. Merging on CI green.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes a reporting/diagnostics bug in TI-VEP rheology where TransverseIsotropicVEPFlowModel.viscosity did not reflect the yield-limited weak-plane coefficient actually used in the stress tensor, and adjusts Nitsche penalty scaling to use the intended stiffness scale (K) rather than the (now yield-sensitive) .viscosity.

Changes:

• Update TransverseIsotropicVEPFlowModel.viscosity to return the yield-limited weak-plane coefficient produced by _eta_for_tensor(...), and fix plastic_fraction to report a meaningful weak-plane plastic fraction.
• Refactor isotropic ViscoElasticPlasticFlowModel.viscosity to return a persistent wrapped container (frozen Picard tangent contract) and drop its custom _build_c_tensor override in favor of the base implementation.
• Switch Nitsche BC penalty scaling in the generic SNES kernels from .viscosity to .K to preserve stiffness scaling for TI models.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated no comments.

File	Description
src/underworld3/constitutive_models.py	Fixes TI-VEP viscosity/plastic_fraction reporting drift; introduces persistent wrapped VEP effective-viscosity container and removes the unwrapped c-tensor override.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Uses constitutive-model stiffness scale K for Nitsche penalties to avoid under-scaling when TI .viscosity yields on the weak plane.
tests/test_1066_stokes_jacobian_layout.py	Adds regression tests for TI-VEP viscosity reporting and isotropic VEP c-tensor coefficient freezing contract.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.