From 31e7fd2f4add82a7ba961f3b303585407cf05f9b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 13 Jun 2026 07:58:59 +1000 Subject: [PATCH] projection: unit-aware smoothing_length with validation and round-tripping MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Reimplements the good parts of #200 on top of current development (which already grew an equivalent smoothing_length API), so it lands cleanly instead of a 145-commit reconcile-rebase. For the scalar, vector and multi-component projectors (tensor inherits the scalar): - a dimensional input (Pint Quantity / UnitAware) is non-dimensionalised and reduced to a plain float — uw.non_dimensionalise() returns a dimensionless UWQuantity which sympify() cannot square (the bug flagged on #200); - negative lengths raise ValueError instead of silently storing nonsense; - _smoothing_is_dimensional is tracked so the getter round-trips a dimensional input as a Pint Quantity and a plain input as a plain float. Adds tests/test_0505_projection_smoothing_length.py (tier_a/level_1, 8 tests) locking the round-trip on all four projector classes, the Pint-quantity path, the negative-raises contract, and the screened-Poisson low-pass behaviour. Supersedes #200. Underworld development team with AI support from Claude Code --- src/underworld3/systems/solvers.py | 95 +++++++---- .../test_0505_projection_smoothing_length.py | 159 ++++++++++++++++++ 2 files changed, 217 insertions(+), 37 deletions(-) create mode 100644 tests/test_0505_projection_smoothing_length.py diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 30f93ee0..5a95a53f 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2487,15 +2487,16 @@ def smoothing_length(self): except (TypeError, ValueError): return sympy.sqrt(s) if sval < 0: - return None + raise ValueError( + f"smoothing is negative ({sval}); smoothing_length is undefined") L_nd = sval ** 0.5 - # Re-dimensionalise to length units if a scaling context - # is set; fall back to the plain ND float otherwise. - try: + # Return a Pint Quantity only if the user set a dimensional value via + # the setter; plain-float input round-trips as a plain float so the + # meaning of the number the user passed is preserved. + if getattr(self, "_smoothing_is_dimensional", False): return uw.scaling.dimensionalise( L_nd, uw.scaling.units.meter) - except Exception: - return L_nd + return L_nd @smoothing_length.setter def smoothing_length(self, L): @@ -2507,18 +2508,20 @@ def smoothing_length(self, L): before being squared and stored as ``self._smoothing``. """ self._needs_function_rewire = True - # Unit-aware: route through non_dimensionalise so the - # caller can pass `2.0 * uw.scaling.units.meter` or a - # plain float interchangeably. - try: + # Unit-aware: a dimensional input (Pint Quantity / UnitAware) is + # non-dimensionalised through the active scaling context and reduced to + # a plain float (uw.non_dimensionalise returns a dimensionless + # UWQuantity, which sympify can't square); a plain number is taken as + # already non-dimensional. + is_dim = hasattr(L, "magnitude") or hasattr(L, "units") + if is_dim: L_nd = uw.non_dimensionalise(L) - except Exception: - # Fall back to magnitude-or-float coercion if the - # value doesn't carry/expect units. - if hasattr(L, "magnitude"): - L_nd = L.magnitude - else: - L_nd = L + L_nd = float(getattr(L_nd, "magnitude", L_nd)) + else: + L_nd = float(L) + if L_nd < 0: + raise ValueError(f"smoothing_length must be ≥ 0, got {L_nd}") + self._smoothing_is_dimensional = is_dim self._smoothing = sympify(L_nd) ** 2 @property @@ -2709,25 +2712,34 @@ def smoothing_length(self): except (TypeError, ValueError): return sympy.sqrt(s) if sval < 0: - return None + raise ValueError( + f"smoothing is negative ({sval}); smoothing_length is undefined") L_nd = sval ** 0.5 - try: + # Return a Pint Quantity only if the user set a dimensional value via + # the setter; plain-float input round-trips as a plain float. + if getattr(self, "_smoothing_is_dimensional", False): return uw.scaling.dimensionalise( L_nd, uw.scaling.units.meter) - except Exception: - return L_nd + return L_nd @smoothing_length.setter def smoothing_length(self, L): """Set the smoothing length scale (unit-aware).""" self._needs_function_rewire = True - try: + # Unit-aware: a dimensional input (Pint Quantity / UnitAware) is + # non-dimensionalised through the active scaling context and reduced to + # a plain float (uw.non_dimensionalise returns a dimensionless + # UWQuantity, which sympify can't square); a plain number is taken as + # already non-dimensional. + is_dim = hasattr(L, "magnitude") or hasattr(L, "units") + if is_dim: L_nd = uw.non_dimensionalise(L) - except Exception: - if hasattr(L, "magnitude"): - L_nd = L.magnitude - else: - L_nd = L + L_nd = float(getattr(L_nd, "magnitude", L_nd)) + else: + L_nd = float(L) + if L_nd < 0: + raise ValueError(f"smoothing_length must be ≥ 0, got {L_nd}") + self._smoothing_is_dimensional = is_dim self._smoothing = sympify(L_nd) ** 2 @property @@ -3052,25 +3064,34 @@ def smoothing_length(self): except (TypeError, ValueError): return sympy.sqrt(s) if sval < 0: - return None + raise ValueError( + f"smoothing is negative ({sval}); smoothing_length is undefined") L_nd = sval ** 0.5 - try: + # Return a Pint Quantity only if the user set a dimensional value via + # the setter; plain-float input round-trips as a plain float. + if getattr(self, "_smoothing_is_dimensional", False): return uw.scaling.dimensionalise( L_nd, uw.scaling.units.meter) - except Exception: - return L_nd + return L_nd @smoothing_length.setter def smoothing_length(self, L): """Set the smoothing length scale (unit-aware).""" self._needs_function_rewire = True - try: + # Unit-aware: a dimensional input (Pint Quantity / UnitAware) is + # non-dimensionalised through the active scaling context and reduced to + # a plain float (uw.non_dimensionalise returns a dimensionless + # UWQuantity, which sympify can't square); a plain number is taken as + # already non-dimensional. + is_dim = hasattr(L, "magnitude") or hasattr(L, "units") + if is_dim: L_nd = uw.non_dimensionalise(L) - except Exception: - if hasattr(L, "magnitude"): - L_nd = L.magnitude - else: - L_nd = L + L_nd = float(getattr(L_nd, "magnitude", L_nd)) + else: + L_nd = float(L) + if L_nd < 0: + raise ValueError(f"smoothing_length must be ≥ 0, got {L_nd}") + self._smoothing_is_dimensional = is_dim self._smoothing = sympify(L_nd) ** 2 @property diff --git a/tests/test_0505_projection_smoothing_length.py b/tests/test_0505_projection_smoothing_length.py new file mode 100644 index 00000000..1fdd4b35 --- /dev/null +++ b/tests/test_0505_projection_smoothing_length.py @@ -0,0 +1,159 @@ +"""Locks the `smoothing_length` API on the four Projection classes. + +The projection solvers form + u − ∇·(α ∇u) = ũ, +i.e. a screened-Poisson smoother whose Green's function decays as +exp(−r/L) with L = √α. The L-valued, unit-aware accessor +`smoothing_length` is a convenience view on `smoothing = α`; this test +locks the round-trip and the Pint-unit handling. +""" +import numpy as np +import pytest + +import underworld3 as uw +import underworld3.systems as sys + + +@pytest.fixture(scope="module") +def mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=0.2, qdegree=2, + ) + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_scalar_smoothing_length_roundtrip(mesh): + """L=0.05 ⇒ α=L²=0.0025; reading back gives L (unit-aware).""" + V = uw.discretisation.MeshVariable( + "V_sl", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + proj = sys.Projection(mesh, V) + proj.smoothing_length = 0.05 + assert float(proj.smoothing) == pytest.approx(0.0025) + L = proj.smoothing_length + # In an active scaling context this is a Pint Quantity in metres; + # without one it's a plain float — both should round-trip. + L_num = getattr(L, "magnitude", L) + assert float(L_num) == pytest.approx(0.05) + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_vector_smoothing_length_roundtrip(mesh): + W = uw.discretisation.MeshVariable( + "W_sl", mesh, vtype=uw.VarType.VECTOR, degree=2, continuous=True) + proj = sys.Vector_Projection(mesh, W) + proj.smoothing_length = 0.1 + assert float(proj.smoothing) == pytest.approx(0.01) + L_num = getattr(proj.smoothing_length, "magnitude", proj.smoothing_length) + assert float(L_num) == pytest.approx(0.1) + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_tensor_smoothing_length_roundtrip(mesh): + T = uw.discretisation.MeshVariable( + "T_sl", mesh, vtype=uw.VarType.SYM_TENSOR, degree=2, continuous=True) + Tw = uw.discretisation.MeshVariable( + "Tw_sl", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + proj = sys.Tensor_Projection(mesh, T, Tw) + proj.smoothing_length = 0.2 + assert float(proj.smoothing) == pytest.approx(0.04) + L_num = getattr(proj.smoothing_length, "magnitude", proj.smoothing_length) + assert float(L_num) == pytest.approx(0.2) + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_multicomponent_smoothing_length_roundtrip(mesh): + proj = sys.solvers.SNES_MultiComponent_Projection( + mesh, n_components=3) + proj.smoothing_length = 0.15 + assert float(proj.smoothing) == pytest.approx(0.0225) + L_num = getattr(proj.smoothing_length, "magnitude", proj.smoothing_length) + assert float(L_num) == pytest.approx(0.15) + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_smoothing_setter_keeps_alpha_view(mesh): + """Setting `smoothing` (α) makes `smoothing_length` = √α.""" + V = uw.discretisation.MeshVariable( + "V_smooth", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + proj = sys.Projection(mesh, V) + proj.smoothing = 0.16 + assert float(proj.smoothing) == pytest.approx(0.16) + L_num = getattr(proj.smoothing_length, "magnitude", proj.smoothing_length) + assert float(L_num) == pytest.approx(0.4) # √0.16 + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_smoothing_length_pint_quantity_roundtrip(mesh): + """Pint Quantity input round-trips back as a Pint Quantity in metres. + + Setting `smoothing_length = 0.05 * meter` under an active scaling + context with a 1000 m reference length stores 5e-5 (non-dim) and + α = (5e-5)² = 2.5e-9. Reading back returns a Pint Quantity that + re-dimensionalises to 0.05 m. + """ + model = uw.Model() + model.set_reference_quantities( + length=uw.quantity(1000, "m"), + nondimensional_scaling=True, + ) + uw.use_nondimensional_scaling(True) + try: + V = uw.discretisation.MeshVariable( + "V_pint", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + proj = sys.Projection(mesh, V) + proj.smoothing_length = uw.quantity(0.05, "m") + + # Non-dim store: 0.05 m / 1000 m = 5e-5; α = 2.5e-9. + assert float(proj.smoothing) == pytest.approx(2.5e-9) + + # Pint-quantity round-trip. + L_out = proj.smoothing_length + assert hasattr(L_out, "magnitude"), \ + f"Pint input should round-trip as a Pint Quantity, got {type(L_out)}" + assert float(L_out.to("m").magnitude) == pytest.approx(0.05) + finally: + uw.use_nondimensional_scaling(False) + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_smoothing_length_negative_raises(mesh): + """Negative `smoothing_length` is ill-posed and must raise ValueError.""" + V = uw.discretisation.MeshVariable( + "V_neg", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + proj = sys.Projection(mesh, V) + with pytest.raises(ValueError, match="smoothing_length must be"): + proj.smoothing_length = -0.1 + + +@pytest.mark.tier_a +@pytest.mark.level_1 +def test_smoothing_length_smoothes_a_step(mesh): + """End-to-end: smoothing a step function at scale L attenuates the + short-wavelength content. We check that the projected field's + peak-to-peak range shrinks as L grows, which is the screened-Poisson + low-pass behaviour the docstring promises.""" + import sympy + x, y = mesh.X + V = uw.discretisation.MeshVariable( + "V_filter", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + # Sharp step in x at x = 0.5 + target = sympy.Piecewise((1.0, x < 0.5), (0.0, True)) + ranges = {} + for L in (0.0, 0.05, 0.2): + proj = sys.Projection(mesh, V) + proj.uw_function = target + proj.smoothing_length = L + proj.solve() + ranges[L] = float(V.data.max() - V.data.min()) + # Pure L2 projection should saturate the step (≈ 1.0) + assert ranges[0.0] > 0.95 + # Larger L should produce a smoother (smaller-range) reconstruction + assert ranges[0.2] < ranges[0.05] < ranges[0.0]