TransverseIsotropicVEPFlowModel.viscosity discards the yield-limited eta_1

[lmoresi]

TransverseIsotropicVEPFlowModel.viscosity (src/underworld3/constitutive_models.py, ~3523-3548) computes the yield-limited fault-plane viscosity and then returns the un-yielded bulk value:

eta_1_eff = inner_self.ve_effective_viscosity
if self.is_viscoplastic:
    vp_eff = self._plastic_effective_viscosity
    eta_1_eff = self._combine_yield(eta_1_eff, vp_eff)
return inner_self.shear_viscosity_0        # <- eta_1_eff is never used

eta_1_eff is a dead local. The docstring says the property "Applies the yield mode (softmin/min/harmonic) to eta_1, leaving eta_0 (bulk) unchanged", which the return value does not do.

Impact

Reporting and derived quantities only. The residual is unaffected: the flux path builds the rank-4 tensor through _eta_for_tensor / _assemble_c_tensor, which apply the yield correctly. But anything reading .viscosity on a TI-VEP model — diagnostics, viscosity renders, generic solver code that asks a constitutive model for "the" viscosity — sees no yielding at all, so a yielding fault looks unyielded in every plot.

Suggested fix

Decide what the property means for an anisotropic model and make the code and docstring agree. Either return the yield-limited eta_1_eff (matching the docstring), or return eta_0 and document that the scalar .viscosity of a TI model is the bulk value with the yield living in the tensor — but then drop the dead computation.

Note K (the Schur preconditioner scaling, next property down) is separately and deliberately shear_viscosity_0; that one is correct and should not change.

Found while auditing the rotated free-slip preconditioner under transverse isotropy.

Underworld development team with AI support from Claude Code