Issue quick-fix batch 2: six small independent fixes (#415 #447 #413 #411 #451 #414)

docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py

@@ -120,7 +120,7 @@
 # calls needed.
 #
 # The Model must be set up **before** mesh creation so that the mesh
-# inherits the correct coordinate units (km).
+# coordinates carry the model's reference length (H = 30 km per ND unit).
 
 # %%
 # Reference quantities that define the scaling (Table 1, Spiegelman et al. 2016)
@@ -160,12 +160,15 @@
 
 
 # %% jupyter={"source_hidden": true}
-### Set up the mesh — geometry in km (Table 1, Spiegelman et al. 2016)
+### Set up the mesh — geometry in units of H (Table 1, Spiegelman et al. 2016)
 #
-# The gmsh geometry is defined in non-dimensional units of H (model height),
-# then scaled to km.  Domain is 4H × H = 120 km × 30 km.
+# The gmsh geometry is defined in NON-DIMENSIONAL units of H (model height):
+# UW3 reads raw plex coordinates as already non-dimensional against the model's
+# reference length (H = 30 km), so the 4H × H = 120 km × 30 km domain is
+# [-2, 2] × [-1, 0]. Building the geometry in km here (S = H.magnitude) made
+# the solved model 30x too large in every direction — issue #451.
 
-S = H.magnitude  # 30.0 — scale factor: non-dimensional → km
+S = 1.0  # geometry is non-dimensional; reference length H supplies the km scale
 
 if problem_size <= 1:
     cl_1 = 0.25 * S
@@ -198,7 +201,7 @@
     gmsh.option.setNumber("General.Verbosity", 0)
     gmsh.model.add("Notch")
 
-    # Domain outline (non-dimensional × S → km)
+    # Domain outline in non-dimensional units of H (S = 1: 4H × H → [-2,2] × [-1,0])
     Point1  = gmsh.model.geo.addPoint(-2 * S,      -1 * S,    0, cl_1)
     Point3  = gmsh.model.geo.addPoint(+2 * S,      -1 * S,    0, cl_1)
     Point4  = gmsh.model.geo.addPoint( 2 * S,    -3/4 * S,    0, cl_1)

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -2844,7 +2844,16 @@ class SolverBaseClass(uw_object):
         traces raise explicitly. Reaction and mass assembly are partition-independent.
         For vector fluxes, supply an analytic ``normal`` when strict partition
         independence of the normal projection is required; geometric facet-normal
-        averaging at partition seams has a small pre-existing partition sensitivity."""
+        averaging at partition seams has a small pre-existing partition sensitivity.
+
+        .. warning::
+           On CURVED boundaries, P2 **vertex** values converge only slowly: the P2
+           vertex basis has zero surface mean, so vertex reactions carry only the
+           O(h) facet-geometry error, which the recovery faithfully reconstructs
+           (measured: 93%→55% error under one refinement, while midpoints go
+           2.8%→0.7%). Pointwise consumers on curved boundaries should use
+           **edge-midpoint values** or integral/fitted quantities, never vertex
+           values (issue #414). Flat boundaries are exact up to solver tolerance."""
         from underworld3.utilities.boundary_flux import boundary_flux as _bf
         return _bf(self, boundary, mass=mass, remove_mean=remove_mean, normal=normal)
 
@@ -5947,7 +5956,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         triangles, and the required consistent surface-mass solve for 3D P2 triangles.
         Explicit ``"lumped"`` and ``"consistent"`` choices remain available where
         mathematically valid. Three-dimensional recovery currently supports triangular
-        P1/P2 traces only."""
+        P1/P2 traces only.
+
+        .. warning::
+           On CURVED boundaries, P2 vertex values of :math:`\sigma_{nn}` converge
+           only slowly (the vertex basis has zero surface mean, so vertex reactions
+           carry only the O(h) facet-geometry error); edge-midpoint values are
+           superconvergent. Pointwise consumers on curved boundaries should use
+           midpoint or integral/fitted quantities (issue #414)."""
         if self._rotated_freeslip_info is None:
             raise RuntimeError(
                 "boundary_normal_traction requires a completed rotated-free-slip solve.")
@@ -5967,7 +5983,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length).
         ``mass="auto"`` selects lumped recovery where valid and the consistent
         surface-mass solve for 3D P2 triangles. Requires a prior
-        :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`."""
+        :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`.
+
+        .. warning::
+           On CURVED boundaries (annulus/spherical free surfaces), the P2 VERTEX
+           values written into ``field`` converge only slowly; edge-midpoint values
+           are superconvergent. Downstream pointwise use of curved-boundary
+           topography should rely on midpoint/fitted quantities (issue #414)."""
         if self._rotated_freeslip_info is None:
             raise RuntimeError(
                 "dynamic_topography requires a completed rotated-free-slip solve.")
@@ -8177,6 +8199,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             dm = self.dm
             xvec = xlocal
             fvec = flocal
+            # Constrained (Dirichlet) DOFs are ABSENT from the global vector, so the
+            # round trip above leaves them ZERO in xlocal: insert the essential
+            # values before integrating, exactly as _assemble_volume_reaction does,
+            # or the residual is garbage wherever g != 0 (issues #407/#411 — this
+            # duplicated variant missed the #407 fix).
+            CHKERRQ(DMPlexInsertBoundaryValues(dm.dm, PETSC_TRUE, xvec.vec,
+                                               residual_time, NULL, NULL, NULL))
             if cell_indices is None:
                 CHKERRQ(DMPlexSNESComputeResidualFEM(dm.dm, xvec.vec, fvec.vec, NULL))
             else:
@@ -8270,6 +8299,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         cdef PetscDS ds
         cdef PetscWeakForm wf
         cdef DMLabel c_label
+        cdef PetscReal residual_time = 0.0
 
         self._build(verbose, False, None)
 
@@ -8288,6 +8318,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             t_nd = self._nondimensional_time(time)
             _time_dm_boundary_residual = self.dm
             UW_DMSetTime(_time_dm_boundary_residual.dm, t_nd)
+            residual_time = <PetscReal>t_nd
 
         self.mesh.update_lvec()
         self.dm.setAuxiliaryVec(self.mesh.lvec, None)
@@ -8312,6 +8343,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             dm = self.dm
             xvec = xlocal
             fvec = flocal
+            # Same insert as _assemble_volume_reaction: constrained DOFs are absent
+            # from the global vector, so without this the boundary residual is
+            # evaluated against zeroed essential values wherever g != 0
+            # (issues #407/#411 — this duplicated variant missed the #407 fix).
+            CHKERRQ(DMPlexInsertBoundaryValues(dm.dm, PETSC_TRUE, xvec.vec,
+                                               residual_time, NULL, NULL, NULL))
             CHKERRQ(DMGetDS(dm.dm, &ds))
             CHKERRQ(UW_PetscDSGetBoundaryWeakForm(
                 ds, <PetscInt>boundary_bc.PETScID, &wf,

src/underworld3/function/expressions.py

@@ -1148,10 +1148,13 @@ def is_number(self):
 
     @property
     def is_comparable(self):
-        """Delegate to wrapped expression."""
-        if self._sym is not None and hasattr(self._sym, 'is_comparable'):
-            return self._sym.is_comparable
-        return True
+        """Never comparable: a UWexpression is a symbolic placeholder until it is
+        unwrapped, even when its current contents are numeric. Advertising the
+        contents' comparability made sympy's Max/Min try an immediate numeric
+        comparison and reach for Float internals (``_prec``) that a Symbol subclass
+        does not have (issue #415). ``False`` matches plain ``sympy.Symbol``
+        semantics: Max/Min stay unevaluated and resolve after unwrap/JIT."""
+        return False
 
     @property
     def is_extended_real(self):

src/underworld3/utilities/_api_tools.py

@@ -1,3 +1,26 @@
+def _reaches_container(sym, target, UWexpression, _seen=None):
+    """True if ``target`` is reachable anywhere inside ``sym``, walking THROUGH the
+    stored contents of nested UWexpressions. A cycle hidden one wrapper down still
+    recurses at ``.value``/unwrap time — and the depth-capped unwrap can hand back a
+    silently garbled result (~2^50 scaling) instead of raising (issue #447). The
+    visited set bounds the walk on legitimate (acyclic) expression DAGs and also
+    catches cross-container cycles (a → b → a)."""
+    if _seen is None:
+        _seen = set()
+    if not hasattr(sym, "atoms"):
+        return False
+    for atom in sym.atoms(UWexpression):
+        if atom is target:
+            return True
+        if id(atom) in _seen:
+            continue
+        _seen.add(id(atom))
+        inner = atom._sym
+        if inner is not None and _reaches_container(inner, target, UWexpression, _seen):
+            return True
+    return False
+
+
 class Stateful:
     """
     This is a mixin class for underworld objects that are stateful.
@@ -299,6 +322,23 @@ def __set__(self, obj, value):
         #   - Unit chain: container.units → value.units → value._sym.units (UWQuantity)
         # At unwrap/JIT time the chain is followed automatically.
         if isinstance(value, UWexpression):
+            # Read-back write (p = Params.x; Params.x = p): __get__ returned THIS
+            # container, so storing it would nest the expression inside itself and
+            # every later .value walk recurses to stack death (issue #447). A
+            # write-back of the same expression is a no-op.
+            if value is expr:
+                return
+            # A cycle hidden inside the assigned expression's contents (wrapper
+            # holding Params.x, or a → b → a across containers) recurses just the
+            # same, or worse: the depth-capped unwrap silently returns a ~2^50-
+            # scaled result. Reject loudly (#447).
+            if _reaches_container(value._sym, expr, UWexpression):
+                raise ValueError(
+                    f"Assigning this expression to '{self.public_name}' would create "
+                    f"a reference cycle: it contains '{self.public_name}' (possibly "
+                    f"nested inside another expression). Snapshot the current value "
+                    f"explicitly instead, e.g. use `param.sym` or `float(param.value)`."
+                )
             expr.sym = value  # Store symbolic reference, not inner value
 
         # Special case: Plain UWQuantity (not UWexpression) - has ._value and ._pint_qty
@@ -320,6 +360,31 @@ def __set__(self, obj, value):
             if hasattr(value, "_sympify_"):
                 value = value._sympify_()
 
+            # Composite self-reference (Params.x = Params.x * 2): the container
+            # appearing in its own new contents is the same cycle as above, one
+            # level down. Snapshot semantics is the only meaning that terminates:
+            # substitute the container's CURRENT contents in its place (#447).
+            if hasattr(value, "atoms") and expr in value.atoms(UWexpression):
+                import sympy as _sympy
+                if isinstance(expr._sym, _sympy.MatrixBase):
+                    # xreplace would substitute a Matrix INTO a matrix entry and
+                    # store malformed nested-matrix sympy silently.
+                    raise ValueError(
+                        f"Cannot snapshot Matrix-valued '{self.public_name}' inside "
+                        f"a composite self-assignment; assign a fresh Matrix built "
+                        f"from `{self.public_name}.sym` instead."
+                    )
+                value = value.xreplace({expr: expr._sym})
+            # Cycles one wrapper down (plain-sympy tree holding an expression
+            # whose contents reach this container) survive the snapshot above —
+            # reject them loudly too (#447).
+            if _reaches_container(value, expr, UWexpression):
+                raise ValueError(
+                    f"Assigning this expression to '{self.public_name}' would create "
+                    f"a reference cycle through a nested expression. Snapshot the "
+                    f"current value explicitly (`param.sym` / `float(param.value)`)."
+                )
+
             # Update the expression's .sym
             expr.sym = value
 

src/underworld3/utilities/boundary_flux.py

@@ -349,22 +349,41 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim]
     n = len(keys); Rg = np.zeros(n)
     for k, i in gi.items():
         Rg[i] = R_by[k]
+    # Trace order from the DATA, not an assumption: a P2 trace has a reaction DOF at
+    # every edge point (the midpoint key is in R_by), a P1 trace has vertices only.
+    # Assembling P2 line masses against a P1 trace died with a bare KeyError on the
+    # missing midpoint (issue #413). A mix means the field layout is inconsistent
+    # with the trace — raise rather than guess.
+    mids_present = [km in R_by for (ka, km, kb) in uniq]
+    if mids_present and any(mids_present) != all(mids_present):
+        raise NotImplementedError(
+            "2D boundary-flux recovery found edge-midpoint reactions on only part "
+            "of the boundary; mixed P1/P2 traces are not supported."
+        )
+    trace_order = 2 if (mids_present and mids_present[0]) else 1
     if mass == "lumped":
         mL = np.zeros(n)
         for (ka, km, kb), h in uniq.items():
-            mL[gi[ka]] += h / 6.0; mL[gi[km]] += 2.0 * h / 3.0; mL[gi[kb]] += h / 6.0
+            if trace_order == 2:
+                mL[gi[ka]] += h / 6.0; mL[gi[km]] += 2.0 * h / 3.0; mL[gi[kb]] += h / 6.0
+            else:
+                mL[gi[ka]] += h / 2.0; mL[gi[kb]] += h / 2.0
         sig = Rg / mL
     else:
-        # consistent P2 line mass — a dense (n×n) solve in the number of boundary nodes
+        # consistent line mass — a dense (n×n) solve in the number of boundary nodes
         # (O(n^3)); fine for a 1D boundary (n ~ resolution) but prefer the default lumped
         # (O(n), monotone) for very large boundaries.
         M = np.zeros((n, n))
-        Me = np.array([[4., 2, -1], [2, 16, 2], [-1, 2, 4]])
+        Me2 = np.array([[4., 2, -1], [2, 16, 2], [-1, 2, 4]])
+        Me1 = np.array([[2., 1], [1, 2]])
         for (ka, km, kb), h in uniq.items():
-            tri = [gi[ka], gi[km], gi[kb]]; Mh = (h / 30.0) * Me
-            for ii in range(3):
-                for jj in range(3):
-                    M[tri[ii], tri[jj]] += Mh[ii, jj]
+            if trace_order == 2:
+                nodes = [gi[ka], gi[km], gi[kb]]; Mh = (h / 30.0) * Me2
+            else:
+                nodes = [gi[ka], gi[kb]]; Mh = (h / 6.0) * Me1
+            for ii in range(len(nodes)):
+                for jj in range(len(nodes)):
+                    M[nodes[ii], nodes[jj]] += Mh[ii, jj]
         sig = np.linalg.solve(M, Rg)
     if remove_mean:
         sig = sig - sig.mean()

tests/test_0104_expression_parameter_guards.py

@@ -0,0 +1,90 @@
+"""Guards on UWexpression construction and constitutive Parameter assignment.
+
+1. ``sympy.Max/Min(UWexpression, number)`` must construct symbolically instead of
+   crashing on Float internals (issue #415 — ``is_comparable`` advertised the wrapped
+   contents' comparability, so sympy attempted an immediate numeric comparison).
+2. Reading a constitutive Parameter and assigning it straight back must be a no-op,
+   not a self-nesting that recurses to stack death on the next ``.value`` access
+   (issue #447); the composite form ``P.x = P.x * 2`` must snapshot, not cycle.
+"""
+
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_max_min_of_expression_and_number_construct():
+    """Issue #415: Max/Min(UWexpression, number) stays symbolic, resolves on unwrap."""
+    k = uw.expression(r"k_{415}", 1.0, "constant")
+    clamped = sympy.Max(k, 0.5)              # crashed with AttributeError '_prec'
+    floored = sympy.Min(k, 2.0)
+    assert clamped.has(k) and floored.has(k)
+    # The lazy contract: comparison resolves once the expression is unwrapped.
+    from underworld3.function.expressions import unwrap_expression
+    assert float(unwrap_expression(clamped)) == 1.0    # max(1.0, 0.5)
+    assert float(unwrap_expression(floored)) == 1.0    # min(1.0, 2.0)
+    k.sym = 0.1
+    assert float(unwrap_expression(clamped)) == 0.5    # max(0.1, 0.5) — still live
+
+
+def _viscoplastic_stokes():
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0, 0), maxCoords=(1, 1), cellSize=0.5
+    )
+    v = uw.discretisation.MeshVariable("v447", mesh, 2, degree=2)
+    p = uw.discretisation.MeshVariable("p447", mesh, 1, degree=1, continuous=True)
+    st = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    st.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
+    st.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    st.constitutive_model.Parameters.yield_stress = 10.0
+    return st
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_parameter_readback_is_noop():
+    """Issue #447: p = Params.x; Params.x = p must not nest the container in itself."""
+    st = _viscoplastic_stokes()
+    params = st.constitutive_model.Parameters
+    saved = params.yield_stress
+    before = params.yield_stress.sym
+    params.yield_stress = saved              # recursed to stack death before the fix
+    assert params.yield_stress.sym == before
+    # .value walks the contents — this is where the recursion previously died.
+    float(params.yield_stress.value)
+    # The flux still builds (smooth_max over the parameter was the reported site).
+    st.constitutive_model.flux
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_parameter_composite_self_reference_snapshots():
+    """Issue #447 (composite): Params.x = Params.x * 2 snapshots the current value."""
+    st = _viscoplastic_stokes()
+    params = st.constitutive_model.Parameters
+    base = float(params.yield_stress.value)
+    params.yield_stress = params.yield_stress * 2
+    assert np.isclose(float(params.yield_stress.value), 2 * base)
+    float(params.yield_stress.value)         # no recursion
+
+
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_parameter_wrapper_cycle_raises():
+    """Issue #447 (review round): a self-reference hidden inside a wrapper
+    expression must be rejected loudly — pre-guard it was accepted and later
+    either recursed (flux) or was silently unrolled ~2^50 by the depth-capped
+    unwrap."""
+    st = _viscoplastic_stokes()
+    params = st.constitutive_model.Parameters
+    wrapper = uw.expression(
+        r"w_{447}", params.yield_stress * 2, "wrapper holding the parameter")
+    before = params.yield_stress.sym
+    with pytest.raises(ValueError, match="cycle"):
+        params.yield_stress = wrapper
+    assert params.yield_stress.sym == before      # rejected, not half-applied
+    st.constitutive_model.flux                    # still builds