From ecc506d21336b8745ae4e1b1eb5259c461ac3289 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 4 Jul 2026 08:58:21 +1000 Subject: [PATCH] fix(maths): CellWiseIntegral evaluates on the mesh DS, not a mis-cloned P1 DM (BF-03) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CellWiseIntegral.evaluate() cloned mesh.dm and attached a fresh single-P1-field discretisation (FE.createDefault + setField + createDS), then handed DMPlexComputeCellwiseIntegralFEM a solution vector packed for mesh.dm's multi-field layout. The layout mismatch made the integral read the wrong DOFs: fn=1.0 on the unit square summed to 2.0 instead of 1.0 (reproduced before the fix; the control Integral gives 1.0). This is audit finding LE-02 / worklist item BF-03 (docs/reviews/2026-07/REMEDIATION-WORKLIST.md). The fix rewrites evaluate() to follow the working Integral pattern in the same file: set the JIT objective on mesh.dm's own DS and integrate on mesh.dm directly, so the compiled objective, the DS field layout, and the packed solution vector all agree. The per-cell results are collected in a plain Vec sized num_local_cells * num_fields (PETSc's implicit cell-major layout when the output vector carries no section) and the field-0 column is returned, one value per rank-local cell. Caveat: setting the objective on the shared mesh DS re-inherits the issue-#171 repeated-call cost-growth behaviour that Integral has. That is accepted here for correctness parity — the clone-DM approach that avoided it is exactly what produced wrong answers (PR #172, reverted in 88807c26). The two xfail markers on tests/test_0501_integrals.py (test_cellwise_integrate_constants, test_cellwise_integrate_meshvar) are removed and the tests now pass, with the global sum reduced across ranks so they also pass under MPI (verified serial and np=2). The serial "level_1 and tier_a" gate is identical before and after the fix (295 passed, 7 skipped, 4 xfailed, 1 xpassed). Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_maths.pyx | 50 ++++++++++++++++---------- tests/test_0501_integrals.py | 31 ++++++++-------- 2 files changed, 46 insertions(+), 35 deletions(-) diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index 4061bf29..effaa951 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -259,14 +259,14 @@ class CellWiseIntegral: super().__init__() @timing.routine_timer_decorator - def evaluate(self) -> float: + def evaluate(self): """ Evaluate the cell-wise integral and return results per cell. Returns ------- ndarray - Array of integral values, one per mesh cell. + Array of integral values, one per (rank-local) mesh cell. Raises ------ @@ -300,27 +300,39 @@ class CellWiseIntegral: cdef Vec cgvec cgvec = a_global - # TODO(BUG): clone+createDefault+createDS gives dmc a single P1 field, - # but cgvec is packed for mesh.dm's multi-field layout — the integral - # reads the wrong DOFs and over-counts by ~2x on the unit square. - # Same bug as PR #172 (reverted in PR #173-followup). Tests in - # tests/test_0501_integrals.py::test_cellwise_integrate_* are xfail - # until this is rewritten to integrate against mesh.dm + getDS() - # directly (the pre-PR-172 Integral pattern). - cdef DM dmc = self.mesh.dm.clone() - cdef FE fec = FE().createDefault(self.mesh.dim, 1, False, -1) - dmc.setField(0, fec) - dmc.createDS() - - cdef DS ds = dmc.getDS() + # Integrate on the mesh DM with its own DS (the `Integral` pattern + # above). The JIT objective is compiled against the mesh's + # multi-field packing, so it must be evaluated with the mesh DS and + # the mesh-packed solution vector. An earlier clone + createDefault + # + createDS variant attached a single P1 field to the cloned DM + # while cgvec remained packed for mesh.dm — the layout mismatch made + # the integral read the wrong DOFs and over-count (~2x for fn=1 on + # the unit square; the same defect that forced the PR #172 revert + # for Integral, see 88807c26). + # + # Note: setting the objective on the shared mesh DS re-inherits the + # issue-#171 repeated-call cost-growth behaviour that Integral has — + # accepted here for correctness parity. + cdef DM dm = self.mesh.dm + cdef DS ds = self.mesh.dm.getDS() CHKERRQ( PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]) ) - cdef Vec rvec = dmc.createGlobalVec() - CHKERRQ( DMPlexComputeCellwiseIntegralFEM(dmc.dm, cgvec.vec, rvec.vec, NULL) ) + # DMPlexComputeCellwiseIntegralFEM writes Nf scalars per cell into a + # flat [cell*Nf + field] layout when the output vector carries no + # DM/section. Only field 0 has an objective; the other field slots + # are left at zero. + num_fields = self.mesh.dm.getNumFields() + cStart, cEnd = self.mesh.dm.getHeightStratum(0) + num_cells = cEnd - cStart + + rvec_py = PETSc.Vec().createMPI((num_cells * num_fields, None), + comm=self.mesh.dm.comm) + cdef Vec rvec = rvec_py + CHKERRQ( DMPlexComputeCellwiseIntegralFEM(dm.dm, cgvec.vec, rvec.vec, NULL) ) self.mesh.dm.restoreGlobalVec(a_global) - results = rvec.array.copy() - rvec.destroy() + results = rvec_py.array.reshape(-1, num_fields)[:, 0].copy() + rvec_py.destroy() return results diff --git a/tests/test_0501_integrals.py b/tests/test_0501_integrals.py index 309e22ca..4c27108b 100644 --- a/tests/test_0501_integrals.py +++ b/tests/test_0501_integrals.py @@ -169,31 +169,26 @@ def test_integrate_swarmvar_deriv_00(): # CellWiseIntegral parity tests — the sister class of Integral that returns a -# per-cell array. These are marked xfail because CellWiseIntegral.evaluate() -# clones mesh.dm and attaches a fresh single-field discretisation -# (FE.createDefault + setField + createDS), but then integrates against -# mesh.dm.getGlobalVec() which is packed for the original multi-field layout. -# The mismatch produces a ~2× over-count on simple cases — the same bug that -# motivated the revert of PR #172 for Integral. -# -# These tests will start passing once CellWiseIntegral.evaluate() is rewritten -# to integrate against mesh.dm + mesh.dm.getDS() directly (matching the -# pre-PR-172 Integral pattern). -@pytest.mark.xfail(reason="CellWiseIntegral has the same field-cloning bug " - "that motivated the PR #172 revert") +# per-cell array. Regression tests for the field-cloning bug (LE-02 / BF-03): +# evaluate() used to clone mesh.dm with a fresh single-P1-field DS while the +# solution vector remained packed for mesh.dm's multi-field layout, reading +# the wrong DOFs and over-counting (~2x for fn=1 on the unit square — the +# same defect that motivated the revert of PR #172 for Integral). def test_cellwise_integrate_constants(): calculator = uw.maths.CellWiseIntegral(mesh, fn=1.0) cell_values = calculator.evaluate() + # evaluate() returns one value per rank-local cell; reduce for the + # global total (no-op in serial). + total = uw.mpi.comm.allreduce(np.sum(cell_values)) + # Sum of per-cell volumes equals the total mesh volume (unit square = 1.0). - assert abs(np.sum(cell_values) - 1.0) < 0.001 + assert abs(total - 1.0) < 0.001 return -@pytest.mark.xfail(reason="CellWiseIntegral has the same field-cloning bug " - "that motivated the PR #172 revert") def test_cellwise_integrate_meshvar(): s_soln.data[:, 0] = np.sin(np.pi * s_soln.coords[:, 0]) @@ -201,7 +196,11 @@ def test_cellwise_integrate_meshvar(): calculator = uw.maths.CellWiseIntegral(mesh, fn=s_soln.sym[0]) cell_values = calculator.evaluate() + # evaluate() returns one value per rank-local cell; reduce for the + # global total (no-op in serial). + total = uw.mpi.comm.allreduce(np.sum(cell_values)) + # Sum matches the global Integral: ∫∫ sin(πx) dx dy = 2/π over the unit square. - assert abs(np.sum(cell_values) - 2 / np.pi) < 0.001 + assert abs(total - 2 / np.pi) < 0.001 return