Request
Expose the inverse of the transversely isotropic viscosity tensor (the compliance) alongside flux, so strain rate can be inferred from a recovered stress rather than differentiated from velocity.
Why it is nearly free
The Mühlhaus–Moresi tensor decomposes as
C = 2 eta_0 I - 2 (eta_0 - eta_1) P
with P the projector onto the on-plane shear modes (the bracketed n-terms in TransverseIsotropicFlowModel). C therefore has two eigenvalues — 2 eta_1 on that subspace, 2 eta_0 off it — and the inverse is the same decomposition with those reciprocated:
C^-1 = (1 / 2 eta_0) I + (1 / 2 eta_1 - 1 / 2 eta_0) P
Same structure, different scale. So _build_c_tensor can build the compliance unchanged under the substitution
eta_0 -> 1 / (4 eta_0)
eta_1 -> 1 / (4 eta_1)
(the quarter absorbing the factor-of-2 convention in both directions).
Why it is wanted
- Plasticity: the return map is derived in compliance form and implemented in stiffness form; the inverse is what maps between them.
- Diagnostics: traction (
sigma . n) is continuous across a fault-zone boundary while the strain rate is not, so the sound route is to recover stress by projection and infer strain rate through the compliance. Differentiating the velocity instead gives a field that jumps at the zone edge and recovers badly. This came up while building stress diagnostics for a TI fault zone (~/+Simulations/listric_extension/fuse_stress_p1.py); the TI case is the one where it has not been done.
Check before trusting it
The substitution assumes the bracketed term is exactly 2P with P idempotent in the implemented form. Build C and C^-1 numerically for several directors and assert C : C^-1 == I to round-off. That either confirms the substitution or exposes a convention factor — much cheaper to settle now than inside a return map.
Underworld development team with AI support from Claude Code
Request
Expose the inverse of the transversely isotropic viscosity tensor (the compliance) alongside
flux, so strain rate can be inferred from a recovered stress rather than differentiated from velocity.Why it is nearly free
The Mühlhaus–Moresi tensor decomposes as
with
Pthe projector onto the on-plane shear modes (the bracketed n-terms inTransverseIsotropicFlowModel).Ctherefore has two eigenvalues —2 eta_1on that subspace,2 eta_0off it — and the inverse is the same decomposition with those reciprocated:Same structure, different scale. So
_build_c_tensorcan build the compliance unchanged under the substitution(the quarter absorbing the factor-of-2 convention in both directions).
Why it is wanted
sigma . n) is continuous across a fault-zone boundary while the strain rate is not, so the sound route is to recover stress by projection and infer strain rate through the compliance. Differentiating the velocity instead gives a field that jumps at the zone edge and recovers badly. This came up while building stress diagnostics for a TI fault zone (~/+Simulations/listric_extension/fuse_stress_p1.py); the TI case is the one where it has not been done.Check before trusting it
The substitution assumes the bracketed term is exactly
2PwithPidempotent in the implemented form. BuildCandC^-1numerically for several directors and assertC : C^-1 == Ito round-off. That either confirms the substitution or exposes a convention factor — much cheaper to settle now than inside a return map.Underworld development team with AI support from Claude Code