From 9cc867a7c1646ff1f07f4aea3980d9e6237fb928 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 13:18:33 +1000 Subject: [PATCH 1/3] fix(cbf): insert essential boundary values before the reaction residual (#407) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Constrained (Dirichlet) DOFs are absent from the global vector, so _assemble_volume_reaction's localToGlobal/globalToLocal round trip left them ZERO in the state it integrated. The residual-based reaction was therefore evaluated against a solution with wrong boundary values wherever g != 0: flux through an inhomogeneous Dirichlet wall was garbage (measured -17.7 for a true unit flux), interior residuals near the wall were O(1) despite a converged solve, and g=0 boundaries were accidentally exact — which is why every existing test was blind. DMPlexInsertBoundaryValues now fills the essential values (with the solver's non-dimensional time) before DMPlexSNESComputeResidualFEM. Regression tests: exact unit flux through a g=1 wall including end nodes (natural sides); and the corner-semantics case — T=x+y with all four walls driven pins that a corner reaction integrates BOTH adjacent walls' fluxes (cancellation/doubling), the method's inherent corner behaviour that consumers must expect. The pre-existing 'serial' heat-flux test is skip-guarded under MPI (0/0 on non-owning ranks); the file now passes at np2 and np4. Not affected (verified separately on #407): the rotated-constraint reaction path — it reads the assembled A*u-b, and rotated constraints keep their DOFs in the global system. Closes #407 Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_extras.pxi | 1 + .../cython/petsc_generic_snes_solvers.pyx | 10 +++ src/underworld3/kdtree.py | 13 ---- tests/test_1019_boundary_flux.py | 64 +++++++++++++++++++ 4 files changed, 75 insertions(+), 13 deletions(-) delete mode 100644 src/underworld3/kdtree.py diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index f864e0d8c..99297f9ee 100644 --- a/src/underworld3/cython/petsc_extras.pxi +++ b/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 *) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index e9283e5af..dcb507c7e 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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, t_nd) + residual_time = 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' diff --git a/src/underworld3/kdtree.py b/src/underworld3/kdtree.py deleted file mode 100644 index 88193d42f..000000000 --- a/src/underworld3/kdtree.py +++ /dev/null @@ -1,13 +0,0 @@ -# KDTree module — re-exports the active backend's KDTree. -# -# The backend is ckdtree (nanoflann). The original pykdtree backend was -# removed from the dependency set entirely (2026-07, WA-26/LE-22): its -# OpenMP runtime caused fatal double-init crashes on macOS when loaded -# alongside other OpenMP-using libraries (PETSc, numpy, etc.). See -# commit 16cddf5 "kdtree swap out (temporary)" for the original swap. -# -# This module exists so that `import underworld3.kdtree` doesn't break -# and any code referencing `underworld3.kdtree.KDTree` gets the working -# implementation. - -from underworld3.ckdtree import KDTree # noqa: F401 diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py index 644a8af70..d2bc64c74 100644 --- a/tests/test_1019_boundary_flux.py +++ b/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) From 6217ea1a777ea4f02e7b783f30ccf2415b3be464 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 13:18:33 +1000 Subject: [PATCH 2/3] =?UTF-8?q?fix(kdtree):=20one=20module=20identity=20?= =?UTF-8?q?=E2=80=94=20sys.modules=20alias=20replaces=20the=20shadowing=20?= =?UTF-8?q?shim=20(#316)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The package __init__ bound underworld3.kdtree to the ckdtree module, but a real kdtree.py shim also existed: the first 'import underworld3.kdtree' anywhere loaded the shim and REBOUND the package attribute to it — silently dropping the memprobe counters (live_count/total_constructed) for every test that ran afterwards. test_0000_imports triggers exactly that import, which is why the level-1 suite failed in test_0780_memprobe while the same file passed in isolation (issue #316's order dependence). The alias is now registered in sys.modules, so 'import underworld3.kdtree' and the package attribute are the SAME module object; the shim is deleted and an identity regression test pins it. Closes #316 Underworld development team with AI support from Claude Code --- src/underworld3/__init__.py | 8 +++++++- tests/test_0000_imports.py | 16 ++++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py index c57a3bb3f..e1d989f66 100644 --- a/src/underworld3/__init__.py +++ b/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 diff --git a/tests/test_0000_imports.py b/tests/test_0000_imports.py index ec379dbc8..3032fdfb6 100644 --- a/tests/test_0000_imports.py +++ b/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") From d02a0956e06633192d90ebdd58eb7d8e35c24866 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 24 Jul 2026 13:18:33 +1000 Subject: [PATCH 3/3] fix(analytic): velocity_error is a GLOBAL relative L2 norm (#370) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The np4 'tolerance miss' in test_parallel_custom_fmg_stokes_constrained was never a solver accuracy problem: velocity_error computed RANK-LOCAL norms, so each rank reported its own partition's error — 0.00027 to 0.00562 for the SAME converged solution at np4, with the rank owning the SolCx viscosity jump reporting 20x its neighbours. The np2/np4 'difference' was which rank's local norm happened to trip the assertion. The error is now reduced globally (identical on every rank; shared seam DOFs counted once per owning rank — a small documented weighting approximation for a benchmark diagnostic). The global np4 error is 0.0006, comfortably inside the unchanged 5e-3 tolerance; the full test_1017 file passes at np2 and np4. Closes #370 Underworld development team with AI support from Claude Code --- src/underworld3/function/analytic.pyx | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/src/underworld3/function/analytic.pyx b/src/underworld3/function/analytic.pyx index 18ad20018..93ddb1c35 100644 --- a/src/underworld3/function/analytic.pyx +++ b/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