diff --git a/docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py b/docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py index 199828edf..611f9068a 100644 --- a/docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py +++ b/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) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index aa48e8553..96703b514 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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 = 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, boundary_bc.PETScID, &wf, diff --git a/src/underworld3/function/expressions.py b/src/underworld3/function/expressions.py index 2ee20a208..48cfb8ae6 100644 --- a/src/underworld3/function/expressions.py +++ b/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): diff --git a/src/underworld3/utilities/_api_tools.py b/src/underworld3/utilities/_api_tools.py index 75a81d78d..e6673cfdb 100644 --- a/src/underworld3/utilities/_api_tools.py +++ b/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 diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 9be243592..4379c9342 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/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() diff --git a/tests/test_0104_expression_parameter_guards.py b/tests/test_0104_expression_parameter_guards.py new file mode 100644 index 000000000..882333206 --- /dev/null +++ b/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 diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py index a940a31d8..9e0ec9867 100644 --- a/tests/test_1019_boundary_flux.py +++ b/tests/test_1019_boundary_flux.py @@ -168,3 +168,69 @@ def test_boundary_flux_corner_semantics_all_walls_driven(): # at Top-Left (opposite outward fluxes), doubling at Top-Right. assert np.allclose(flux[left_corner], 0.0, atol=1e-3) assert np.allclose(flux[right_corner], 2.0, atol=1e-3) + + +@pytest.mark.parametrize("mass", ("auto", "consistent")) +def test_boundary_flux_p1_trace_2d(mass): + """#413: a 2D P1 trace has no edge-midpoint DOFs — recovery previously died + with a bare KeyError assembling P2 line masses. T = 1 - y is exact in P1, so + every wall node must read the exact unit flux with the P1 line mass.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)) + T = uw.discretisation.MeshVariable("T413", mesh, 1, degree=1) + poisson = uw.systems.Poisson(mesh, u_Field=T) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 0.0 + poisson.add_dirichlet_bc(1.0, "Bottom") + poisson.add_dirichlet_bc(0.0, "Top") + poisson.solve() + + for wall, sign in (("Top", -1.0), ("Bottom", +1.0)): + _xs, flux = poisson.boundary_flux(wall, mass=mass) + flux = np.asarray(flux) + assert np.allclose(flux, sign, atol=1e-3), ( + f"{wall} (mass={mass}): flux range [{flux.min()}, {flux.max()}]" + ) + + +def test_volume_residual_fields_insert_essential_values(): + """#411: compute_volume_residual_fields (Stokes-only diagnostic) missed the + #407 insert — its residual must now match _assemble_volume_reaction (the + fixed, validated core) exactly, including on g != 0 Dirichlet walls where the + un-inserted version returned garbage at constrained rows.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)) + v = uw.discretisation.MeshVariable("v411", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("p411", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + x, y = mesh.X + # v = (x, -y): divergence-free, exact in P2, g != 0 on every wall. + for wall in ("Top", "Bottom", "Left", "Right"): + stokes.add_dirichlet_bc(sympy.Matrix([x, -y]).T, wall) + stokes.petsc_use_pressure_nullspace = True + stokes.solve() + + reference = np.asarray(stokes._assemble_volume_reaction()) + out = stokes.compute_volume_residual_fields() + + # Extract each field's slice of the reference using the same local-section + # walk the method uses, and require exact agreement. + local_section = stokes.dm.getLocalSection() + pStart, pEnd = local_section.getChart() + for name, var in stokes.fields.items(): + field_id = getattr(var, "_solver_field_id", None) + if field_id is None: + field_id = var.field_id + indices = [] + for point in range(pStart, pEnd): + dof = local_section.getFieldDof(point, field_id) + if dof > 0: + offset = local_section.getFieldOffset(point, field_id) + indices.extend(range(offset, offset + dof)) + assert np.allclose(out[name], reference[indices], rtol=1e-12, atol=1e-14), ( + f"{name}: compute_volume_residual_fields diverges from " + f"_assemble_volume_reaction on g != 0 walls — essential values not inserted" + )