From 57d156c69009a3e1209d42bf5fd251eda5cac187 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Apr 2026 10:34:35 +1000 Subject: [PATCH 1/2] Fix bare-variable composition under sympy Matrix dispatch (#137) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit MathematicalMixin advertises (in CLAUDE.md and MATHEMATICAL_MIXIN_DESIGN.md) that mesh and swarm variables can be used directly in sympy arithmetic without explicit .sym access. In practice it only worked when the bare variable was the *innermost* operand. As soon as a bare variable appeared on the right of an already-sympified subexpression — for example as the second factor inside a sympy.exp() of a product — composition raised: TypeError: Incompatible classes , Cause ----- .sym on a scalar variable returns a 1×1 sympy Matrix. SymPy's Matrix.__mul__ raises TypeError directly instead of returning NotImplemented, so Python's normal fall-through to the right operand's __rmul__ never fires. The mixin already had a working __rmul__ — it just was never being called. Fix --- SymPy provides _op_priority as the documented opt-in escape hatch: any class with _op_priority strictly greater than the LHS's wins dispatch. Matrix._op_priority is 10.01; Symbol/Mul/Expr are 10.0. Setting MathematicalMixin._op_priority = 11.5 makes sympy delegate "Matrix * " and similar mixed-form expressions to our existing reverse dunders, which sympify self via .sym and re-do the operation cleanly as Matrix * Matrix. This is a one-line class attribute change. The reverse dunders (__rmul__, __radd__, __rsub__, __rtruediv__, __rpow__) were already in place and correct; they just needed to be reachable. Verification ------------ - The four-case reproducer from issue #137 (A_sym_both, B_bareC_symT, C_symC_bareT, D_bare_both): pre-fix C and D raise TypeError, post-fix all four return ImmutableDenseMatrix and simplify to identical expressions. - tests/test_0726_bare_variable_composition_137.py covers all four forms plus an equivalence assertion. 5/5 pass in 2.3s. - pytest -m "level_1 and tier_a" on amr-dev: 56 passed, 3 skipped, 0 failed. No regressions. Closes #137. Underworld development team with AI support from Claude Code --- .../utilities/mathematical_mixin.py | 13 +++ ...test_0726_bare_variable_composition_137.py | 95 +++++++++++++++++++ 2 files changed, 108 insertions(+) create mode 100644 tests/test_0726_bare_variable_composition_137.py diff --git a/src/underworld3/utilities/mathematical_mixin.py b/src/underworld3/utilities/mathematical_mixin.py index 31f572b5b..2bb87ba26 100644 --- a/src/underworld3/utilities/mathematical_mixin.py +++ b/src/underworld3/utilities/mathematical_mixin.py @@ -26,6 +26,19 @@ class MathematicalMixin: - Consistent behavior across different variable types """ + # BUGFIX(#137): SymPy's binary operators on Matrix raise TypeError directly + # instead of returning NotImplemented, so Python's normal fall-through to + # the right operand's __rmul__/__radd__/etc. never fires. SymPy provides an + # opt-in escape hatch: any class with _op_priority strictly greater than + # the LHS's wins dispatch (Matrix._op_priority = 10.01, Symbol/Expr = 10.0). + # Setting this above all sympy core priorities makes sympy delegate + # `Matrix * ` and similar mixed-form expressions to our + # reverse dunders, which then sympify self via .sym and re-do the + # operation cleanly. This is what makes the "no .sym needed" promise in + # CLAUDE.md and MATHEMATICAL_MIXIN_DESIGN.md actually hold for the + # bare-variable-on-the-right composition case (issue #137 cases C and D). + _op_priority = 11.5 + def _validate_sym(self): """Validate that sym property is available and valid.""" try: diff --git a/tests/test_0726_bare_variable_composition_137.py b/tests/test_0726_bare_variable_composition_137.py new file mode 100644 index 000000000..fae3d4388 --- /dev/null +++ b/tests/test_0726_bare_variable_composition_137.py @@ -0,0 +1,95 @@ +"""Regression test for issue #137 — bare-variable composition asymmetry. + +MathematicalMixin advertises that mesh / swarm variables can be used directly +in sympy arithmetic without explicit ``.sym`` access. Pre-fix this only worked +when the bare variable was the *innermost* operand; as soon as a bare variable +appeared on the right of an already-sympified subexpression (e.g. inside a +sympy ``exp()`` of a product), composition raised: + + TypeError: Incompatible classes + , + +Cause: ``.sym`` on a scalar returns a 1×1 sympy ``Matrix``. SymPy's +``Matrix.__mul__`` raises TypeError directly instead of returning +NotImplemented, so Python's normal fall-through to the right operand's +``__rmul__`` never fires. + +Fix: ``MathematicalMixin._op_priority = 11.5`` (above sympy's +``Matrix._op_priority = 10.01``) makes sympy delegate the operation to our +reverse dunder, which then sympifies ``self`` via ``.sym`` and re-runs the +multiplication cleanly as ``Matrix * Matrix``. +""" + +import pytest +import sympy + +import underworld3 as uw + + +pytestmark = pytest.mark.level_1 + + +@pytest.fixture(scope="module") +def vars_TC(): + """Two scalar variables of different kinds: a MeshVariable and a SwarmVariable.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5, + ) + T = uw.discretisation.MeshVariable("T_137", mesh, 1, degree=1) + + swarm = uw.swarm.Swarm(mesh) + C = uw.swarm.SwarmVariable("C_137", swarm, size=1, proxy_degree=1) + swarm.populate(fill_param=2) + return T, C + + +@pytest.mark.parametrize( + "label,build", + [ + ("A_sym_both", lambda T, C: sympy.exp(-C.sym * T.sym)), + ("B_bareC_symT", lambda T, C: sympy.exp(-C * T.sym)), + ("C_symC_bareT", lambda T, C: sympy.exp(-C.sym * T )), + ("D_bare_both", lambda T, C: sympy.exp(-C * T )), + ], +) +def test_bare_variable_composition_under_sympy_function(vars_TC, label, build): + """All four mixed-form combinations must compose cleanly under sympy.exp. + + Pre-fix, cases C and D raised TypeError. Post-fix all four return a + sympy expression of identical structure. + """ + T, C = vars_TC + eta_0 = sympy.symbols("eta_0_137") + result = eta_0 * build(T, C) + # Result should be a sympy object (Matrix or Expr) — type may differ + # between cases but the structure is equivalent. + assert result is not None + # The four results should all simplify to the same canonical form. + # We can't easily compare across the parametrize boundary in a + # parametrised test, but we can at least check it's sympifiable. + assert hasattr(result, "free_symbols") or hasattr(result, "shape") + + +def test_bare_variable_composition_all_forms_agree(vars_TC): + """The four mixed-form expressions should be mathematically equivalent.""" + T, C = vars_TC + eta_0 = sympy.symbols("eta_0_137") + + forms = [ + eta_0 * sympy.exp(-C.sym * T.sym), + eta_0 * sympy.exp(-C * T.sym), + eta_0 * sympy.exp(-C.sym * T ), + eta_0 * sympy.exp(-C * T ), + ] + + # Reduce to a comparable scalar by extracting the [0,0] element if needed. + def scalarise(x): + if hasattr(x, "shape") and x.shape == (1, 1): + return x[0, 0] + return x + + canonical = [sympy.simplify(scalarise(f) - scalarise(forms[0])) for f in forms] + for i, diff in enumerate(canonical): + assert diff == 0, ( + f"form {i} differs from form 0 after simplification: residual={diff}" + ) From aba85fd7d42e9a7e7e4e72260c8a923afb3eabdb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 1 May 2026 13:32:12 +1000 Subject: [PATCH 2/2] Pin UWexpression _op_priority back to sympy default (avoid Matrix/UWexpression regression) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous commit set MathematicalMixin._op_priority = 11.5 to win sympy dispatch over Matrix (10.01) for the bare-variable composition case. UWexpression inherits from MathematicalMixin (first in MRO), so that change inadvertently bumped UWexpression's priority too — and UWexpression has its OWN __rtruediv__ (and __rmul__, etc.) for sympy interop. Its __rmul__ explicitly handles MutableDenseMatrix; its __rtruediv__ does not, falling through to Symbol.__rtruediv__ which fails on MatrixBase. Result: solvers that compute Matrix / UWexpression (e.g. SNES_Diffusion's F0 expression `self.DuDt.bdf(0) / self.delta_t` at solvers.py:2499) broke with TypeError on every transient/Darcy/AdvDiff path. CI surfaced this on test_1005, test_1006, test_1100, test_1110. Fix: explicitly set UWexpression._op_priority = 10.0 to opt out of the mixin's bump. UWexpression already handles Matrix dispatch correctly through its own dunders (and via sympy's standard machinery for fall-through cases); it doesn't need the priority bump that pure MathematicalMixin classes (EnhancedMeshVariable, SwarmVariable) need. Verified: - #137 4-case reproducer still passes (the original fix is intact). - Previously-failing tests all pass: test_1005_TransientDarcyCartesian, test_1006_RichardsCartesian, test_1100_AdvDiffCartesian, test_1110_advDiffAnnulus — 8/8 pass. - pytest -m "level_1 and tier_a" on amr-dev: 56 passed, 3 skipped, 0 failed. Underworld development team with AI support from Claude Code --- src/underworld3/function/expressions.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/underworld3/function/expressions.py b/src/underworld3/function/expressions.py index d6e6f8cb8..2f24914b3 100644 --- a/src/underworld3/function/expressions.py +++ b/src/underworld3/function/expressions.py @@ -629,6 +629,18 @@ class UWexpression(MathematicalMixin, uw_object, Symbol): # Slot for unique ID used in _hashable_content (like sympy.Dummy) __slots__ = ('_uw_id',) + # Override the MathematicalMixin priority bump back to sympy's default. + # MathematicalMixin sets _op_priority = 11.5 to win dispatch over + # sympy.Matrix (10.01) for the bare-variable composition case (#137 — + # MeshVariable / SwarmVariable on the right of a sympified subexpression). + # UWexpression is itself a sympy.Symbol subclass with its own __rmul__ / + # __rtruediv__ that already handle the Matrix case; inheriting the high + # priority would route Matrix / UWexpression through UWexpression's + # __rtruediv__ (which falls back to Symbol.__rtruediv__ → fails on + # MutableDenseMatrix). Pin it back to 10.0 so sympy's standard + # Matrix-dispatch path keeps handling these. + _op_priority = 10.0 + def __new__( cls, name,