CBF flux on driven boundaries, kdtree module identity, global velocity_error

src/underworld3/__init__.py

@@ -228,8 +228,14 @@ def view():
 # Unit utilities (top-level convenience for user code)
 from .function.unit_conversion import _extract_value
 
-# KDTree backend is ckdtree (nanoflann); see src/underworld3/kdtree.py
+# KDTree backend is ckdtree (nanoflann). Register the submodule alias in
+# sys.modules so `import underworld3.kdtree` and the `underworld3.kdtree`
+# attribute are the SAME module object. A separate shim file used to shadow
+# this attribute on the first `import underworld3.kdtree` (submodule import
+# rebinds the package attribute), silently dropping the memprobe counters
+# for every later test in the run — issue #316.
 import underworld3.ckdtree as kdtree
+sys.modules["underworld3.kdtree"] = kdtree
 import underworld3.cython
 import underworld3.scaling
 import underworld3.visualisation

src/underworld3/cython/petsc_extras.pxi

@@ -63,6 +63,7 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[])
     PetscErrorCode DMPlexSNESComputeBoundaryFEM( PetscDM, void *, void *)
     PetscErrorCode DMPlexSNESComputeResidualFEM( PetscDM, PetscVec, PetscVec, void *)
+    PetscErrorCode DMPlexInsertBoundaryValues( PetscDM, PetscBool, PetscVec, PetscReal, PetscVec, PetscVec, PetscVec )
     PetscErrorCode DMPlexComputeResidualByKey( PetscDM, PetscFormKey, PetscIS, PetscReal, PetscVec, PetscVec, PetscReal, PetscVec, void *)
     PetscErrorCode DMPlexComputeBdResidualSingle( PetscDM, PetscWeakForm, PetscFormKey, PetscVec, PetscVec, PetscReal, PetscVec )
     # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, void *, void *, void *)

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -2250,6 +2250,7 @@ class SolverBaseClass(uw_object):
             t_nd = self._nondimensional_time(time)
             _time_dm_reaction = self.dm
             UW_DMSetTime(_time_dm_reaction.dm, <PetscReal>t_nd)
+            residual_time = <PetscReal>t_nd
 
         self.mesh.update_lvec()
         self.dm.setAuxiliaryVec(self.mesh.lvec, None)
@@ -2280,6 +2281,15 @@ class SolverBaseClass(uw_object):
             dm = self.dm
             xvec = xlocal
             fvec = flocal
+            # Constrained (Dirichlet) DOFs are ABSENT from the global vector,
+            # so the localToGlobal/globalToLocal round trip above leaves them
+            # ZERO in xlocal. Insert the essential boundary values before
+            # integrating: without this the residual is evaluated against a
+            # state whose boundary values are wrong wherever g != 0, and the
+            # 'reaction' on inhomogeneous Dirichlet boundaries is garbage
+            # (issue #407 — g=0 boundaries were accidentally correct).
+            CHKERRQ(DMPlexInsertBoundaryValues(dm.dm, PETSC_TRUE, xvec.vec,
+                                               residual_time, NULL, NULL, NULL))
             CHKERRQ(DMPlexSNESComputeResidualFEM(dm.dm, xvec.vec, fvec.vec, NULL))
 
             # Return the RAW local residual: each rank has computed its OWNED cells'

src/underworld3/function/analytic.pyx

@@ -220,10 +220,27 @@ class SolCx:
         return out
 
     def velocity_error(self, velocity_var):
-        """Relative L2 error of a velocity MeshVariable against the analytic."""
+        """GLOBAL relative L2 error of a velocity MeshVariable against the
+        analytic solution — the same value on every rank.
+
+        The previous rank-local norms measured each rank's own partition:
+        the rank owning the hardest region (e.g. the SolCx viscosity jump)
+        reported an error 10-20x larger than its neighbours, so tests
+        asserting the value were partition-dependent (issue #370: the np4
+        'tolerance miss' was one rank's local norm, not the solution).
+
+        Shared partition-boundary DOFs contribute once per owning rank to
+        the reductions — a small seam-weighting approximation, acceptable
+        for a benchmark diagnostic.
+        """
         import numpy as _np
+        import underworld3 as _uw
+        from mpi4py import MPI as _MPI
         ua = self.evaluate_velocity(velocity_var.coords)
-        return float(_np.linalg.norm(velocity_var.data - ua) / _np.linalg.norm(ua))
+        diff = _np.asarray(velocity_var.data) - ua
+        e2 = _uw.mpi.comm.allreduce(float((diff * diff).sum()), op=_MPI.SUM)
+        a2 = _uw.mpi.comm.allreduce(float((ua * ua).sum()), op=_MPI.SUM)
+        return float(_np.sqrt(e2 / a2))
 
     def evaluate_stress(self, coords):
         """Exact total (Cauchy) stress at ``coords`` (N×2 array), via the

tests/test_0000_imports.py

@@ -59,3 +59,19 @@ def test_underworld_kdtree_import():
 
 def test_underworld_mpi_import():
     import underworld3.mpi
+
+
+def test_kdtree_module_identity():
+    """#316: `underworld3.kdtree` must be ONE module however it is reached.
+
+    A shim submodule used to shadow the package attribute on first
+    `import underworld3.kdtree`, stripping the memprobe counters for every
+    test that ran afterwards (order-dependent level-1 failures).
+    """
+    import underworld3
+    import underworld3.kdtree
+    import underworld3.ckdtree
+
+    assert underworld3.kdtree is underworld3.ckdtree
+    assert hasattr(underworld3.kdtree, "live_count")
+    assert hasattr(underworld3.kdtree, "total_constructed")

tests/test_1019_boundary_flux.py

@@ -44,6 +44,7 @@ def _heatflux_diagnostics(res=48):
     return np.asarray(flux), q_an, bd_q
 
 
+@pytest.mark.skipif(uw.mpi.size > 1, reason="serial diagnostic: rank-local flux norms are 0/0 on non-owning ranks")
 def test_boundary_flux_scalar_heatflux_serial():
     """Surface heat flux reproduces the analytic flux to high accuracy, and its mean is
     the (analytic) Nusselt number — NOT removed."""
@@ -63,3 +64,66 @@ def test_boundary_flux_scalar_heatflux_serial():
     _f, _a, _b = _heatflux_diagnostics()
     c = np.dot(_f, _a) / (np.linalg.norm(_f) * np.linalg.norm(_a))
     print(f"corr={abs(c):.4f} relL2={np.linalg.norm((_f if c>=0 else -_f)-_a)/np.linalg.norm(_a):.4f}")
+
+
+def test_boundary_flux_inhomogeneous_dirichlet_wall():
+    """#407: flux through a g != 0 Dirichlet wall must be exact.
+
+    The reaction was previously evaluated against a state whose constrained
+    DOFs were zero-filled (the global vector carries no constrained DOFs and
+    the round-trip stripped them), so the g=1 wall read ~-17.7 for a true
+    unit flux while the g=0 wall was accidentally exact.
+    """
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0))
+    T = uw.discretisation.MeshVariable("T407", mesh, 1, degree=2)
+    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()
+
+    # Side walls are natural (zero flux), so every node — end nodes
+    # included — must read the exact unit flux, outward-normal signed.
+    for wall, sign in (("Top", -1.0), ("Bottom", +1.0)):
+        _xs, flux = poisson.boundary_flux(wall)
+        flux = np.asarray(flux)
+        assert np.allclose(flux, sign, atol=1e-3), (
+            f"{wall}: flux range [{flux.min()}, {flux.max()}], expected {sign}"
+        )
+
+
+def test_boundary_flux_corner_semantics_all_walls_driven():
+    """Corner reactions integrate flux through BOTH adjacent driven walls.
+
+    T = x + y (exact in P2) with all four walls Dirichlet: interior wall
+    nodes read the exact unit flux; a corner node reads the SUM of the two
+    walls' contributions (0 on one diagonal, 2 on the other) divided by the
+    queried wall's mass. This pins the method's documented corner
+    semantics — consumers reading pointwise flux at a two-Dirichlet-wall
+    corner must expect the mixture (see issue #407 discussion)."""
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0))
+    T = uw.discretisation.MeshVariable("T407c", mesh, 1, degree=2)
+    x, y = mesh.X
+    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
+    for wall in ("Top", "Bottom", "Left", "Right"):
+        poisson.add_dirichlet_bc(sympy.Matrix([x + y]), wall)
+    poisson.solve()
+
+    xs, flux = poisson.boundary_flux("Top")
+    xs = np.asarray(xs)
+    flux = np.asarray(flux)
+    interior = (xs[:, 0] > 1e-6) & (xs[:, 0] < 1.0 - 1e-6)
+    assert np.allclose(flux[interior], 1.0, atol=1e-3)
+    left_corner = np.isclose(xs[:, 0], 0.0)
+    right_corner = np.isclose(xs[:, 0], 1.0)
+    # The Top corner reaction mixes the adjacent wall's flux: cancellation
+    # 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)